Statistical methods for assessing drug interactions using observational data

With advances in medicine, many drugs and treatments become available. On the one hand, polydrug use (i.e. using more than one drug at a time) has been used to treat patients with multiple morbid conditions, and polydrug use may cause severe side effects. On the other hand, combination treatments have been successfully developed to treat severe diseases such as cancer and chronic diseases. Observational data, such as electronic health record data, may provide useful information for assessing drug interactions. In this article, we propose using marginal structural models to assess the average treatment effect and causal interaction of two drugs by controlling confounding variables. The causal effect and the interaction of two drugs are assessed using the weighted likelihood approach, with weights being the inverse probability of the treatment assigned. Simulation studies were conducted to examine the performance of the proposed method, which showed that the proposed method was able to estimate the causal parameters consistently. Case studies were conducted to examine the joint effect of metformin and glyburide use on reducing the hospital readmission for type 2 diabetic patients, and to examine the joint effect of antecedent statins and opioids use on the immune and inflammatory biomarkers for COVID-19 hospitalized patients.


Introduction
With advances in medicine, many drugs have become available to treat patients.Drug interactions could be concerns when multiple drugs are applied to the same patients.For example, several drugs are prescribed to a patient by health providers when the patient has multiple morbid conditions.Polydrug use is referred as using more than one drug or treatment during the same time period to treat different underlying morbid conditions [29,38,41].Polydrug use may cause adverse side effects, and has increasingly caused concerns [38].On the other hand, the combination of multiple treatments are developed to treat the same condition, such as cancer [27].Combination therapy may be synergistic, thus more effective than each single drug (i.e.monotherapy) [21].In developing combination therapy for a certain disease, drug interaction often refers to the situation where the effect of two drugs is more or less than the predicted additive effect.When the effect of the combination of two drugs is more than their predicted additive effect, the two drugs are said to be synergistic [20,21].When the effect of the combination of two drugs is less than their predicted additive effect, the two drugs are said to be antagonistic [20,21].However, drug interaction can also refer to adverse effects when polydrug is used.When the outcome of interest is treatment efficacy, the synergistic treatment effect is preferred.In the case that the outcome is the adverse side effect, the synergistic effect would amplify the adverse side effect, and should be avoided.
Although the topic on assessing drug interaction is not new, many statistical models have been developed for assessing drug interaction using in vitro cell culture data [20,21] and in vivo human tumor xenograft data [10].These models have the form of , where Y is the outcome variable, h is the dose-response curve for drug 1, d 1 and d 2 are, respectively, the dose level for drug 1 and drug 2, ρ captures the relative potency of drug 2 versus drug 1, and τ captures the interaction of two drugs.τ > 0 indicates that the combination of two drugs is synergistic, and τ < 0 indicates that the combination of the two drugs is antagonistic.These approaches on assessing drug interactions are in line with isobologram and combination index [20,21,30,44].However, these methods do not control confounding variables, and they are not applicable to observational data from clinical settings.Drug-drug interaction and drug-herb interaction in clinical studies have drawn much attention [3,25,27].However, much of these drug interactions in clinical settings were discovered through clinical experience.As multiple drug use increases, drug interaction has been recognized as an important problem recently, yet the rigorous statistical methods assessing drug interaction based on clinical observational data are lacking.The data depositories from routine clinical practice provide a great opportunity to study the treatment effect of combination therapy or the adverse effect of polydrug use.Observational data from electronic health records or claims data are examples of such data resources which can be used to examine treatment efficacy due to combination therapy or adverse outcomes due to polydrug use.
Recently several data mining algorithms were proposed to identify drug interaction quantitatively using clinical data [36,39,40].A commonly used multivariate regression model has the form of g(E(Y|X, d 1 , where g is a known link function, and X is a vector of covariates.However, this type of model adjusts the confounding variables in regression, and assumes that the treatment effects given X are the same for different values of covariates.The parameters (τ * 1 , τ * 2 , τ * 12 ) capture the average treatment effects (ATEs) and drug interactions only under a certain link function g.The generalized marginal structural models (MSMs) directly model the relationship between the treatment (or combination treatment) and its potential outcome.The generalized MSMs can be written as g(E(Y (d 1 ,d 2 ) )) = τ 0 + τ 1 d 1 + τ 2 d 2 + τ 12 d 1 d 2 , where Y (d 1 ,d 2 ) describes the potential outcome if a patient had received treatment (d 1 , d 2 ) [34].Note that not everyone receives the treatment (d 1 , d 2 ).Under certain assumptions, the inverse probability of treatment weighting (IPTW) method weights each subject with the inverse of the probability of the treatment received, and the resulting weighted sample for patients receiving treatment (d 1 , d 2 ) has a similar distribution to the entire study cohort.
That is, in the weighted sample, the distribution of each covariate across different treatment groups are similar, the parameters in the MSMs thus are obtained from the weighted sample and have causal interpretation [5,14,32].The IPTW weights are often obtained from generalized propensity score (GPS) models [17,42].Literature has indicated that the variances of the causal parameters from IPTW method are optimal if propensity score models include confounding variables (i.e.variables associated with both outcome and treatment selection) and predictor variables (i.e.variables associated with outcome but not treatment selection) [4,6].In Section 2, we propose a generalized MSM to assess ATEs and drug interactions, where the confounding variables and predictor variables are obtained from fitting a multivariate regression model using the elastic net variable selection method [45], and the selected variables are included in the propensity score estimation.The proposed statistical methods in this article can be used to investigate drug interaction on treatment efficacy as well as drug interaction on the adverse event, which depends on the outcome of interest.In Section 3, extensive simulation studies are conducted to examine the performance of the proposed method.In Section 4, one case study is conducted to examine the joint effect of metformin and glyburide in reducing hospital readmission in type 2 diabetic patients, and another case study is conducted to examine the joint effect of statins and opioids on the immune and inflammatory biomarkers in hospitalized COVID-19 patients.Section 5 is devoted to conclusions and discussions.
The major contributions of this article include: (1) we establish a connection between the MSMs (suitable for clinical data) and the drug interaction model initiated for preclinical data; (2) we investigate the impact of variable selection on estimating ATEs and drug interactions; (3) we compare the performance of the proposed methods with traditional multivariate regression models in estimating ATEs and drug interactions; and (4) we illustrate the application of the proposed methods in estimating ATEs and drug interactions based on clinical data sets.
Y denotes the outcome variables from an exponential family, which include Gaussian distribution for a continuous outcome and binomial distribution for a binary outcome.X denotes a vector of the confounding variables which impact both the treatment choice and the outcome variable.One example of a confounding variable is comorbid conditions, which impact both the patient's treatment choices and health outcomes.The outcome variable Y could be continuous (e.g. the reading value for a biomarker) or binary (e.g.death, hospital readmission).In the literature, the additive effect of drug 1 and drug 2 was obtained from two dose-response curves when drug 1 and drug 2 are applied alone [20,21].If the treatment effect of the combination treatment of drug 1 and drug 2 is more than the additive effect, the combination treatment is said to be synergistic; and alternatively, if the treatment effect of the combination treatment of drug 1 and drug 2 is less than the additive effect, the combination treatment is said to be antagonistic.The model in the form of h(d 1 + ρd 2 + τ d 1 d 2 ) has been used to capture drug interaction in in vitro study [20,21].The generalized MSM proposed in this section is on the line with the form of E(Y) = h(d 1 + ρd 2 + τ d 1 d 2 ) but for potential outcome.Potential outcomes can help us understand what the ATEs and drug interaction parameters represent [34].In the case of two drugs, there are four potential outcomes (say, Y (0,0) , Y (1,0) , Y (0, 1) , and Y (1,1) ) for each patient.The potential outcome Y (d 1 ,d 2 ) for a patient would be the outcome when the patient had received treatment combination (d 1 , d 2 ), where d 1 and d 2 take values 0 or 1.In practice, only one potential outcome is observed.The observed outcome, say Y, is the potential outcome corresponding to the treatment the patient receives.Suppose that a patient receives treatment and 0 otherwise.We assume the outcome variable follows an exponential family distribution, and we propose the following generalized MSM to assess ATEs and drug interaction: ( Here g is a known monotonic link function.When the outcome is continuous, we may take g as the identity link function, and model (1) becomes When the outcome is binary, we may take g as the logit link function, and model ( 1) becomes When the outcome is count data, the log-link function may be applied.The generalized MSM (1) and its special cases (2) and ( 3) are distinct from the traditional marginal regression model g(E{Y|d 1 , that (i) the generalized MSMs are structural and they model the relationship between potential outcome and treatment (d 1 , d 2 ) as if everyone had received the treatment (d 1 , d 2 ), thus the parameters in the MSM have a causal interpretation and capture the ATEs and drug interactions; (ii) the multivariate regression model could capture the ATEs and drug interactions if the multivariate regression model has an identity link function; (iii) when the multivariate regression model has a logit link function, the multivariate model cannot be used to capture the ATEs and drug interactions; (iv) the traditional marginal regression model does not capture the ATEs and drug interactions.Heuristic proofs for these statements are provided in the Appendix of this article.The proposed generalized MSM (1) is quite similar to the response surface models for assessing drug interactions proposed by Kong and Lee [20] and Lee et al. [21].However, the models proposed in [20,21] are used to assess drug interactions based on in vitro cell culture data, where there are no confounding variables involved.Here we extend these models to assess drug interactions based on clinical observational data, where confounding variables often exist and thus must be controlled for any valid inferences for a causal relationship between combination treatment and outcome variable.
Let us first consider that the link function g is the identity function with the resulting model (2).The ATE due to drug 1 in the absence of drug 2 is defined as E(Y (1,0) ) − E(Y (0,0) ), which is captured by τ 1 based on model (2).The ATE due to drug 2 is defined as E(Y (0,1) ) − E(Y (0,0) ), which is captured by τ 2 .We can also see that the ATE due to drug 1 in the presence of drug 2 is E(Y (1,1) ) − E(Y (0,1) ), which equals τ 1 + τ 12 .That is, when τ 12 = 0, the ATE due to drug 1 varies depending on whether drug 2 is administered or not.When τ 12 > 0, the ATE due to drug 1 in the presence of drug 2 is larger than that in the absence of drug 2, indicating that drug 1 and drug 2 are synergistic.Vice versa, when τ 12 < 0, the ATE due to drug 1 in the presence of drug 2 is less than that in the absence of drug 2, indicating that drug 1 and drug 2 are antagonistic.This definition is consistent with the definition of drug interaction on the additive scale [39], which states that the sum of ATE of each individual treatment is different from the ATE when the two drugs are applied together.That is, (1,1) ) − E(Y (0,0) ). ( Note that E(Y (1,0) ) − E(Y (0,0) ) = τ 1 and E(Y (0,1) ) − E(Y (0,0) ) = τ 2 , while E(Y (1,1) ) − E(Y (0,0) ) = τ 1 + τ 2 + τ 12 .Equation (4) above is equivalent to τ 12 = 0, thus the parameter τ 12 captures drug interaction on the additive scale.When the outcome is binary, the logistic MSM (3) is often used.We illustrate that the causal parameter τ 12 captures drug interaction in the odds ratio scale [39].That is, τ 12 captures whether the product of causal odds ratios of each individual treatment is different from the odds ratio when two treatments are applied together.That is, OR (1,0) vs (0,0) × OR (0,1) vs (0,0) = OR (1,1) vs (0,0) .
When the risk ratio is the quantity of interest, the following MSM with log-link function could be applied: The interaction on the multiplication scale is defined when the product of causal risk ratios of individual treatment is different from the causal risk ratio when two treatments are applied together [39].That is, Pr(Y (0,0) =1) , which is equivalent to τ 12 = 0. Therefore, τ 12 captures drug interaction in the multiplication scale [39].
In summary, the parameter τ 12 in the generalized MSM (1) captures drug interaction.τ 12 > 0 implies that the effect of the combination treatments is more than expected, indicating a synergistic effect; τ 12 < 0 implies that the effect of the combination treatments is less than expected, indicating an antagonistic effect.

Estimation of the parameters in MSMs
Note that the aforementioned generalized MSMs are models for potential outcome Y (d 1 ,d 2 ) when the subject had received combination treatment (d 1 , d 2 ).However, the potential outcome Y (d 1 ,d 2 ) is not observed if the subject does not receive the combination treatment (d 1 , d 2 ).Although the parameters τ = (τ 0 , τ 1 , τ 2 , τ 12 ) in the generalized MSM (1) have a causal interpretation, an appropriate estimating method must control for confounding variables.Following the causal inference literature, we apply the IPTW method to estimate the causal parameters τ in the generalized MSMs [5,32].The IPTW method essentially creates a weighted sample where the distributions of each confounding variable across different treatment groups are similar, thus removing the confounding effect between the treatment assignment and the outcome variable.The consistency of the estimates for τ provided by the MSMs with IPTW holds under the following assumptions [5,[12][13][14]17,33]: (i) Weak ignorability (i.e., weak unconfoundness): That is, given X, the potential outcome is independent of the treatment received; (ii) Positivity: , that is, the observed outcome is the same as the potential outcome corresponding the treatment received; and (iv) Correctly specification of propensity score model.The first three assumptions are key in the causal inference literature.The fourth assumption is specifically required for the IPTW method [5,19,32].Violation of any one of the assumptions may result in biased estimates for ATEs and drug interactions when our proposed IPTW method is applied.
Let us denote (d i1 , d i2 , x i , y i ) as the observed quadruplets for i th subject (i = 1, 2, . . ., n).The IPTW weight for i th subject is: .
Thus, we form a weighted sample where i th subject has w i copies of (d i1 , d i2 , x i , y i ) instead of 1 copy.Let us denote the distribution of X in the weighted sample as f * , and in the original sample as f.The weight for an observation with covariate X = x and treatment Thus, in the weighted sample, the distribution of X in each treatment group is proportional to the marginal distribution of X in the population where the original sample comes from.That is, X is not associated with treatment assignment anymore in the weighted sample.In practice, the probability of treatment assignment, say Pr( ), needs to be estimated.In the literature, parametric methods (e.g.multinomial regression [17], covariate balance propensity score (CBPS) method [18]) and non-parametric method (e.g.generalized boosting method (GBM) [26]) have been proposed to estimate the GPSs [42].
In this article, we apply the multinomial logistic regression to estimate GPSs, which is easy to implement and fast to compute.To obtain the GPSs, say Pr(D 1 . Here ) T for (d 1 , d 2 ) ∈ D, and β (0,0) takes a vector of zero values.The GPS can be obtained as: ) The parameters β (d 1 ,d 2 ) with (d 1 , d 2 ) ∈ D − (0, 0) can be estimated from the maximum likelihood (ML) method using the observed treatment assignments (d i1 , d i2 ) and confounding variables are estimated, the weight for i th subject is obtained as ŵi = , where the Pr(D 1 = by their ML estimates in Equation (6).
Recent studies indicate that a more efficient estimate for ATEs can be obtained by including only the confounding variables and predictor variables in the GPS estimation [4,6].Note that the confounding variables and predictor variables are associated with the outcome variable.We assume that the multivariate outcome regression model includes treatment D 1 , D 2 , D 1 D 2 , and all covariates.We propose applying the LASSO method (i.e. using L 1 penalty) [45] to the multivariate outcome regression model to select the covariates which are associated with the outcome variable.The selected covariates are included in the GPS model.
We should also assess the balance of the covariates in the application.GPSs are also called balance score [33].A metric to evaluate whether a covariate is balanced is the standardized mean difference (SMD) among all groups.The SMD between two treatment groups, say between treatment (d 1 , d 2 ) and (d 1 , d 2 ), is defined as: , where X(d 1 ,d 2 ) and X(d 1 ,d 2 ) are, respectively, the sample means of the covariate in the is the standard deviation of the covariate based on the observations in the (d 1 , d 2 ) group.The SMDs in the weighted sample are defined similarly but with X(d 1 ,d 2 ) being the weighted sample mean in the (d 1 , d 2 ) group.The summarized metric of covariate balance across all groups for a covariate is obtained by the mean of SMDs over all pairs of different treatment groups.An SMD greater than 0.1 is considered a sign of imbalance of the covariate [43].SMD allows researchers to quantitatively compare balance in measured baseline covariates between two groups in the IPTW weighted sample [2].However, SMD also can be used to assess the balance of confounding variables in the original sample, which helps determine whether or not weighting is needed to correct for confounding variables.Furthermore, SMD can be used to assess the balance of confounding variables under different GPS models.The GPS model should be correctly specified and should result in balanced confounding variables.In case that the confounding variables in the weighted sample are not balanced, more advanced models such as CBPS model [18] and GBM [26] can be applied to estimate the GPSs.The balance of the confounding variables in the weighted sample should be assessed and achieved.
Under the four assumptions (i.e.weak ignorability, positivity, consistency, and correct specification of GPS model), there is no confounding anymore in the weighted sample.The parameters τ = (τ 0 , τ 1 , τ 2 , τ 12 ) in the generalized MSM (1) can be obtained by maximizing the weighted log-likelihood function.That is, Here, l(τ ; d i1 , d i2 , y i ) is the log-likelihood function for i th observation.When the outcome is continuous and the identity link function is used, the log-likelihood function for i th observation is When the outcome is binary and the logit link function is used, the log-likelihood function for i th observation is The seminal work [32] on MSMs indicates that the weighted ML results in consistent estimator for the causal parameter for τ under the four assumptions for causal inference.The weighted ML estimate for τ can be obtained by using the R package survey [23], where the weights are obtained by the GPS model which achieves the balance of confounding variables.Although a robust variance estimator for τ can be obtained from the survey package, it does not incorporate the uncertainty in estimating the GPSs in the IPTW method.Instead, we use the bootstrap sampling techniques to estimate the variance of τ .That is, we obtain B (say, 100) bootstrap samples from the original sample.For the b th bootstrap sample (b = 1, . . ., B), we repeat the same estimating process as outlined to obtain an estimate τ (b) for τ .Var( τ ), the estimate of the variance of τ , is obtained as the variance of the B bootstrap estimates τ (b) (b = 1, . . ., B) [28].In the following section, extensive simulation studies were carried out to examine the performance of the proposed methods in estimating ATEs and drug interactions.

Simulation studies
In this section, we carried out extensive simulation studies to examine the performance of the proposed method in estimating ATEs and drug interactions using the generalized MSMs with IPTW method.

Design of simulation studies
In our simulation studies, we examined the performance of our proposed method under three different sample sizes (say n=500, 1000, and 5000) for continuous responses as well as binary responses (say Y).For continuous responses Y, we considered the following two regression models, the first one assumed homogeneous treatment effects (i.e. the conditional treatment effect given X was the same over different covariates X): and the second one assumed heterogeneous treatment effects (i.e. the conditional treatment effect given X varied over different covariates X): Here we set τ * 1 = τ * 2 = 1, and τ * 12 was taking values from -1 to 1 by a step of 0.5, Here X is a vector of p covariates (p = 10) with each covariate being independently normally distributed with mean zero and variance 1. was a random error with normal distribution of mean 0 and variance 0.5 2 .γ was set as (1, 2, −2, 2, −2, 2, 0 5 ) T so that the outcome models only depended on the first 5 covariates.Here 0 q represented a vector of zero with q components.
For binary responses Y, we considered the following two logistic regression models, the first one assumed homogeneous treatment effects: and the second one assumed heterogeneous treatment effects: Here τ * 1 = τ * 2 = 2, and τ * 12 was taking values from -1 to 1 by a step of 0.5, δ 1 = δ 2 = δ 3 = 1, X were set the same as those for the continuous outcomes, and γ was set as Given X, the treatment assignment (say, (D 1 , D 2 )) was generated using the multinomial distribution with the following probabilities: where 1, 0.5, 0 3 ) T .Thus, in our simulation setting, the first three covariates were confounding variables, which were related to both the treatment assignment and the outcome; X 4 and X 5 were predictor variables, which were only related to the outcome but not the treatment assignment; X 6 and X 7 were instrumental variables, which were only related to the treatment assignment; and all the other three variables X 8 , X 9 and X 10 were spurious variables which were neither related to the treatment assignment nor to the outcome.Based on the generated covariates and treatments, we generated the continuous responses using models (7) and (8), respectively.We also generated the binary responses using models ( 9) and (10), respectively.
For each outcome model (four models in total), each τ * 12 (five specifications), and each sample size (n=500, 1000, and 5000), we generated 1000 simulated data sets.The data generating and estimating procedures were carried out in the following steps.
Step 1: A data set was generated by (i) first generating n observations for the p covariates, say x i (i = 1, . . ., n), which followed MVN(0, I 2 ); (ii) generating treatment assignment (d i1 , d i2 ) based on the multinomial distribution with probability of treatment assignment model (11) and the covariate x i (i = 1, . . ., n); (iii) generating the outcome y i based on the outcome model along with the covariate x i generated in (i) and treatment assignment Step 2: For each simulated data set, the four potential outcomes {y (0,0) i , y (1,0)   i , y (0,1) i , y (1,1)   i } corresponding to the covariate x i were generated using the outcome model (i = 1, . . ., n).For continuous outcome, the true sample ATEs for drug 1 used alone versus control, for drug 2 used alone versus control, and drug interaction were, respectively, obtained as ), and ATE (1,1) vs. (0,1) −ATE (1,0) For binary outcome, the true sample odds ratios and drug interaction in terms of odds ratio were obtained as OR (1,1) vs. (0,1) Step 3: For each dataset, we estimated the parameters τ = ( τ0 , τ1 , τ2 , τ12 ) in the generalized MSM without using the IPTW and with using the IPTW, respectively, and we also estimated τ using multivariate homogeneous regression model as specified in model (7) for continuous outcome and model ( 9) for binary outcome.In the IPTW approach, the following five sets of covariates were used to estimate GPSs and weights: (i) all available covariates, i.e.X 1 ∼ X 10 ; (ii) true confounding variables only, i.e.X 1 ∼ X 3 ; (iii) confounding variables and instrumental variables, i.e.X 1 ∼ X 3 and X 6 ∼ X 7 ; (iv) confounding variables and predictors, i.e.X 1 ∼ X 5 ; and (v) the set of covariates selected based the LASSO method.We calculated SMDs for each covariate in the unweighted dataset and in the weighted datasets using the weights obtained from the five different estimated GPSs.
Step 4: Bootstrap resampling method was used to estimate the variance for each estimated causal parameter, say τ1 , τ2 , and τ12 , resulted from Step 3.
Step 6: The mean of the 1000 estimates (Est), mean squared errors (MSE), the mean of 1000 standard errors (SE) based on the bootstrap method in Step 4, the empirical standard deviation for the 1000 estimates (E.SD) for each causal parameter, and the true coverage rate (CR) of the 95% confidence intervals for the causal parameters (τ 1 , τ 2 , τ 12 ) (i.e. the ground truth approximations from Step 5) were obtained.
Step 7: Repeat Steps 1-6 for each of the 60 combined settings (four outcome models, five specifications for τ * 12 , and three sample sizes.) We also carried out additional simulation studies by letting β and β (1,1) 1 = (1, −1, 0.5, 0 2 , 1, 1.5, 0 3 ) T in the GPS model (11) to simulate the cases where the positivity assumption was violated [24], where all the other conditions specified in Steps 1-7 remained the same.To examine the performance for correlated covariates, we sampled x i from a multivariate normal with mean 0 and variance matrix as , where = (1 − ρ)I p×p + ρJ × J T with J being a vector of p ones.We set ρ = 0.2 for the moderately correlated case and ρ = 0.5 for the strongly correlated case.

Simulation results
We assessed the balance of the covariates in our simulation study.It has been well known that the propensity scores are used to balance covariates among different treatment groups [33].For each simulated data set, we calculated the average of SMDs for each covariate to assess the balance of covariate among the four treatment groups.Figure S1 in the supplement shows the boxplots of the average of SMDs among all pairs of the four treatment groups for each covariate without and with IPTW for the 1000 simulated data sets with sample size n = 5000, where we presented three sets of covariates in the GPS models: (i) all covariates; (ii) confounding variables and predictors; and (iii) the selected variables via LASSO.A dashed horizontal line at the height 0.1 indicated the threshold on whether the covariate was balanced or not.From Figure S1, we can see that (i) the confounding variables (i.e.X 1 ∼ X 3 ) and instrumental variables (i.e.X 6 and X 7 ) were unbalanced in the original samples, and the predictors (i.e.X 4 and X 5 ) and the spurious variables (i.e.X 8 ∼ X 10 ) were less unbalanced in the original sample; (ii) the GPS model included all covariates balanced all covariates; (iii) the GPS model included confounding variables and predictors only balanced these variables; (iv) the GPS model included the selected variables via LASSO balanced confounding variables and predictors.All GPS models resulted in balanced confounding variables but not necessary instrumental variables.To estimate ATEs and drug interaction appropriately, the confounding variables were important to be controlled.It was clear that different GPS models balanced all the confounding variables across the four groups.
The boxplots of the quantities based on the 1000 simulated data with heterogeneous treatment effects for the three sample sizes and different τ * 12 are shown in Figure 1 for continuous outcomes and Figure 2 for binary outcomes, and the boxplots of these quantities based on the 1000 simulated data sets with homogeneous treatment effects are shown in Figure S2 for continuous outcomes and Figure S3 for binary outcomes in the supplement.Based on these simulation results, we conclude that (i) the estimated ATEs and drug interactions based on the proposed weighting method were close to their corresponding ground truth approximations, while the estimated ATEs and drug interactions without IPTW were far from their corresponding ground truth approximations; (ii) as the sample size increased (n=500, 1000, 5000 in the first, second, and third column, respectively, in each one of Figures 1 and 2, S2 and S3), the estimated ATEs and drug interactions based on the proposed methods became more accurate and variations became smaller; and (iii) although the IPTW-based estimates for ATEs and drug interactions using different sets of covariates were close to the ground truth approximations, the variations of the estimated ATEs and drug interactions with GPS including confounding variables and predictors were the smallest one, followed by the GPS model with covariates selected by LASSO.
The summarized metrics, including Est, MSE, SE, E.SD, and the CR for the true causal parameter are presented in Table 1 for continuous outcomes and Table 2 for binary outcomes for τ * 12 = 0 and the sample size n = 1000.The summarized simulation results with sample size n = 500 and n = 5000 are reported in Tables S1 and S2 in the supplement for continuous outcomes, and in Tables S3 and S4 in the supplement for binary outcomes.Based on the results presented in these tables, we concluded that (i) the SEs based on the bootstrap method were close to the E.SD for each causal parameter, and the true CRs were close to the nominal rate of 0.95, indicating that the bootstrap method performed well in estimating the variances of the causal parameters; (ii) the multivariate regression models could provide unbiased estimates for ATEs and drug interactions when the underlying link function was the identity link function, and the multivariate regression models may provide biased estimates for ATEs and drug interactions (in terms of log(OR)) when the underlying link function was the logit link function.We concluded that the proposed method performed well in estimating ATEs and drug interactions in these simulation studies.12 .The third row shows the estimated τ 12 to capture drug interactions.In each block (i.e. for a fixed τ * 12 ), the first boxplot ("Truth") shows the true sample ATEs or drug interaction, the second boxplot ("UW") shows the estimates according to MSMs without IPTW, and the third to the seventh boxplots show the estimates according to MSMs with IPTW being obtained from the following five different sets of covariates: (i) all covariates ("W_ALL"), (ii) confounding variables only ("W_C"), (iii) confounding variables and instrumental variables ("W_CI"), (iv) confounding variables and predictors ("W_CP"), and (v) covariates selected by LASSO ("W_Lasso").
We would like to clarify that the first boxplot in each block in Figures 1 and 2, S2 and S3 was the boxplot of the true sample ATEs and drug interactions based on the potential outcomes obtained from the underlying outcome models in Step 2 in the simulation algorithm.The four potential outcomes for each patient in a simulated data set were obtained from the underlying models with the generated covariates and four different treatment combinations.The ground truth approximations to ATEs and drug interaction reported in Tables 1 and 2 in this article and Tables S1 and S4 in the Supplement were the mean of the 1000 true sample ATEs and drug interactions.The true sample ATEs and drug interactions were not estimates from a sample but were obtained with a known underlying outcome model with plugging in the generated covariates and each one of the four treatment combinations.While the other quantities were obtained from certain estimation procedures as outlined in the simulation study.This may explain why the variations for the true sample ATEs were smaller than the other quantities.12 .The third row shows the estimated τ 12 to capture drug interaction.In each block (i.e. for a fixed τ * 12 ), the first boxplot ("Truth") shows the true sample ATEs or drug interactions, the second boxplot ("UW") shows the estimates according to MSMs without IPTW, and the third to the seventh boxplots show the estimates according to MSMs with IPTW being obtained from the following five different sets of covariates: (i) all covariates ("W_ALL"), (ii) confounding variables only ("W_C"), (iii) confounding variables and instrumental variables ("W_CI"), (iv) confounding variables and predictors ("W_CP"), and (v) covariates selected by LASSO ("W_Lasso").
The positivity assumption for the GPS approach is important, and a violation of the positivity assumption may result in biased estimates for ATEs and drug interactions.The histograms of the GPSs for patients assigned to each group for a sample of size 5000 based on the simulation settings for the GPS model are shown in Figure S4 in the supplement, which clearly indicated that the probability assigned to each group was positive (i.e. the positivity assumption held).The additional simulation studies were carried out by enlarging β in the GPS model (11) to simulate the cases where the positivity assumption was violated [24].The histograms of the GPSs for patients assigned to each group for a sample of size 5000 were shown in Figure S5 in the supplement, which clearly indicated the violation of the positivity assumption.The boxplots for the simulation results for binary heterogeneous outcome model (10) with independent covariates were presented in Figure S6, indicating that the MSMs model based on IPTW method can result in biased estimates for ATEs and drug interactions when the positivity assumption is violated.In the case of moderately and strongly correlated covariates x i (i = 1, . . ., n) drawn from the multivariate normal distribution, the simulation results (not shown) indicated that the correlated covariates may impact the variance of the estimated ATEs and drug interactions, and the biases of the estimated ATEs and drug interactions were largely impacted by the GPS specifications.

Case study 1: glyburide/metformin on readmission rates for diabetes
Recent evidence suggests that the management of hyperglycemia in hospitalized patients has a critical impact on the clinical outcomes for both morbidity and mortality [35].The combination of glyburide and metformin significantly reduced fasting plasma glucose and 2 hour postprandial glucose values compared with either monotherapy, allowing more patients to achieve American Diabetes Association (ADA) treatment goals with lower component doses in drug-naive patients with type 2 diabetes [11].In our case study, we studied the joint effect of glyburide and metformin on readmission rates in individuals admitted to a hospital obtained from electronic health records data.The data set consisted of patients who were hospitalized with diabetic encounter and had recorded diabetes medications.The data set was created from the health fact database which included patients' comprehensive clinical records from 130 US hospitals for years 1999-2008, and the data set was available in the University of California Irvine (UCI) machine learning repository database.The details of the data set can be obtained from Strack et al. [35].We obtained data for a total of 1233 patients who had type 2 diabetes with a HbA1c test result greater than 8% and did not have any other diabetes medications other than glyburide and metformin.Per the ADA guidelines, a HbA1c level 7% is the recommended target, thus the patients in our case study would be considered patients with inadequate glucose control, unless they are older frail adults.Among the patients, 738 patients did not use any diabetes medication (say, control group), 150 patients used glyburide only, 198 patients used metformin only, and 147 patients used both glyburide and metformin (see Table 3(A)).The observed readmission rates were 43.8%, 40.7%, 30.8%, and 32.7%, respectively, for control group, glyburide only group, metformin only group, and the combination of glyburide and metformin group Table 3(A).The covariates included race, gender, age, weight, the comorbidity conditions prior encounter hospital admission such as circulatory diseases, diabetes, digestive diseases, genitourinary diseases, musculoskeletal disease, respiratory diseases, injury, neoplasms, and other diseases.The SMD for some of the covariates in the original sample was larger than 0.1, indicating unbalance of the covariates in the original sample (Figure S7 in the supplement).Thus, the IPTW method, which adjusts confounding variables, should be applied to estimate the ATEs and drug interaction.We applied LASSO to the outcome model, which included all the covariates, treatment indicator variables and the interaction of the two treatments.The selected covariates were race, gender, age, and the comorbidity conditions, which were included in the GPS model.We first applied the multinomial regression model to estimate the GPSs ; however, the balance of the selected covariates among different treatment groups was not achieved.Instead, we applied the CBPS method [18] to estimate the GPS and balance the selected variables.The covariate balance in terms of SMD for this case study in the original observed sample and the IPTW weighted sample are presented in Figure S7 in the supplement.It clearly shows that the SMD for the covariates across the four different groups are small in the weighted sample, indicating balance of the covariates and comorbid conditions.The histograms of the GPSs for the four groups (see Figure S8 in the supplement) indicate that there are not zero or one probabilities, which implies that the positivity assumption is not violated.The MSM with logit link function was used to assess treatment effect and drug interaction.The estimated readmission rates in the weighted samples were 44.0%, 32.2%, 33.2%, and 31.8%,respectively, for control group, glyburide only group, metformin only group, and the combination of glyburide and metformin group (see the column under "Readmission with IPTW (LASSO)" in Table 3(A)).Comparing with the observed readmission rates in Table 3(A), there was a significant difference in readmission rates between the original sample and the weighted sample.For example, the readmission rate for the glyburide only group in the original sample was substantially lower than that in the weighted sample (32.2% vs. 40.7%).The resulting parameters for MSM without IPTW and with IPTW are reported in Table 3(B).Due to the presence of confounding variables, the inference for ATEs and drug interaction should be made based on the estimates applying MSM with IPTW (LASSO).Based on Table 3(B), both glyburide and metformin used alone significantly reduce the readmission rate, with OR as e −0.503 = 0.604 (p=0.019) for glyburide used alone, and OR as e −0.458 = 0.633 (p=0.015) for metformin used alone.Their combination also significantly reduces the readmission rate with an OR=e −0.524 = 0.592 (p=0.012).However, the interaction parameter is not significantly different from zero (p=0.193),indicating that the two drugs do not have interacting effect on hospital readmission based on this case study.
The model proposed can potentially detect a causal relationship between the medications and outcome of re-hospitalization, if confounding variables are controlled.The covariates included in this analysis seek to address variables that make a patient more vulnerable to hospitalization when added glucose-lowering treatments are assigned, such as age and disease burden.This approach helps diminish the likelihood of including rehospitalization for hypoglycemia which is a problematic reason for re-hospitalization in older and frail patients (e.g.older age and higher diagnosis burden).In this example, the medications available in the data base excluded medications other than diabetes medications.Clinically, there are potentially other confounding variables, such as medications that increase glucose levels and impair diabetic control, such as prednisone, which could probably explain the inability to detect a significant interaction effect of metformin and glyburide on re-hospitalization.Another confounder and limitation to the data set is that how the patient actually takes a medication and manages diabetes at home is very difficult to account for.Skipping doses, in addition to poor eating choices or lack of exercise, can negate the impact of medications in any individual.Lastly, although re-hospitalization is the ultimate outcome that most health systems seek to avoid, the breadth of reasons a patient could be re-hospitalized for is vast.In this analysis, all re-hospitalizations were considered the outcome of interest.Re-hospitalizations due to poor glucose control such as urinary frequency, dizziness, fainting, blurry vision and others, were not identifiable.

Case study 2: effect of antecedent statins and opioids use on inflammatory biomarkers in hospitalized COVID-19 patients
Statins are a class of drugs that lower the level of cholesterol in the blood by reducing the production of cholesterol by the liver, and statins also have pleiotropic effects (nonlipid, often beneficial effects).Cholesterol is critical to the normal function of every cell in the body, and it also contributes to the development of atherosclerosis.The impact of statins on coronavirus disease 2019 (COVID-19) severity and recovery is important given their high prevalence of use among individuals at risk for severe COVID-19 [7].Statins inhibit the production of pro-inflammatory cytokines and reactive oxygen species, and diminish platelet reactivity [31].In other words, statins appear to diminish the harms of inflammation and clotting that persons infected with COVID-19 experience, which may lead to lower morbidity and mortality.Some studies indicate that some persons receiving statins during COVID-19 infection have lower morbidity and mortality compared to some not receiving statins [9].Countering the theory that statins diminish COVID-19 infection morbidity and mortality is that statins may up-regulate the ACE2 receptor activity in the lungs, which is how SARS-Co-V2 enters cells, implying that statins may increase infection risk [37].
Opioids are a class of drugs that include the illegal drug heroin, synthetic opioids such as fentanyl, and pain relievers available legally by prescriptions (https://www.drugabuse.gov/drug-topics/opioids).There is increasing recognition that persons engaging in chronic opioid misuse have higher rates of viral and other infections.Opioid use can make people more vulnerable to infection via suppression of immune surveillance.Opioid use is associated with elevated risk of infection due to immune suppression, and paradoxically immune suppression may play a protective role during COVID-19 infection due to mitigated inflammatory response [1].Louisville Kentucky is an epidemic area for opioid crisis.Statins and opoids are therefore theorized to impact COVID-19 outcomes in unconfirmed and variable ways.One approach is to study their impact on biomarkers of inflammation and infection.We used a large database of hospitalized COVID-19 patients established by the University of Louisville Center of Excellence for Research in Infectious Disease (CERID) to study the impact of opioids and statins on the immune and inflammatory biomarkers among patients hospitalized for COVID-19.The data used in this study consisted of 685 adult inpatients hospitalized with COVID-19 infections at nine different hospitals within the Louisville metropolitan area from March 9, 2020 to June 20, 2020.To assess the effect of antecedent statins and opioids use on the immune and inflammatory biomarkers, the patients were formed into four treatment groups according to their antecedent statins and opioids use: control group (i.e.neither statin nor opioids used, n=402), statins only group (n=214), opioids only group (n=34), and statin plus opioids group (n=35).The target immune and inflammatory biomarkers of interest regarding acute COVID-19 infection outcomes, the outcome variables, included laboratory Ct values (i.e. the cycle threshold value in RT-PCR tests for the coronavirus, the smaller the Ct value, the more severe the infection), neutrophils percentage (a special type of white blood cells that are involved in the fight against infection), lymphocytes percentage (white blood cells that are one of the body's main types of immune cells), activated partial thromboplastin time (aPTT) (characterizing coagulation of the blood), and procalcitonin (a higher level of procalcitonin indicates a response to a pro-inflammatory stimulus).Several chronic co-morbidities and factors are associated with elevated pro-inflammatory states and elevated inflammatory biomarkers [8].Race, gender, age, weight, body mass index (BMI), and different comorbidity conditions are possible confounding variables because they are associated with treatment choices as well as outcome risks.
Generalized propensity scores were obtained from multinomial regression models and were employed to balance the covariates in each group.Figure S9 in the supplement presents the SMDs for the original sample and the weighted sample.It shows that the SMDs in the original sample are much larger than the threshold value 0.1, and the SMDs in the weighted sample are smaller or close to the threshold value 0.1, indicating that covariates are balanced among the four groups in the weighted sample.We then applied IPTW to estimate the ATEs and drug interactions on the different outcomes.and the weighted samples (see Table 4(B)).From Table 4, there were some substantial differences between the original observed sample and the weighted sample.For example, the lymphocytes percentages in the original sample versus in the weighted sample were 19.3% versus 18.8% in control group, and 19.6% versus 17.8% in opioids only group.We report the estimated ATE and drug interaction parameters using MSM without IPTW as well as with IPTW (see Table 5).Due to the unbalanced covariates in the original sample, the statistical inference for ATEs and drug interactions should be based on the weighted sample, i.e. the estimates from MSMs with IPTW.The variances for the estimated parameters in the MSMs were obtained from the bootstrap method.Since there were five biomarkers in this investigation, Bonferroni corrections were applied to adjust the p-values for statistical inferences.Based on Table 5, statins significantly increased neutrophils percentage and aPTT, and statins significantly decreased lymphocytes percentage.Opioids used alone did not change the aforementioned biomarkers significantly (although the p-value for the procalcitonin level was 0.046, it became insignificant after Bonferroni corrections).In particular, statins used alone increased aPTT significantly, and opioids used alone did not change aPTT significantly.However, the effect of statins on aPTT in the presence of opioids use was significantly decreased comparing with the effect of statins on aPTT in the absence of opioids by a magnitude of −9.5 unit (adjusted p=0.004×5=0.02).That is, statins and opioids did interact on aPTT significantly.The multiple pathways associated with immune response and sequent inflammatory cascade may be reflected in this analysis.For example, an increase in aPTT in statin use alone may indicate the impact of statin use on higher viral load and intensity of infection, which in turn is associated with higher inflammatory levels.This finding supports the increased ACE2 receptor up-regulation hypothesis with statin use.And, the presence of opioids mitigating aPTT levels in statin users supports the hypothesis that opioids may mitigate inflammatory response.This result could also reflect the unknown impacts of opioids and statins on innate immune function and inflammation.The mechanism of the interaction of statins and opioids may be further investigated, which is beyond the scope of this work.We note that the histogram of GPSs for this case study (Figure S10 in the Supplement) indicated that some propensity scores were close to zero and some were close to one, which implies that the positivity assumption might be violated.We performed sensitivity analysis by trimming weights [22]] to the range of ( 1 0.95 , 1 0.05 ).The parameters in the generalized MSM with the trimmed weights were presented in Table S5 in the Supplement.The results were similar to those in Table 5.However, the violation of the positivity assumption was a concern, and a confirmatory conclusion requires a further study.

Conclusion and discussion
In this article, we propose the generalized MSMs and provide the procedures for estimating ATEs and drug interactions using observational data, where the confounding variables are controlled via the IPTW method.This proposed method, paired with strong clinical modeling and the appropriate data set, provides a novel approach to studying medication signals embedded in complex clinical scenarios.Our extensive simulation studies illustrate that the proposed models and estimating algorithms provide consistent estimates for the causal parameters in the MSMs and capture the true ATEs and drug interactions under the assumptions of weak ignorability, positivity, consistency and correct specification of GPS.However, in practice, there could be many factors that can confound the translation from the medication list to outcomes.If there are unmeasured confounders, the underlying assumptions for the proposed method do not hold, thus the proposed method may not be suitable to assess ATEs and drug interactions.However, the focus of this article is the analytic methodology of drug-drug interaction under these common assumptions.The cases presented illustrate drug-drug interactions that may be detected from observational data for potential prospective clinical study.The two case studies illustrate the method described in this paper, with acknowledgement that the case studies are not based on comprehensive clinical models aimed at controlling all confounding variables.The case studies illustrate the MSM modeling and identify drug-drug interactions associated with outcomes of consequence, providing a potential starting point for clinical studies.Presentation of these case studies illustrates the potential utility of the model in mining observational retrospective data to identify potential drug-drug synergies for future clinical studies.
Once the generalized propensity scores are obtained, one can use the Horvitz-Thomp son survey sampling estimator [16] to estimate the mean of the potential outcomes, which can be used to further estimate ATEs and drug interactions.Never the less, the Horvitz-Thompson method neither directly estimates ATEs and drug interactions nor estimates their variances.On the contrary, the generalized MSMs based on the weighted sample directly estimate ATEs and drug interactions and their variances.Thus, the generalized MSMs along with IPTW provide a rigorous statistical method for assessing drug interactions.With robust interdisciplinary collaboration and clinical modeling, the impact of drug interactions detected using this method may identify significant leads for a biomedical clinical study.
This article presents the MSMs and algorithms for estimating ATEs and drug interactions when two drugs are used together and each drug has only two levels (present or not).However, the doses of the medicine sometimes vary among patients and may be adjusted based on patients' conditions, where the treatments could be continuous variables.The proposed approach could be extended to the situation when each drug has multiple levels or in continuous scale using the GPS [15], which will be investigated in our future research.

Figure 1 .
Figure 1.Boxplots of 1000 estimated ATEs and drug interactions for continuous outcome with heterogeneous treatment effects and different sample sizes.The first row and the second row, respectively, show the estimated ATEs for drug 1 (i.e.τ1 ) and drug 2 (i.e.τ2 ), with different specifications of τ *12 .The third row shows the estimated τ 12 to capture drug interactions.In each block (i.e. for a fixed τ * 12 ), the first boxplot ("Truth") shows the true sample ATEs or drug interaction, the second boxplot ("UW") shows the estimates according to MSMs without IPTW, and the third to the seventh boxplots show the estimates according to MSMs with IPTW being obtained from the following five different sets of covariates: (i) all covariates ("W_ALL"), (ii) confounding variables only ("W_C"), (iii) confounding variables and instrumental variables ("W_CI"), (iv) confounding variables and predictors ("W_CP"), and (v) covariates selected by LASSO ("W_Lasso").

Figure 2 .
Figure 2. Boxplots of 1000 estimated ATEs and drug interactions for binary outcome with heterogeneous treatment effects and different sample sizes.The first row and the second row, respectively, show the estimated ATEs (OR in log scale) for drug 1 (i.e.τ1 ) and drug 2 (i.e.τ2 ), with different specification of τ *12 .The third row shows the estimated τ 12 to capture drug interaction.In each block (i.e. for a fixed τ * 12 ), the first boxplot ("Truth") shows the true sample ATEs or drug interactions, the second boxplot ("UW") shows the estimates according to MSMs without IPTW, and the third to the seventh boxplots show the estimates according to MSMs with IPTW being obtained from the following five different sets of covariates: (i) all covariates ("W_ALL"), (ii) confounding variables only ("W_C"), (iii) confounding variables and instrumental variables ("W_CI"), (iv) confounding variables and predictors ("W_CP"), and (v) covariates selected by LASSO ("W_Lasso").
A: Sample size and estimated readmission rate Table 4 presents the group means of the outcome variables based on the original observed sample (see Table 4(A))

Table 4 .
Summarized statistics for each of the five outcome variables (mean ± standard deviation) across different treatment groups: control, statins use only, opioids use only, and opioid and statins use group for hospitalized COVID-19 patients.: Summarized mean and standard deviation based on the observed sample A

Table 5 .
The estimated ATEs and drug interactions for opioids and statins use.Statins vs Control; τ 2 , Opioids vs Control; τ 12 , Drug interaction; τ 1 +τ 2 +τ 12 , Statins+Opioids.vs Control.[ * ]− The p-values were obtained from the Wald tests without Bonferroni correction.The Bonferroni corrected p-values would be the minimal value of 1 and 5 times of the p-values reported.