REVIEW published: 18 September 2019 doi: 10.3389/fphar.2019.00973 Propensity Score Methods in Health Technology Assessment: Principles, Extended Applications, and Recent Advances M Sanni Ali 1,2,3*, Daniel Prieto-Alhambra 2,4, Luciane Cruz Lopes 5, Dandara Ramos 3, Nivea Bispo 3, Maria Y. Ichihara 3,6, Julia M. Pescarini 3, Elizabeth Williamson 1, Rosemeire L. Fiaccone 3,6,7, Mauricio L. Barreto 3,6 and Liam Smeeth 1,3 1 Faculty of Epidemiology and Population Health, London School of Hygiene and Tropical Medicine, London, United Kingdom, 2 Nuffield Department of Orthopaedics, Rheumatology and Musculoskeletal Sciences (NDORMS), Center for Statistics in Medicine (CSM), University of Oxford, Oxford, United Kingdom, 3 Centre for Data and Knowledge Integration for Health (CIDACS), Instituto Gonçalo Muniz, Fundação Osvaldo Cruz, Salvador, Brazil, 4 GREMPAL Research Group (Idiap Jordi Gol) and Musculoskeletal Research Unit (Fundació IMIM-Parc Salut Mar), Universitat Autònoma de Barcelona, Edited by: Barcelona, Spain, 5 University of Sorocaba–UNISO, Sorocaba, São Paulo, Brazil, 6 Institute of Public Health, Federal Ileana Mardare, University of Bahia (UFBA), Salvador, Brazil, 7 Department of Statistics, Federal University of Bahia (UFBA), Carol Davila University of Medicine Salvador, Brazil and Pharmacy, Romania Reviewed by: Bogdan Ileanu, Randomized clinical trials (RCT) are accepted as the gold-standard approaches to Bucharest Academy of Economic measure effects of intervention or treatment on outcomes. They are also the designs of Studies, Romania choice for health technology assessment (HTA). Randomization ensures comparability, in Douglas Faries, Eli Lilly, United States both measured and unmeasured pretreatment characteristics, of individuals assigned to *Correspondence: treatment and control or comparator. However, even adequately powered RCTs are not M Sanni Ali always feasible for several reasons such as cost, time, practical and ethical constraints, Sanni.ali@lshtm.ac.uk; sanni.ali@ndorms.ox.ac.uk and limited generalizability. RCTs rely on data collected on selected, homogeneous population under highly controlled conditions; hence, they provide evidence on efficacy of Specialty section: interventions rather than on effectiveness. Alternatively, observational studies can provide This article was submitted to evidence on the relative effectiveness or safety of a health technology compared to one Pharmaceutical Medicine and Outcomes Research, or more alternatives when provided under the setting of routine health care practice. In a section of the journal observational studies, however, treatment assignment is a non-random process based on Frontiers in Pharmacology an individual’s baseline characteristics; hence, treatment groups may not be comparable Received: 18 April 2019 Accepted: 31 July 2019 in their pretreatment characteristics. As a result, direct comparison of outcomes between Published: 18 September 2019 treatment groups might lead to biased estimate of the treatment effect. Propensity score Citation: approaches have been used to achieve balance or comparability of treatment groups in Ali MS, Prieto-Alhambra D, terms of their measured pretreatment covariates thereby controlling for confounding bias Lopes LC, Ramos D, Bispo N, Ichihara MY, Pescarini JM, in estimating treatment effects. Despite the popularity of propensity scores methods and Williamson E, Fiaccone RL, recent important methodological advances, misunderstandings on their applications and Barreto ML and Smeeth L (2019) Propensity Score Methods in Health limitations are all too common. In this article, we present a review of the propensity scores Technology Assessment: Principles, methods, extended applications, recent advances, and their strengths and limitations. Extended Applications, and Recent Advances. Front. Pharmacol. 10:973. Keywords: bias, confounding, effectiveness, health technology assessment, propensity score, safety, secondary doi: 10.3389/fphar.2019.00973 data, observational study Frontiers in Pharmacology | www.frontiersin.org 1 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA INTRODUCTION using routinely collected data (Schneeweiss et al., 2011, Ali et al., 2016, Bärnighausen et al., 2017). In observational Randomized clinical trials (RCTs) are generally accepted as the studies, however, treatment selection is mainly influenced by gold-standard approaches for measuring the “causal” effects of the patient, the physician, and, to a certain extent, the health treatments on outcomes (Sibbald and Roland, 1998; Concato system characteristics. Hence, treated and untreated groups et al., 2000) and the design of choice for health technology differ not only in receiving the treatment but also in other assessment (HTA). In causal inference terminology using pretreatment characteristics, leading to non-comparability or Rubin’s potential outcomes framework (Rubin, 2005), the effect non-exchangeability, a phenomenon leading to confounding bias of a certain treatment (Z = 1) versus a control or comparator (Greenland and Morgenstern, 2001). This means that differences (Z = 0) on an outcome (Y) involves comparison of potential in outcomes between the two groups, treated versus untreated, outcomes under treatment (Y1)) and an alternative treatment could be explained by either the treatment, or other pretreatment (Y0)). In RCT, with sufficient numbers of participants and variables, or both. In other words, direct comparison of outcome adequate concealment of allocation, randomization ensures that events between the two groups leads to biased estimate of the individuals assigned to treatment and control or comparator treatment effect. Hence, any systematic difference in pretreatment groups are comparable in all pretreatment characteristics, both characteristics between treatment should be accounted for by measured and unmeasured (Sibbald and Roland, 1998). The only design, or analysis, or both (Rubin, 1997). Over the years, several difference is that one group received the treatment (Z = 1) and methodologies have been developed to control for confounding the other received no treatment or the alternative treatment (Z = bias in observational studies (Figure 1); the propensity score 0); hence, any difference in outcomes between the two groups methods (Rosenbaum and Rubin, 1983) are among the popular can be attributable to the effect of the treatment. In other words, approaches in pharmacoepidemiology and health technology the “causal” effect of treatment in the study population (the evaluations (Ali et al., 2015). average treatment effect, ATE) on outcomes can be estimated Propensity score approaches were first introduced by by a direct comparison of the outcomes between the treatment Rosenbaum and Rubin in 1983 (Rosenbaum and Rubin, 1983), and the comparator groups (Equation 1) (Concato et al., 2000). and their use to control for confounding has been increasing in However, even adequately powered RCT may not always be the previous decade. Propensity score (PS) is a scalar summary of feasible for reasons such as cost, time, ethical, and practical all measured pretreatment characteristics (often called potential constraints (Sibbald and Roland, 1998). RCTs also rely on data confounders); stated formally, the propensity score e(X) is the collected on selected, homogeneous population under highly conditional probability of receiving a certain treatment, versus controlled conditions; hence, they provide evidence on efficacy a comparator or no treatment, given the measured pretreatment rather than on effectiveness of interventions or treatments characteristics (Rosenbaum and Rubin, 1983), X, denoted as (Eichler et al., 2011). e( X ) = pr ( Z = 1|X ), (2) ATE = E[Y1 − Y0 ] = E[Y1 ] − E[Y0 ] (1) where Z = 1 for individuals in the treatment group and Z = 0 With steadily increasing costs of health care and the for individuals in the comparison group (Rosenbaum and Rubin, introduction of novel, yet very expensive, pharmaceutical 1983; Rosenbaum and Rubin, 1984). Treated and untreated products and diagnostics, HTA agencies such as the UK National individuals with similar propensity scores have, on average, Institute for Health and Care Excellence (NCIE) are inquiring similar or comparable pretreatment characteristics, a situation robust methods for evaluation of relative effectiveness and safety similar to an RCT. However, this comparability, conditional on of medications, devices, and diagnostics in daily clinical practice. the propensity score, of the treatment groups is limited only to In contrast to efficacy, relative effectiveness of an intervention measured pretreatment characteristics included in the propensity or treatment is “the extent to which an intervention does more score model and may not hold for unmeasured ones (Rosenbaum good than harm, when compared to one or more alternative and Rubin, 1983). Hence, balancing these pretreatment potential intervention(s)” when used under the routine setting of health confounders through propensity scores enables researchers to care practice” (Eichler et al., 2011; Schneeweiss et al., 2011). obtain a “quasi-randomization” of treatment groups to reduce In addition, for medical devices and diagnostics, waiting for confounding and hence to get a better estimate of the treatment evidence from RCTs when the health technology is diffusing effect. Implicitly, researchers assume “Strongly Ignorable in the clinical practice could be costly for the payers, inefficient Treatment Assignment” (SITA) given the measured covariates; from policy perspective, and methodologically questionable this comprises “unconfoundedness” and “positivity” (Rosenbaum (Tarricone et al., 2016). On the other hand, regulators’ and HTA and Rubin, 1983). Unconfoundedness implies that all relevant agencies’ perception of the importance of real-world data in pretreatment characteristics are measured and included in the complementing evidence on the relative effectiveness of health propensity score model; hence, given these measured covariates technologies has been steadily increasing (Makady et al., 2017; are included in the propensity score, there is no unmeasured Yuan et al., 2018). confounding. Positivity, on the other hand, implies that every The effect of a particular health technology, e.g., a medication, individual has a non-zero (positive) probability of receiving all on a certain outcome event could also be investigated using non- values of the treatment variable: 0 < P(Z = 1|X) < 1 for all values randomized studies (i.e., observational or quasi-experimental) of Z (Rosenbaum and Rubin, 1983). Frontiers in Pharmacology | www.frontiersin.org 2 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA FIGURE 1 | Methods to control for confounding in observational studies.*Multiple imputation is valid when the assumption of Missing at Random (MAR) holds;**if time-varying confounder is affected by previous treatment, all PS-based methods except marginal structural model (MSM) using inverse probability of treatment weight (IPTW) will give biased estimate;***self-controlled case-series design; ♠)stratification using effect modifier and adjustment within the strata to account for other covariates; ♠♠)Disease risk score (prognostic score) method; ♠♠♠)restriction or choosing an active comparison group vs non-user group; ♣)G-formula and ♣♣)G-estimation of structural nested models, which rely on specification of the outcome model; ♣♣♣)instrumental variable methods. (Adapted in part from Schneeweiss (2006), Uddin et al. (2016), and Zhang et al. (2018)). In the last decade, the propensity score methods have approaches; the section Covariate Balance Assessment describes been popular among clinical researchers, their use in methods for assessment of covariate balance in propensity score pharmacoepidemiology and HTAs has been ubiquitous, and methods; the section Propensity Score Methods summarizes the they have undergone substantial methodological advances. On different types of propensity score methods; the section Extended the other hand, confusions and misunderstandings on what a Applications describes extended applications of propensity propensity score method can and cannot do as well as errors in scores; the section Advantages and Limitations of Propensity Score the design, analysis, interpretation, and reporting of propensity Methods summarizes strengths and limitations of the propensity score-based analyses are unfortunately all too common (Ali score methodology; the section Reporting highlights on reporting et al., 2015). With increasing availability of routinely collected of propensity score based analysis; and the section Conclusion electronic medical records for evaluation of effects (both concludes the discussion. comparative effectiveness and safety) of health technologies, and relatively rapid development of the methods, an up-to-date review of the methods and their characteristics is necessary. In VARIABLE SELECTION AND PROPENSITY this article, we aim to introduce propensity score methods with SCORE ESTIMATION an emphasis on important aspects of the methods; describe their extended applications and recent developments; and discuss Observational studies using administrative or clinical databases their strengths and limitations. often involve high dimensionality with respect to the number The manuscript, including the introduction, is organized of pretreatment covariates available for analysis including into eight sections: the section Introduction has introduced RCT, socioeconomic characteristics, demographics, comorbidities, observational studies, and propensity score in relation to HTA; comedications, and health system characteristics, among the section Variable Selection and Propensity Score Estimation others. The inclusion of a large number of covariates in discusses variable selection and propensity score estimation conventional regression models, particularly in nonlinear Frontiers in Pharmacology | www.frontiersin.org 3 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA models such as logistic regression and Cox regression models, confounders and the outcome) can lead to increased bias in the requires sufficient number of outcome events (approximately estimated treatment effect (Pearl, 2011; Myers et al., 2011a, Myers 10 outcome events per covariate) (Peduzzi et al., 1995; et al., 2011b; Pearl, 2012; Ali et al., 2016). Peduzzi et al., 1996; Cepeda et al., 2003). For example, to It is important to emphasize that, similar to conventional adjust for 5 confounders using logistic regression model, regression modelling, intermediates (variables on the causal one would need to have 5*10 = 50 outcome events. However, pathway between the treatment and the outcome) and colliders many practical settings in pharmacoepidemiology and other should not be included in the propensity score model (Greenland HTAs involve relatively few or rare outcome events; hence, and Morgenstern, 2001) since including these variables will confounding adjustment using regression methods requires tend to increase (rather than reduce) bias. In addition, strong selection of a limited number of covariates to avoid problems instruments should also be excluded, particularly when strong such as over-fitting (Peduzzi et al., 1995). Alternatively, the unmeasured confounding is a concern thereby avoiding any use of propensity score methods to summarize a large pool amplification of the residual bias (Pearl, 2011; Myers et al., 2011a; of covariates into a single score, the propensity score, avoids Myers et al., 2011b; Pearl, 2012; Ali et al., 2016). However, it is not over-fitting and collinearity issues in estimating treatment common to come across with such a scenario; the use of propensity effects (Cepeda et al., 2003). When the number of covariates score method is meaningful when the assumption of “Strongly available in the study dataset is relatively small, it is common Ignorable Treatment Assignment”, SITA, is met (i.e., there is no practice to include all the pretreatment covariates in the unmeasured confounding given the measured covariates and also propensity score model; however, covariate selection might there is positivity) (Rosenbaum and Rubin, 1983). Compared be required when researchers are presented with very large to residual confounding by unmeasured characteristics, bias number (several hundreds) of covariates and limited number amplification should be considered a secondary concern; hence, of outcome events (Schneeweiss et al., 2009). researchers should be cautious and are advised to err on the side Covariates selection in propensity score is often based on of including rather than excluding any potential confounder prior subject-matter knowledge on the relationships underlying (Myers et al., 2011b; Ali et al., 2017c). Alternatively, when a the covariates in the study data, statistical tests on the association strong instrument—essentially a proxy measure of difference in between the covariates and the outcome event (using p-values treatment—is identified that is independent of confounders and or change in effect estimates) (Brookhart et al., 2006; Patrick outcome, instrumental variable analysis can be a powerful tool to et al., 2011; Ali et al., 2015; Adelson et al., 2017), strength of account for any unmeasured confounding (Angrist et al., 1996). associations with the treatment and/or the outcome event A common question asked by clinical researchers who have not (Patrick et al., 2011; Ali et al., 2015; Adelson et al., 2017), and used propensity score methods is “why do we need to estimate the machine learning methods such as generalized boosted models probability that an individual receives a certain treatment versus (McCaffrey et al., 2004). Each approach has its own strengths a comparator while we certainly know from the data whether that and limitations; however, emphasis should be given to achieve particular individual has received the treatment?” A brief answer balance on important prognostic pretreatment characteristics to this important question is as follows: propensity score exists (Rosenbaum and Rubin, 1983) and not to improve model fit both in RCT and in observational studies (Joffe and Rosenbaum, or to predict treatment as well as possible. Hence, the use of 1999; Rubin, 2004b; Ali et al., 2016). In RCT, the true propensity p-values, goodness-of-fit tests, and model discrimination tests score is known by design or the treatment allocation mechanism, such as c-statistics should be avoided (Weitzen et al., 2005; i.e., randomization. For example, consider a simple two-arm Patrick et al., 2011; Westreich et al., 2011). The iterative approach RCT in which individuals are assigned to a treatment versus a of model fitting, by including interactions and square terms of comparison group by flipping of a fair coin (also assume that the the covariates, and subsequent balance assessment, which was sample sizes are equal in both treatment groups). The propensity recommended in the seminal paper by Rosenbaum and Rubin score for every individual, the probability of being assigned to (1983), is still a more robust approach. This application helps the treatment group versus the comparator group, is equal to to achieve the goal of propensity score modelling, “improving 0.5, apart from chance variations. In contrast, in observational balance” of potential confounders between treatment groups so studies, the true propensity score for individuals is unknown that the groups are comparable or exchangeable conditional on and is dependent on several pretreatment characteristics, both the propensity score. clinical and nonclinical, under consideration by the physician. One of the greatest strengths of propensity score approaches is As a result, the propensity score should be—and can often be— the separation of design from analysis, i.e., propensity score methods estimated using the study data (Joffe and Rosenbaum, 1999; purposefully disregard outcome information at this stage of the Rubin, 2004b; D’Agostino, 2007; Ali et al., 2016). Estimation of design (Rubin, 2004b; Leacy and Stuart, 2014). That would also the propensity score is needed to create a “quasi-randomized mean, as in the classical implementation of the methods, association experiment” by using the individual’s probability of receiving between the covariates and the outcome event in the study data is not the treatment as a summary score of all measured pretreatment assessed for selection of covariates while constructing the propensity covariates. It enables appropriate adjustment for measured score model. However, this approach is not without disadvantages: potential confounders to estimate the effect of the treatment. This failure to exclude colliders (variables that are common effects of explains one of the key properties of the propensity score method: the treatment and the outcome) and strong instruments (variables if we find two individuals with the same propensity score, one in that are strongly related to treatment but independent of both the the treated group and one in the untreated group, we can assume Frontiers in Pharmacology | www.frontiersin.org 4 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA that these two individuals are more or less “randomly assigned” It is helpful to start propensity score analysis by examining to one of the treatment groups in the sense of being equally likely the distribution of propensity scores using histograms or density to be treated or not, with respect to measured pretreatment plots. This facilitates subjective judgment on whether there is characteristics (Ali et al., 2015; Ali et al., 2016). sufficient overlap, also called “the common support,” between In practice, the propensity score is often estimated using propensity score distributions of treated and untreated groups ordinary logistic regression model, in which treatment status (Dehejia and Wahba, 2002). However, such plots should not is regressed on measured pretreatment characteristics (Austin, be considered as proper measures of covariate balance; they 2008a; Ali et al., 2015). The estimated propensity score is the can guide the choice of matching algorithms in propensity predicted probability of receiving the treatment derived from the score matching and the number of strata in propensity score fitted logistic regression model. Logistic regression has several stratification (Ali et al., 2015; Ali et al., 2016). For example, when advantages: it is a familiar and well-understood statistical tool there is very little overlap in the propensity score distributions, for researchers as well as easy to implement using standard matching treated and untreated individuals with replacement, statistical software packages (Setoguchi et al., 2008; Westreich with or without caliper, can be a better option because it will be et al., 2010; Ali et al., 2016). However, logistic regression is not challenging to find sufficient number of untreated individuals the only approach; other methods have also been used including for the treated individuals (Ali et al., 2016). Inadequate overlap recursive partitioning (D’Agostino, 2007) and several machine in the propensity score distributions, which can be quantified learning methods, for example, classification and regression using overlapping coefficient (Ali et al., 2014), should also warn trees (CARTs), neural networks, and random forests (Setoguchi researchers that the dataset, no matter how large, could not et al., 2008; Lee et al., 2010; Westreich et al., 2010; Lee et al., support any causal conclusion about the effect of the treatment 2011). Comparative simulation studies favor the use of machine on the outcome of interest without relyng o untrustworthy model learning methods over logistic regression when there is moderate assumptions (Rubin, 1997; Ali et al., 2016). or high nonlinearity (square or cubic terms of covariates) and To assess covariate-specific balance, several metrics have been non-additivity (interactions between pretreatment covariates) in proposed in the literature (Austin, 2009; Belitser et al., 2011; Ali the propensity score models. This could be explained by the fact et al., 2014). Each balance metric has its own advantages and that machine learning methods include interactions and square limitations; the absolute standardized difference in means or terms by default (Setoguchi et al., 2008), compared to logistic proportions (ASMD) (Austin, 2009) is more robust in terms of regression where the researcher should “manually” include sample size and covariate distribution requirements in comparison interactions and square terms. When important interaction and to other balance diagnostics, such as overlapping coefficients (Ali square terms are included, the performance of logistic regression et al., 2014; Ali et al., 2015; Ali et al., 2016). The ASMD is also is as good as other machine learning methods (Ali et al., 2017b). a familiar, easy-to-calculate and present, and well-understood statistical tool (Austin, 2009; Ali et al., 2015; Ali et al., 2016). Hence, it is recommended for checking and reporting covariate COVARIATE BALANCE ASSESSMENT balances in propensity score methods (Austin, 2009; Belitser et al., 2011; Ali et al., 2014; Ali et al., 2015; Ali et al., 2016). The ASMD The aim of propensity score methods is to balance covariates is calculated for each covariate and can be averaged to compute an between treatment groups and hence control for measured overall covariate balance and to compare propensity score models confounding (Rosenbaum and Rubin, 1983). Therefore, the (Belitser et al., 2011; Ali et al., 2014). The covariate-specific ASMD quality of propensity score model should be assessed primarily is useful to identify the variable that is still imbalanced and to on the covariate balance achieved. It should not be evaluated modify the propensity score model with squares and interaction based on how well the propensity score model discriminates terms of the variable to improve its balance. Although there is no between treated and untreated individuals, i.e., whether the universal threshold below which the level of covariate imbalance treatment assignment is correctly modeled (Rubin, 2004b; is always acceptable (Imai and Van Dyk, 2004; Ali et al., 2016), Westreich et al., 2011; Ali et al., 2015, Ali et al., 2016) or whether the use of arbitrary cutoffs for balance diagnostics (e.g., < 10% for the subsequent estimates of treatment effect are smaller or larger the ASMD) is common in the medical literature (Ali et al., 2015; than expected (Rosenbaum and Rubin, 1984; Hansen, 2004). Ali et al., 2016). Covariate balance is not only a property of the Hence, propensity score modelling can be considered as an sample means but also of the overall distribution of the covariate; iterative step where the propensity score model is updated by hence, higher-order sample moments of the covariate distribution adding different covariates, interactions between covariates, or such as variance should also be evaluated (Rosenbaum and Rubin, higher-order terms of continuous covariates until an acceptable 1985; Rubin, 2001; Ho et al., 2007; Austin, 2009; Linden and level of balance on important confounding variables is achieved Samuels, 2013). Rubin (2001) proposed the ratio of variances of (Rosenbaum and Rubin, 1984). It is also important to underline treated and untreated groups as an additional check on balance; that variable selection and covariate balance are inseparably a variance ratio of 1.0 in the propensity score matched sample linked; however, covariate balance is often checked on a indicates a good matching and acceptable balance, and a variance preselected list of pretreatment covariates (Ali et al., 2015). On ratio below 2 is generally considered acceptable balance (Rubin, the other hand, there are propensity score modelling techniques 2001; Linden and Samuels, 2013). that optimize covariate balance while estimating the propensity In addition to numerical quantification of the covariate score (Imai and Ratkovic, 2014; Austin, 2019). balance achieved by the specified propensity score model, Frontiers in Pharmacology | www.frontiersin.org 5 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA graphical methods such as (weighted) side-by-side box plots, Once a matched sample has been formed, covariate balance quintile-quintile (Q-Q) plots, plots of ASMD, and empirical can be easily checked between the matched groups using one of density plots of continuous pretreatment covariates provide the balance diagnostics, preferably ASMD, and then treatment a simplified overview on whether balance on individual effect can be estimated by directly comparing outcomes between pretreatment covariates has improved, compared to pre- treated and untreated individuals in the matched sample matching, pre-stratification, or pre-weighting (Rosenbaum (Rosenbaum and Rubin, 1983; Rubin and Thomas, 1996). With and Rubin, 1983; Ali et al., 2016). dichotomous or binary outcomes such as the presence or absence of a disease (“Yes” or “No”), the effect of the treatment can be estimated as the difference or the ratio between the proportion PROPENSITY SCORE METHODS of individuals experiencing the outcome event in each of the two treatment groups (treated vs. untreated) in the matched Once the propensity score has been estimated, researchers sample. If the outcome is continuous, for example blood pressure have several options of using the propensity score in the measurement or HBA1c level, the effect of the treatment is design or analyses, including matching, stratification (also estimated as the difference between the mean outcome for treated called subclassification), covariate adjustment using the and the mean outcome for untreated individuals in the matched propensity score, inverse probability of treatment weighting, sample (Rosenbaum and Rubin, 1983). and combinations of these methods (Rosenbaum and Rubin, If matching is done with replacement or in one-to-many 1983; Rosenbaum and Rubin, 1984; Rubin and Thomas, 2000; matching, weights should be incorporated to account for the Hirano and Imbens, 2001; Johnson et al., 2018). Each method multiple use of the same untreated individual to match with has its own advantages and disadvantages; the choice of a several treated individuals or the multiple use of the same specific propensity score method is in part determined by the treated individual to match with several untreated individuals, inferential goal of the research (i.e., the type of treatment effect respectively (Stuart, 2010). Whether or not to account for the estimand: the average treatment effect in the entire population, matched nature of the data in estimating the variance of the ATE, versus the average treatment effect in the treated treatment effect, for example, using paired t-test for continuous population, ATT) (Imbens, 2000; Stuart, 2008; Ali et al., 2016). outcome or McNemar’s test for binary outcome, is an ongoing Although it is possible to estimate both ATT and ATE using all discussion (Schafer and Kang, 2008; Stuart, 2008; Austin, 2008a; of the four propensity score methods, for example, by assigning Austin, 2011). different weights for the treated and untreated individuals, the The most appealing feature of propensity score matching default approach in each method might give slightly different is that the analysis can partly mimic that of an RCT, meaning estimand. For example, propensity score matching primarily that the distribution of measured pretreatment covariates will estimates the treatment effect in the treated group, ATT be, on average, similar between treatment groups. Hence, direct (Imbens, 2004; Stuart, 2008). Therefore, to get an estimate of comparison of outcomes between treated and untreated groups the average treatment effect in the entire population, ATE, within the propensity score matched sample has the potential one has to use either full matching (Hansen, 2004) or different to give unbiased estimate of the treatment effect, depending on weighting (Stuart, 2008, Stuart, 2010; Ali et al., 2015; Ali et al., the extent to which the measured variables have captured the 2016). The use of a specific propensity score method has potential confounding factors (Rosenbaum and Rubin, 1983). also direct implication on the covariate balance assessment However, RCT, on average, guarantees balance on both measured (Rosenbaum and Rubin, 1983; Rosenbaum and Rubin, 1984; an unmeasured confounders, whereas propensity score improves Ali et al., 2016) and interpretation of the estimated treatment balance on measured confounders but those of unmeasured effect (Stuart, 2008; Ali et al., 2015; Ali et al., 2016). confounders only to the extent that they are related to the measured confounders included in the propensity score model (Rubin, 2004b; Austin, 2011). Other useful features include: Propensity Score Matching separation of the design from analysis via preprocessing of the Propensity score matching, the most common application of data to improve covariate balance without using outcome data, propensity score (Ali et al., 2015), entails forming matched thereby a minimal reliance on model specification; relatively groups of treated and untreated individuals having a similar easy assessment, visualization, and communication of covariate value of the propensity score (Rosenbaum and Rubin, 1983; balance using simple statistics or plots; and qualitative indication Rubin and Thomas, 1996). The matching could be done in many of whether the dataset at hand is good enough to address the ways: one-to-one or one-to-many (1:n, where n is the number causal question without relying on untrustworthy “model- of untreated individuals often up to five), exact or caliper dependent” extrapolations (Rubin, 2004b; Ho et al., 2007; Ali matching, matching with or without replacement, stratified et al., 2016). matching, and full matching (Hansen, 2004). However, one-to- Recently, the use of propensity score for matching has been one caliper matching without replacement is the most common criticized on the basis of an argument that propensity score matching implementation of propensity score matching (Ali et al., 2015; approximates complete randomization and not completely blocked Ali et al., 2016). For detailed discussion on different matching randomization; hence, it engages in random pruning or exclusion approaches, we refer to the literature (Rosenbaum and Rubin, of individuals during matching. “Unlike completely blocked 1985; Hansen, 2004; Stuart, 2010). randomization, random exclusion of individuals in propensity score Frontiers in Pharmacology | www.frontiersin.org 6 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA matching, as in complete randomization, means a decrease in sample individuals will have a similar distribution of measured covariates; size leading to covariate imbalance and more model dependence, hence, the effect of the treatment can be estimated by direct so called the ‘propensity score paradox’” (King and Nielsen, 2016). comparison of outcomes between treated and untreated groups At first this might seem a valid argument; however, the practical within each strata (Rosenbaum and Rubin, 1984; D’Agostino, 2007; implication of this paradox is very limited, if any (Ali et al., 2017a). Ali et al., 2017a). The stratum-specific treatment effects can then This is partly due to the fact that propensity score matching could be aggregated across subclasses to obtain an overall measure of the do better than complete randomization with respect to the balance treatment effect (Rosenbaum and Rubin, 1984). of measured covariates if variables related to treatment are included Rosenbaum and Rubin (1983, 1984) proposed quintile in the propensity score model (Joffe and Rosenbaum, 1999). In stratification on the propensity score based on their finding that addition, the use of matching algorithms such as caliper matching five equal-size propensity score strata removed over 90% of the or matching with replacement retains the best matches thereby bias due to each of the pretreatment covariates used to construct avoiding random pruning or exclusion, and hence the paradox is the propensity score. However, it is recommended that researchers not a big concern. Furthermore, it is currently a standard practice examine the sensitivity of their results to the number of subclasses to check covariate balance in the propensity score matched sample by repeating the analysis using different quantiles of the propensity before estimating the treatment effect, further minimizing any risk score (Imai and Van Dyk, 2004; Adelson et al., 2017). Similar to of exacerbating covariate imbalance (Ali et al., 2015). matching, residual imbalances after stratification can be accounted Similar to RCT, when there are residual differences in pretreatment for using regression adjustment within each stratum (Rosenbaum characteristics between treatment groups in propensity score and Rubin, 1984; Rubin, 2001). Alternatively, the propensity matched sample, regression adjustment can be used on the matched score, defined as quintiles and deciles, can be used as a categorical sample to reduce bias due to residual differences in important variable in a model-based adjustment to estimate treatment effects prognostic factors (Rubin and Thomas, 2000; Imai and Van Dyk, (Rosenbaum and Rubin, 1984; Ali et al., 2016). 2004; Schafer and Kang, 2008). This method has been described as a Propensity score stratification can estimate the stratum-specific doubly robust (DR) approach, i.e., correct specification of either the ATT, or the overall ATT across strata, or the ATE, depending on matching or the regression adjustment, but not necessarily both, is how the subclass treatment effect estimates are weighted. Weighting required to obtain unbiased estimate of the treatment effect (Schafer stratum-specific estimates by the total number of individuals and Kang, 2008; Funk et al., 2011; Nguyen et al., 2017). Propensity (treated and untreated) in each stratum yields the ATE. On the score matching primarily estimates the effect of treatment in other hand, weighting stratum-specific estimates by the proportion the treated individuals (ATT), not the effect of treatment in the of treated individuals in each stratum provides ATT (Stuart, 2010; population (treated and untreated individuals, ATE) (Imbens, Ali et al., 2016). Similarly, pooling stratum-specific variances 2004; Stuart, 2008). This is because the closest untreated and treated provides pooled estimates of the variance for the pooled ATT or individuals are matched and the remaining untreated individuals ATE estimate (Imbens, 2004; Ali et al., 2016). Pooling the stratum- that were not matched are often excluded from the analysis (Stuart, specific treatment effect is straightforward when treatment effect is 2008; Ali et al., 2016). It is important to emphasize that exclusion homogeneous among the propensity score strata (Ali et al., 2016). of unmatched individuals from the analysis not only affects the When there is heterogeneity of treatment effect among the strata precision of the treatment effect estimate but also could have even after automated iterations of the number and boundaries of consequences for the generalizability of the findings, even for the propensity score strata (Imbens, 2004; Imbens and Rubin, 2015; ATT (Lunt, 2013; Ali et al., 2016). For example, exclusion of treated Ali et al., 2016), pooling the stratum-specific treatment effect individuals due to a lack of closer untreated matches could change might complicate interpretation of the treatment effect estimate the estimand from the effect of treatment in the treated (ATT) to the (Ali et al., 2014; Ali et al., 2016). In the presence of treatment effect effect of treatment in those treated individuals for whom we can find modification regardless of the presence of confounding, Mantel- untreated matches (ATT) (Lunt, 2013; Ali et al., 2016). However, it Haenszel methods do not estimate a population parameter (ATE); is possible to estimate the ATE in the matched sample with slight hence, estimating the effect of treatment in the treated (ATT) rather modifications of the matching algorithms. For example, using full than the whole population (ATE), for example, using propensity matching that retains all the treated and untreated individuals in score matching is preferable (Stürmer et al., 2006b). Alternatively, the study data, one can estimate either the ATE or ATT (Hansen, one could standardize the stratum-specific estimates to a specified 2004; Stuart, 2010). Generally, matching discards some data (often distribution of propensity scores, for example, to calculate a unmatched untreated individuals); however, it may increase the standardized mortality ratio (AMR) from the stratum-specific efficiency, reducing the estimated standard error, of the treatment estimates (Stürmer et al., 2006b; Lunt et al., 2009) effect estimate by reducing heterogeneity of observations (Ho et al., Stratification has several advantages: it is an easy and well- 2007; Ali et al., 2016). understood method to implement; it is straightforward to evaluate and communicate covariate balance, and to interpret particularly to non-technical audiences; it separates the design of the study from Propensity Score Stratification the analysis, like propensity score matching, hence less dependent Propensity score stratification, also called propensity score on parametric models (Rosenbaum and Rubin, 1984); it is less subclassification, involves grouping individuals into strata based sensitive to nonlinearities in the relationship between propensity on their propensity scores (often 5 groups using quintiles or 10 scores and outcomes; and it can accommodate additional model- groups using percentiles). Within these strata, treated and untreated based adjustments (Rosenbaum and Rubin, 1983; Rosenbaum and Frontiers in Pharmacology | www.frontiersin.org 7 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA Rubin, 1984). However, this propensity score approach is prone to treated and untreated individuals, respectively (Hernán et al., 2000; residual confounding, which might be an issue due to propensity Ali et al., 2016). score heterogeneity within the strata. Alternative approaches such as weight trimming and weight truncation have been suggested (Cole and Hernán, 2008; Lee et al., 2011). Weight trimming involves removing individuals in the tails Regression Adjustment Using Propensity of the propensity score distributions using percentile cut-points Score (Cole and Hernán, 2008; Lee et al., 2011), i.e., individuals who have The propensity score, as a single summary of all covariates included extreme values of the propensity score—both very high and very low in the propensity score model, can be included as a covariate in a are excluded. On the other hand, weight truncation involves setting regression model of the treatment, i.e., the outcome variable is a maximum allowable weight, Wma), such that individuals with regressed on the treatment variable and the estimated propensity a weight greater than Wma) will be assigned Wma) instead of their score (Rosenbaum and Rubin, 1983; Ali et al., 2016). Although this actual weights. Both approaches may help stabilize weights, reduce approach is very easy to implement, it is generally considered to be a the impact of extreme observations, and can improve the accuracy sub-optimal application of the propensity score for several reasons: and precision of parameter estimates; however, both involve bias- 1) The treatment effect estimation is highly model-dependent variance trade-offs (Lee et al., 2011). For example, trimming the tails because it mixes the study design and data analysis steps; hence, it excludes some individuals with extreme values and hence changes requires correct specification of the propensity score model (Rubin, the population, which might introduce bias depending on the cut- 2004b; Johnson et al., 2018). 2) It also makes additional assumptions off (Cole and Hernán, 2008). Recently, Li et al. (2018) suggested a unique to regression adjustment; the relationship between the different set of weights called “overlapping weights” which weight estimated propensity score and the outcome must be linear and each individual proportional to its probability of receiving the there should be no interaction between treatment status and the alternative treatment, i.e., the overlap weight is defined as 1-PS propensity score (Rosenbaum and Rubin, 1983; Austin, 2011; Ali for a treated individual and PS for an untreated individual. Unlike et al., 2016). However, both assumptions can be checked with the standard IPW, the overlap weights are bounded between 0 and 1; data, and can be relaxed if necessary, for example, by combining with hence, they are less sensitive to extreme weights. It also means that propensity score stratification. 3) It enables estimation of the ATE; there is no need for arbitrary choice of a cut-off for inclusion in the however, its interpretation is complicated particularly in nonlinear analysis as well as exclusion of individuals, unlike weight trimming models such as logistic regression or Cox regression where the (Li et al., 2018). estimand of interest is non-collapsible. Non-collapsibility refers to In the weighted population, weighted standardized difference can a phenomenon in which, in the presence of a non-null treatment be used to compare means, proportions, higher-order moments, and effect, the marginal (overall) treatment effect estimate is different interactions between treated and untreated individuals. In addition, from the conditional (stratum-specific) treatment effect estimate, graphical methods can be employed to compare the distribution of even in the absence of confounding (Greenland et al., 1999; Austin, continuous covariates between treated and untreated individuals 2008b). In addition, assessment and communication of covariate (Austin and Stuart, 2015). Once sufficient covariate balance is balance are not straightforward (Ali et al., 2016). achieved, the effect of the treatment can be estimated by direct comparison of outcomes between treated and untreated groups. The weights can also be used in weighted regression models to Inverse Probability Treatment Weighting estimate the effect of the treatment; and adjustment can be made for Inverse probability weights (IPW) calculated from propensity score covariates that are not sufficiently balanced in the weighted sample. can also be used to create a weighted “artificial” population, also This method focuses on estimating the average treatment effect in called a “pseudo-population” in which treatment and measured the entire population (ATE); modification of the weights allows pretreatment characteristics included in the propensity score are to estimate the average treatment effect in the treated population independent (Hernán et al., 2000; Robins et al., 2000; Cole and (ATT) (Stuart, 2010; Ali et al., 2016). Most importantly, the variance Hernán, 2008; Ali et al., 2016). Hence, treated individuals will be estimation should take into account the weighted nature of the assigned weights equal to the inverse of their propensity scores (1/ “pseudo-population” since some observations can have weights that PS, as they have received the treatment) and untreated individuals are unequal to one another (hence, potentially inducing a within- will be assigned weights equal to the inverse of one minus their individual correlation in outcomes), for example, by using the sample propensity scores [1/(1 – PS)] (D’Agostino, 2007). A particular weights in robust variance estimation (Hernán et al., 2000; Cole and diagnostic concern in using propensity score weighting is that Hernán, 2008; Austin and Stuart, 2015). Alternatively, bootstrapping individuals with extremely large weights may disproportionately could be used to construct 95% confidence intervals, which also influence results and yield estimates with high variance (Lee et al., takes into account the estimation of the propensity score, in addition 2011). When some individuals have probabilities of receiving to the lack of independence between duplicate observations in the the treatment close to 0 or 1, the weights for such individuals weighted sample (Hernán et al., 2000; Austin and Stuart, 2015; Ali et become extremely high or extremely low, respectively (Ali et al., al., 2016; Ali et al., 2017b). 2016). Weight stabilization to “normalize” the range of the inverse Inverse probability of treatment weights (IPTW) can be also probabilities is often considered: the “1” in the numerator of the be used to estimate parameters of marginal structural models inverse probability weights can be replaced with the proportion of (MSMs) to deal with time-varying confounding (Hernán et treated individuals and the proportion of untreated individuals for al., 2000), time-modified confounding (Platt et al., 2009), and Frontiers in Pharmacology | www.frontiersin.org 8 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA competing risks (Hernán et al., 2000; Ali et al., 2017b). Hence, the plan, at each subsequent visit (t = 1,2,…), is time-varying: the implementation of propensity scores as inverse probability weights clinician in consultation with the individual decides treatment is often referred to as MSM using IPTW. All other propensity score AZTt based on the changing values of the individual’s clinical and approaches can only be extended to time-varying confounding demographic history recorded during the previous and current and treatment settings under certain conditions as described in visits. These include prior treatment history, current CD4 count, Figure 2. Comparison of the four propensity score approaches is and other confounders, which are not included in this discussion summarized in Table 1. and ignored for now for the sake of simplicity. The relationships between treatment, confounder, and outcome are presented using directed acyclic graphs (DAGs) for clarity. EXTENDED APPLICATIONS In Figure 2, we considered two time points or visits t = 0 (baseline/pretreatment) and t = 1; hence, CD40 refers to Time-Varying Treatments baseline CD4 count and AZT0 refers to treatment at the first In clinical practice, it is common for patients to start on a visit. Treatment decision at the first visit AZT0 is influenced certain medication, stop or switch to another one (for example, by pretreatment CD4 count (CD40), represented in Figure due to intolerance or lack of adequate response); in such cases, 2A by the arrow from CD40 to AZT0. In the second visit treatment might be treated as a time-varying exposure. Consider (t = 1), treatment decision AZT1 is based on previous a cohort study to estimate the effect of antiretroviral zidovudine treatment (AZT0) and CD4 count at the current visit (CD41), treatment (AZT) in HIV (human immunodeficiency virus) represented in Figure 2A by the arrows from AZT0 and CD41 positive individuals, on progression to AIDS (acquired immune to AZT1. deficiency syndrome), where CD4 count is a confounder. In settings such as DAG of Figure 2A, where there is no arrow Assuming individuals show up for clinical visits at baseline/ from AZT0 to CD41 implying previous treatment does not affect pretreatment (t = 0) and then every 6 months (t = 1, 2, 3,…), current CD4 count, all the standard propensity score approaches and CD4 counts are recorded at these visits (CD4t), represented can deal with the time-varying confounder CD4 count by matching, as CD40, CD41, CD42,…). If AZT is a time-varying dichotomous conditioning, stratification, or weighting, for example, by combining treatment variable indicating whether the individual is on with time-varying Cox models to estimate the treatment effect. antiretroviral treatment at each of the visits (AZTt, represented However, this is not biologically plausible; RCTs have proved that as AZT0, AZT1, AZT2,…), this means, an individual’s treatment antiretroviral treatment indeed affects CD4 count. It is important FIGURE 2 | Causal diagrams representing time-varying treatment (AZT), outcome (progression to AIDS, AIDS), and time-varying confounding (CD4 count). Time-varying confounding is not affected by prior treatment (A), time-varying confounding is affected by prior treatment (B), time-varying confounding affected by unmeasured factor U, which is also associated with the outcome (C), and conditioning or stratifying on time-varying confounder, indicated by box around CD41), creates association between time-varying confounder CD40) and unmeasured factor U (D). Frontiers in Pharmacology | www.frontiersin.org 9 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA TABLE 1 | Comparison of the different propensity score methods. Characteristics Matchinga Stratificationb Regressionc IPTWd Model dependence Minimum Minimum High Minimum Application1 Easy Easy Easy Complex Overall transparency High High Low Medium Easy to communicate Yes Yes Not always Not always Design and analysis Separated Separated Mixed Separated Easy to check balance Yes Yes No Yes Requires unique assumption2 No No Yes No Excluded individuals from analysis3 Yes No No Yes-No Variance estimation Not clear Easy Easy Complex Easy to interpret4 Not always Yes No Often ”Propensity score paradox” Sensitive No No No Estimand5 Often ATT ATE, ATT ATE ATE, ATT Time-varying confounding6 No No No Yes Multiple treatments Possible Complex Complex Easier Multi-level treatment applications Exist Exist None Exist Treatment effect modification Easier Complex Easier Complex aConstructs treated and untreated matched groups with similar propensity scores. bConstructs subgroups of treated and untreated individuals, often quintiles or deciles of PS. cPS is used, as a single summary of all covariates included in PS model, in regression model. dPSs are used as weights to create a pseudo-population in which exposure and measured covariates included in the treatment (PS) model are independent (Ali et al., 2016). 1Estimation of stabilized weights as well as extension to time-varying treatment and confounding setting in MSMs framework can be complex (Ali et al., 2016). 2Requires correct specification of PS and outcome model, apart from the basic assumptions that there is positivity and no unmeasured confounding (Ali et al., 2016). 3Weight trimming excludes some individuals in the tails of the propensity score distribution. 4In PSM, when treated individuals are excluded, interpretation of the treatment effect may change, not just ATT and in Stratification, when there is treatment effect modification by the PS, in regression adjustment using PS, when non-collapsible effect measures such as odds ratios are used. 5Modification of the matching or weighting method enable to estimate either ATT or ATE. 6When time- varying confounder is affected by previous treatment, all the propensity score based methods fail to correctly control for the confounding bias including standard IPWs; however, MSMs using IPWs. to mention that there are many practical examples where both in the DAG of Figure 2D, leading to collider-stratification bias treatment and confounders are time-varying or dynamic, but (Hernán et al., 2000; Cole et al., 2009; Ali et al., 2013). previous treatment does not affect time-varying confounder; hence, In such settings, MSM using inverse probability weighting the DAG in Figure 2A may still be valid in other situations. is the method of choice; unlike conditioning or stratification, When a time-varying confounder (such as CD4 count in our weighting creates a “pseudo-population” in which the association example, CD41) is affected by previous treatment (AZT0) as in the between the time-varying confounder and treatment is removed DAG of Figure 2B, the time-varying confounder (CD41) is also an (Hernán et al., 2000; Robins et al., 2000). Additional methods “intermediate” for the effect of previous treatment (AZT0) on the are also available to deal with time-varying treatment and outcome (progression to AIDS), represented by the path AZT0 → confounding including other classes of marginal structural CD41 → AIDS. Furthermore, if there is an unmeasured common models (g-formula and g-estimation of structural nested models) cause (U) of both the time-varying confounder (CD41) and the (Hernán et al., 2000; Robins et al., 2000). outcome (progression to AIDS) as in DAG of Figure 2C, the time- It is straightforward to hypothesize that such a time-varying varying confounder (CD41) is also a “collider” on the path AZT0 → confounding can also be time-modified, which means not only the CD41 ← U → AIDS (the arrows from U and CD40 collide on CD41). confounder (CD4 count) change over time but also its association Hence, the path AZT0 → CD41 ← U → AIDS is a closed or non- with the treatment and its impact on the outcome (progression causal path because it is blocked at CD41 (using DAG terminologies). to AIDS) varies during these times. The effects of the confounder It also means that there is no association between AZT0 and U unless change over time mean that the strength of association between one conditions, matches, or stratifies on this collider, CD41 (Hernán CD40 and AIDS (CD40 → AIDS) is different from that of CD41 and et al., 2000; Robins et al., 2000). Such a time-dependent variable is a AIDS (CD41 → AIDS) (Platt et al., 2009). However, time-modified confounder, an intermediate, and also a collider all at the same time; confounding might still exist in longitudinal treatment settings hence, adjustment requires careful consideration. where the confounder is time-invariant or fixed. Standard methods Conventional statistical approaches including propensity score are sufficient to deal with time-modified confounding unless the methods (matching, stratification, and regression adjustment) confounder is both time-varying and affected by previous treatment, that condition or stratify on such a covariate will result in a biased which requires the implementation of marginal structural models, estimate of the treatment effect (Hernán et al., 2000; Robins et al., such as using inverse probability weighting. 2000). This happens because conditioning or stratifying on an intermediate will adjust away the indirect effect of the treatment mediated by the cofounder, in this case CD41; and conditioning or Multiple Treatments stratifying on a collider creates a spurious association between the Propensity score methods are often used to estimate the effect of treatment and the unmeasured common cause that did not exist a binary treatment (whether treatment is received: Yes = 1 or No before conditioning (creating an open backdoor path AZT0 → = 0) in observational data. However, with more than two levels of CD41,... U→ AIDS), which is indicated by using dotted lines treatment, which is common in pharmacoepidemiology such as Frontiers in Pharmacology | www.frontiersin.org 10 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA comparison of three or more statins (e.g., simvastatin, atorvastatin, safety are rare. However, simulations studies have shown that fluvastatin, lovastatin, pravastatin, and rosuvastatin) or of multiple multilevel propensity score matching (Arpino and Mealli, 2011) doses of a certain medication (e.g., low, medium and high doses), and weighting approaches (Li et al., 2013), without imposing a estimation of treatment effects requires additional assumptions and within-cluster matching or weighting requirement, reduce bias modelling techniques (Imbens, 2000; McCaffrey et al., 2004). These due to unmeasured cluster-level confounders. include the use of multinomial logistic and multinomial probit models for nominal treatments and ordinal logistic regression or the proportional odds model for ordinal treatments (Imbens, Missing Data 2000). Alternatively, generalized boosted model, a machine learning Missing data is a common problem in the estimation of approach involving an iterative process using multiple regression treatment effects using routinely collected data. The impact trees to capture complex, nonlinear, and non-additive relationships of such missing data on the results of the treatment effect between treatment assignment and pretreatment covariates without estimation depends on the mechanism that caused the data the risk of over-fitting the data, can be used to fit inverse probability to be missing and the way missing data are handled. Missing weighting for multiple treatments (McCaffrey et al., 2004). However, data can be categorized into three distinct classes based on applications in pharmacoepidemiology using observational data the relationship between the missing data mechanism and the are infrequent partly due to methodological complexities in fitting missing and observed values: i) Missing Completely at Random the models and understanding their assumptions as well as limited (MCAR), when the missing data mechanism is unrelated to the availability of guidance documents on these methods. values of any variable, whether missing or observed. Hence, the observed values are representative of the entire sample without missing values. ii) Missing at Random (MAR), when Multilevel Treatments the missing data mechanism is unrelated to the missing values Propensity score methods have been extensively studied and but may be related to the observed values of other variables. widely applied in a single-level treatment (no clustering among iii) Missing Not at Random (MNAR), when the missing data participants); however, most healthcare data have a multilevel mechanism is related not only to the observed values of other structure such that individuals are grouped into clusters such as variables but also to the missing values (Rubin, 1996). For each geographical areas, treatment centers (hospital or physician), or of the missing data patterns, different statistical techniques are insurance plans (Goldstein et al., 2002). The unknown mechanisms used to correct for its impact on the quality of the inference. It that assign individuals to clusters may be associated with individual- is important to emphasize that MCAR, MAR, and MNAR could level measured confounders (such as race, age, and clinical exist for different variables in a specific data. However, if one characteristics) and unmeasured confounders (such as unmeasured variable is MAR or MNAR, generally, the dataset is considered severity of disease, aggressiveness in seeking treatment) (Li et al., MAR or MNAR, respectively. 2013). These measured and unmeasured confounders might also Complete case analysis, including only those individuals create a cluster-level variation in treatment and/or outcome. If this who have no missing data in any of the variables that are variation is correlated with group assignment at the group or cluster required for the analysis, performs well when data are MCAR level, it might lead to confounding (Greenland, 2000; Li et al., 2013). and may be valid under some MAR and MNAR conditions. Hence, the use of standard regression or propensity score methods However, it often results in biased estimate of the treatment ignoring the cluster structure should be avoided. This is because effect if missing is at random (MAR) (Rubin, 1996; Sterne ignoring the cluster structure often leads to invalid inferences: not et al., 2009). In MAR, as stated before, any systematic only the standard errors are inaccurate but also the cluster-level difference between the missing values of a variable and the effects could be confounded with individual-level effects. observed values of the variable can be explained by differences Propensity score matching and weighting are often used in in observed data (Sterne et al., 2009). Furthermore, missing such settings (Arpino and Mealli, 2011; Li et al., 2013). One data in several variables often lead to exclusion of a substantial might consider the use of within-cluster PSM (of treated and proportion of the original sample, which leads to a substantial untreated individuals), which automatically achieves perfect loss of precision (i.e., power) and hence estimates with wider balance on all the measured cluster characteristics. However, it confidence intervals (Cummings, 2013). Other approaches to is very unlikely, particularly in small clusters, to find a sufficient deal with missing data include: 1) replacing missing values number of untreated matches to treated individuals in the same with values imputed from the observed data (for example, cluster. Alternatively, PSM could be performed across clusters using the mean of the observed values); 2) using a missing taking into account the cluster structure in the propensity score category indicator; and 3) using the last observed value to estimation model. Preferably, cluster structure should be taken replace missing values particularly in longitudinal studies into account in estimation of both the propensity score and the [often called “last observation carried forward” (LOCF)]. treatment effect (Li et al., 2013). These three approaches are generally statistically invalid, Multilevel regression models that include fixed effects and/ except under certain conditions, and they might lead to or random effects have been developed (Greenland, 2000; serious bias (Rubin, 1996; Sterne et al., 2009). Missing category Goldstein et al., 2002), and extended to propensity scores indicator and LOCF approaches require specific assumptions approaches (Arpino and Mealli, 2011). Empirical applications for validity that are distinct from the MCAR, MAR, and of such methods in medication and device effectiveness and MNAR categorization. On the other hand, single imputation of Frontiers in Pharmacology | www.frontiersin.org 11 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA missing values (mean imputation) usually results in too small ADVANTAGES AND LIMITATIONS OF standard errors, because it fails to account for the uncertainty PROPENSITY SCORE METHODS about the missing values (Sterne et al., 2009). A relatively flexible approach to allow for the uncertainty in Previous literature reviews of observational studies have found the missing data is multiple imputation. Multiple imputation that results from both traditional regression and propensity scores involves creating multiple different copies of the dataset analyses are similar (Shah et al., 2005; Stürmer et al., 2006a). These with the missing values replaced by imputed values (Step 1); findings may be in part due to sub-optimal implementations of estimating treatment effects in each copy of the data (Step propensity score methods (Shah et al., 2005; Austin, 2008a; Ali et al., 2); averaging the estimated treatment effects to give overall 2015); however, similarity of findings has been used to question the estimated measure of association and calculating standard need for propensity score methods if they do not provide better ways errors using Rubin’s rules (Step 3) (Rubin, 1996; Rubin, to improve confounding control. Despite these findings, propensity 2004a). Applications of propensity score methods in data score methods will remain advantageous for several reasons with missing values involve a similar approach: 1) creation of compared to covariate-adjustment techniques, which correct for multiple copies of imputed data; 2) estimation of propensity covariate imbalances between treatment groups by conditioning scores and treatment effects in each of the imputed copies of them in the regression model for the outcome. the dataset (Qu and Lipkovich, 2009; Leyrat et al., 2019); and 3) pooling of treatment effects by averaging across the multiple datasets and estimation of standard errors using Ruben’s Transparency rule (Crowe et al., 2010; Leyrat et al., 2019) (Figure 3A). An Propensity score methods primarily aim at balancing treatment alternative approach is pooling the propensity scores from groups with respect to covariate distributions; when sufficient the multiple copies of data, in step 2, and conducting the covariate balance is achieved, it is relatively easy to check and analysis in the pooled data (Figure 3B); however, this method communicate the balance (Ali et al., 2015, Ali et al., 2016) by has been proved sub-optimal in terms of bias reduction using simple graphical tools or quantitative statistics. In addition, (Leyrat et al., 2019). propensity score methods, unlike regression adjustment, can FIGURE 3 | Multiple imputation in propensity score methods; multiple copies of imputed data are created and propensity score is estimated using these datasets. Treatment effects are estimated in several datasets (A) and propensity scores from multiple datasets are pooled and treatment effect estimated in a single dataset (B). *Other PS methods, stratification, IPTW, and covariate adjustment using PS could also be used instead of matching. Frontiers in Pharmacology | www.frontiersin.org 12 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA give investigators an insight into the quality of the data at hand. score when there are fewer than eight outcomes per included Inadequate overlap in propensity score distributions (also called covariate to effectively improve estimation. poor “common support”) between treatment groups should be considered as a warning that the data set at hand may not be sufficient to reliably address the causal question without “model- Doubly Robust Estimations Generally, doubly roubst estimations (DR) estimation methods dependent” extrapolations based on untrustworthy assumptions apply different procedures or models simultaneously and (Dehejia and Wahba, 2002; Rubin, 2004b; Rubin, 2007; Ali et al., produce a consistent estimate of the parameter if either of the two 2016). In some cases, the researcher might decide to focus on models, not necessarily both, has been correctly specified (Imai individuals only in the overlapping regions using propensity score and Ratkovic, 2014). Several applications of propensity scores matching or trimming; as a consequence, the conclusions of the have been described as DR in terms of estimating the effect of a findings should be restricted to individuals that are sufficiently certain treatment, including: represented in the overlapping regions of the propensity score distributions (Ali et al., 2016). Conventional regression 1) The combined use of propensity score methods (matching, methods do not provide the researcher with these possibilities. regression, or weighting) with regression adjustments. These Furthermore, covariate balance in regression methods is a approaches use non-parametric pre-processing of the data “black-box” and, irrespective of inadequate overlap (i.e., when to minimize imbalances in measured covariates and, if there the treated and untreated groups are disparate on pretreatment are still residual differences, the covariates can be adjusted covariates), conventional models use extrapolations to estimate in the outcome model (Rubin and Thomas, 2000; Nguyen treatment effects that may not be generalizable to the entire et al., 2017). population in the data set. 2) The combined use of propensity and prognostic score methods (Leacy and Stuart, 2014; Ali et al., 2018b); a Design Tools prognostic score is any function of a set of covariates that when Similar to RCTs, propensity score methods can be considered as conditioned on creates independence between the potential design tools for pre-processing of the data (matching, stratification, outcome under the control (no treatment) condition and the and weighting) without using any outcome information at this stage. unreduced covariates (Hansen, 2008). Hence, differences in As a result, formal causal inference models (also called the potential outcomes between treated and untreated individuals can outcomes framework) (Rubin, 2005) can be applied to clearly be attributed to the effect of the treatment under study. The specify the causal question without conflating with the modeling two approaches could be combined in several ways such approach (Vandenbroucke et al., 2016); hence, it allows for a simple as full matching on a Mahalanobis distance combining the and transparent analysis. In addition, this approach minimizes estimated propensity and prognostic scores; full matching bias from potential misspecification of the outcome model (Rubin, on the estimated prognostic score within propensity score 2004b). Furthermore, matched, stratified, and weighted analyses calipers; and subclassification on an estimated propensity do not make strong assumptions of linearity in the relationship and prognostic score grid with five subclasses, among of propensity score with the outcome. If a non-parametric pre- others (Leacy and Stuart, 2014; Ali et al., 2018b). Methods processing of the data using propensity score methods does not combining propensity and prognostic scores were no less reduce model dependence, it is reasonable to accept that the data robust to model misspecification than single-score methods do not have enough information to reliably support the causal even when both prognostic and propensity score models inference by any other statistical method. In fact, this knowledge in were incorrectly specified in simulation and empirical itself should still be useful and the conclusion may be correct (Rubin, studies (Leacy and Stuart, 2014). 2004b; Ho et al., 2007; Rubin, 2007; Ali et al., 2016). 3) The use of covariate balancing propensity score (CBPS) introduced by Imai and Ratkovic (2014) involves estimation Dimension Reduction of the propensity score such that the resulting covariate Propensity score typically summarizes a large number of balance is optimized. This approach utilizes the dual measured pretreatment covariates to a single score; hence, characteristics of the propensity score as a covariate it is called a “summary score.” This is particularly useful in balancing score and the conditional probability of treatment high-dimensional data with a substantially large number of assignment. Specifically, “the covariate balancing property pretreatment covariates compared to the number of outcome (i.e., mean independence between the treatment status and events including rare events, typical of most medication measured covariates after inverse propensity score weighting) safety studies in pharmacoepidemiology (Glynn et al., 2006). is used as condition to imply estimation of the propensity In this setting, maximum likelihood estimations used in score while also incorporating the standard estimation conventional regression techniques such as logistic and Cox procedure” (Imai and Ratkovic, 2014). Unlike other covariate regression require several outcome events for each parameter balancing methods, a single model determines the treatment included in the regression model; the rule of thumb is that assignment mechanism and the covariate balancing weights. ≥ 10 outcome events are required per confounder included Once CBPS is estimated, various propensity score methods in a model (Peduzzi et al., 1995; Peduzzi et al., 1996). On the such as matching and weighting can be implemented without other hand, Cepeda et al. (2003) suggested using propensity modification (Imai and Ratkovic, 2014). The basic idea of Frontiers in Pharmacology | www.frontiersin.org 13 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA CBPS is optimizing covariate balance so that even when the variables that are available to the researcher. The only way propensity score model is misspecified, there will still be a to convince a critical reader that the study is not subject to reasonable balance of the covariates between the treatment unmeasured confounding is to have a rich set of covariates and comparator groups. Unlike standard DR estimators, for constructing the propensity score model. Therefore, it however, the CBPS approach does not require estimation of is important to provide a detailed account of the variables the outcome model. collected and included in the propensity score model (Ali 4) Calculation of DR estimators using different approaches, et al., 2015). for example, using the propensity score, predicted, and Modifications of the standard propensity score applications observed outcome (Ŷ and Y, respectively). This approach have been suggested to further reduce the risk of unmeasured involves specifying regression models for the treatment confounding including the use of high-dimensional propensity (Z) and the outcome (Y) as a function of covariates (X) score and propensity score calibration. High-dimensional and combining these subject-specific values to calculate propensity score refers to the use of a large number (in the the DR estimate for each individual. First, treatment is range of several hundreds) of covariates to improve control of modelled as a function of covariates to estimate propensity confounding; the underlying assumption is that the variables scores for each individual using the observed data. Second, may collectively be proxies for unobserved confounding factors the relationships between measured confounders and the (Schneeweiss et al., 2009; Rassen et al., 2011). Propensity score outcome are modelled within treated and untreated groups calibration refers to the use of a “gold standard” propensity score separately. The resulting parameter estimates are then used estimated in a separate validation study, with more detailed to calculate predicted outcomes (Ŷ1, Ŷ0) for each individual covariate information unmeasured in the main study, to correct in the population that is treated (setting Z = 1) and not the main-study effect of the drug on the outcome (Stürmer et al., treated (setting Z = 0) given covariate values. Third, 2005; Stürmer et al., 2007). the DR estimates of the outcome are calculated for each Furthermore, sensitivity analyses (Rosenbaum and Rubin, individual both in the presence and absence of treatment 1983; Rosenbaum, 2005) are useful to assess the plausibility of (DR1 and DR0), respectively) using the subject-specific the assumptions underlying the propensity score methods and predicted (Ŷ) and observed (Y) outcomes weighted by the how violations of them might affect the conclusions drawn propensity score. Finally, the means of DR1 and DR0 are (Stuart, 2010). Methods to deal with unmeasured confounding calculated across the entire study population and these are summarized in Figure 1. means will be used to calculate the effect of the treatment (Funk et al., 2011). Effect Modification In estimating treatment effects, there is often an interest to explore if the effect of treatment varies among different Unmeasured Confounding subgroups (for example, men versus women) of the population Propensity score methods, like other conventional regression under study, often called “treatment effect modification.” methods, can account for only measured confounding factors There are many ways to utilize propensity score methods and not unmeasured factors (Rosenbaum and Rubin, 1983). to adjust for confounding in a subgroup analysis; however, Therefore, propensity score analyses are only as good as common implementation of propensity score matching in the the completeness and quality of the potential confounding medical literature is sub-optimal (Wang et al., 2017; Ali et al., FIGURE 4 | Methods to assess treatment effect modification in propensity score matching. Frontiers in Pharmacology | www.frontiersin.org 14 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA TABLE 2 | Summary of considerations when planning, conducting, and reporting propensity score analysis. Characteristics What to consider Methods available to deal with What should or should not be done Missing data Missing data mechanism Multiple imputation if missing at Avoid complete case analysis and missing random (MAR) indicator category, the later may be biased even when MCAR assumption holds. Variable selection Potential confounders, intermediates, Clinical knowledge/expert opinion. Avoid adjusting for intermediates, colliders, and colliders strong instrumental variables the later (only when sure or suspect strong unmeasured confounding). Association between variables with Avoid the use of p-values, or step-wise variable outcome (and treatment). selection methods. Balance diagnostics. Propensity score estimation Variables included, interactions and Logistic regression, Recursive Report on the method used for estimation higher order terms. partitioning, Neural network, and variables included in the propensity score Classification and regression trees, method. Random forest, and Boosting regression. Propensity score methods The research question, the treatment Density plots of propensity scores. Report the density plots or histograms in effect estimand, and the extent of the propensity score distribution (preferably overlap. overlapping coefficients of the density plots). Propensity score matching Matching algorithm, matching Exact (coarsened) matching, nearest Report on the number of starting population, with or with our replacement, and neighbor matching (with or without number matched, and number excluded (with matching ratio caliper), stratified matching, and their pre-treatment characteristics). full matching. Matching ratio can be: 1-to-1 matching, 1-to- many matching, variable ratio matching, and full matching. Propensity score stratification Number of strata Deciles and quintiles of propensity Report on the number of strata used and the scores. covariate balance between treatment groups in each strata. Regression adjustment using Linear relationship between the Report on whether linear relationship between propensity score outcome and the propensity core. the outcome and propensity core is checked and is fulfilled. Inverse probability of treatment Whether there is sufficient overlap Weighted regression. Robust variance Report on how weights are calculated, if weighting (positivity). estimation or Bootstrapping for weights are stabilized, the mean weights in both constructing confidence intervals. treatment groups, if trimming has been done. Time-varying exposure Whether there is time-varying Marginal Structural models using If previous treatment affect time-varying confounding, and if any, whether it is IPTW, G-formula and G-estimation of confounding avoid matching, stratification and affected by previous treatment. structural nested models. regression adjustment; apply MSM using IPTW. Treatment effect modification Identify potential effect modifier. Matching on PS within strata of effect Avoid the use of stratified analysis using the PSM modifier, among others. data without adjustment for covariates. Multilevel treatment Whether multilevel structure exists Multilevel propensity score methods. Avoid use of single-level propensity score in the data, the number of clusters/ applications. Include multilevel structure at levels least in propensity score estimation or outcome analysis, preferably in both. Multiple treatments Number of treatment groups, whether Multiple matching and weighting: there is order in the treatment multinomial logistic regression, ordinal categories (such as dosage). logistic regression, or generalized boosted model. Residual Confounding Whether there is imbalance in Doubly robust methods, propensity Report on which method was used and why? covariates. score calibration (PSC), high dimensional propensity score (HDPS) method. Unmeasured confounding Whether there is potential Alternative methods such as Report on the method used and the sensitivity unmeasured confounding, or instrumental variable methods, PSC, analysis conducted. whether the data contain proxies for HDPS, or consider sensitivity analysis. unmeasured confounding. 2018a). The use of propensity score matched (PSM) cohort for covariates in the outcome model or ii) re-matching within the subgroup analysis breaks the matched sets and might result subgroups either using the propensity score estimated in the full in imbalance of covariates (Ali et al., 2018a). Depending on cohort or fitting new propensity score within subgroups (Figure 4) the frequency of treatment or outcome, small changes in the (Rassen et al., 2012; Wang et al., 2017). The choice of a specific matched cohort might lead to large fluctuations for measures method should take into account several factors: prevalence of of association (Rassen et al., 2012). the treatment and the outcome; strength of association between To account for covariate imbalances, subgroup analyses pretreatment covariates and the treatment; the true effect size of propensity score matched cohorts involve: i) adjusting for Frontiers in Pharmacology | www.frontiersin.org 15 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA within subgroups; and the amount of confounding within the CONCLUSION subgroups (Wang et al., 2018). Propensity score methods will remain important design and analytic tools to estimate effects of treatment from observational REPORTING data. Preferably, they should be utilized in the design stage as tools for preprocessing of the data and they should be considered The credibility of any research depends on a critical assessment complementary tools, and not replacements, to conventional by others of the strengths and weaknesses in study design, regression adjustments. In fact, when appropriate, propensity score conduct, and analysis. Hence, transparent and adequate methods should be used in combination with other model-based reporting of critical aspects of propensity score-based analysis regression techniques. In addition, propensity score methods (Ali et al., 2015), like other observational studies, helps readers should not be regarded as magical remedies for the inadequacies of follow “what was planned, what was done, what was found, observational studies such as residual or unmeasured confounding and what conclusions were drawn” (Von Elm et al., 2007). It (Rubin and Thomas, 2000; Ali et al., 2016). The ability of propensity also makes it easier for other researchers to replicate the study score methods to overcome confounding is entirely dependent findings using other data sources and to judge whether and on the extent to which measured variables capture potential how results can be included in systematic reviews (Von Elm confounding. Taking full advantage of these methods requires et al., 2007). Despite substantial methodological developments explicit definition of the research question and appropriate and common applications of the propensity score methods, choice of the propensity score method, transparent and detailed in general, reporting on important features of the propensity description of all subsequent statistical analyses to be conducted, score analysis is poor, incomplete, and inconsistent in the and adequate reporting of the important aspects of the propensity medical literature (Austin, 2008a; Ali et al., 2015; Ali et al., score analyses (Ali et al., 2016). 2016; Wang et al., 2017). This could in part be due to a lack of standards for the conduct and reporting of propensity score based studies in guidelines. Therefore, critical items AUTHOR CONTRIBUTIONS relevant to propensity score analyses should be incorporated MA, DP-A, RF, MB, and LS contributed to the conception and in guidelines on the conduct and reporting of observational design of the study. MA wrote the first draft of the manuscript. studies, such as the STROBE statement (Von Elm et al., 2007; DR and NB wrote sections of the manuscript. All authors Ali et al., 2015) and the ENCePP guide on methodological contributed to manuscript revision and read and approved the standards in pharmacoepidemiology (Blake et al., 2012; Ali submitted version. et al., 2015) to improve the quality of the conduct and reporting of propensity score based studies (Ali et al., 2015; Ali et al., 2016). Table 2 summarizes important consideration when FUNDING planning, conducting, and reporting propensity score analysis and list of items that should be reported are summarized by This work is part of The 100 Million Brazilian Cohort project Ali et al. (2016). funded by the Wellcome Trust. Grant code: 202912/B/16/Z. REFERENCES Ali, M. S., Groenwold, R. H., and Klungel, O. H. (2016). Best (but oft-forgotten) practices: propensity score methods in clinical nutrition research–3. Am. J. Adelson, J. L., McCoach, D., Rogers, H., Adelson, J. A., and Sauer, T. M. (2017). Clin. Nutr. 104, 247–258. doi: 10.3945/ajcn.115.125914 Developing and applying the propensity score to make causal inferences: Ali, M. S., Groenwold, R. H., Pestman, W. R., Belitser, S. V., Hoes, A. W., De Boer, variable selection and stratification. Front. Psychol. 8, 1413. doi: 10.3389/ A., et al. (2013). Time-dependent propensity score and collider-stratification fpsyg.2017.01413 bias: an example of beta 2-agonist use and the risk of coronary heart disease. Ali, M. S., Collins, G., and Prieto-Alhambra, D. (2017a). The “propensity score Eur. J Epidemiol. 28, 291–299. doi: 10.1007/s10654-013-9766-2 paradox”: a threat to pharmaco-epidemiological studies? Pharmacoepidemiol. Ali, M. S., Groenwold, R. H., Pestman, W. R., Belitser, S. V., Roes, K. C., Hoes, A. W., Drug Saf. 26:(Suppl.2):3–636. doi:10.1002/pds.4275 et al. (2014). Propensity score balance measures in pharmacoepidemiology: Ali, M. S., Douglas, I. J., Williamson, E., Prieto-Alhambra, D., and Smeeth, L., a simulation study. Pharmacoepidemiol. Drug Saf. 23, 802–811. doi: 10.1002/ (2018a). “Evaluation of treatment effect modification in propensity score pds.3574 matching: An empirical example,” in Pharmacoepidemiology and drug safety Ali, M. S., Khalid, S., Collins, G., and Prieto-Alhambra, D. (2017b). The comparative vol. 27. (NJ USA: Wiley 111 River St, Hoboken 07030-5774), 25–25. performance of logistic regression and random forest in propensity score Ali, M. S., Douglas, I. J., Williamson, E., Prieto-Alhambra, D., and Smeeth, L. methods: A simulation study. Pharmacoepidemiol. Drug Saf. 26:(Suppl.2):3– (2018b). “A joint application of disease risk score and propensity score to 636. doi:10.1002/pds.4275 control for confounding: A clinical example,” in Pharmacoepidemiology and Ali, M. S., Khalid, S., Groenwold, R., Collins, G. S., Klungel, O., and Prieto- drug safety vol. 27. (NJ USA: Wiley 111 River St, Hoboken 07030-5774), Alhambra, D. (2017c). Instrumental variables to test for unmeasured 27–27. confounding: a precautionary note. Pharmacoepidemiol. Drug Saf. Ali, M. S., Groenwold, R. H., Belitser, S. V., Pestman, W. R., Hoes, A. W., Roes, 26:(Suppl.2):3–636. doi:10.1002/pds.4275 K. C., et al. (2015). Reporting of covariate selection and balance assessment Angrist, J. D., Imbens, G. W., and Rubin, D. B. (1996). Identification of causal in propensity score analysis is suboptimal: a systematic review. J. Clin. effects using instrumental variables. J. Am. Stat. Assoc. 91, 444–455. doi: Epidemiol. 68, 122–131. doi: 10.1016/j.jclinepi.2014.08.011 10.1080/01621459.1996.10476902 Frontiers in Pharmacology | www.frontiersin.org 16 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA Arpino, B., and Mealli, F. (2011). The specification of the propensity score in Glynn, R. J., Schneeweiss, S., and Stürmer, T. (2006). Indications for propensity multilevel observational studies. Comput. Stat. Data Anal. 55, 1770–1780. doi: scores and review of their use in pharmacoepidemiology. Basic Clin. Pharmacol. 10.1016/j.csda.2010.11.008 Toxicol. 98, 253–259. doi: 10.1111/j.1742-7843.2006.pto_293.x Austin, P. C. (2008a). A critical appraisal of propensity-score matching in Goldstein, H., Browne, W., and Rasbash, J. (2002). Multilevel modelling of medical the medical literature between 1996 and 2003. Stat. Med. 27, 2037–2049. data. Stat. Med. 21, 3291–3315. doi: 10.1002/sim.1264 doi: 10.1002/sim.3150 Greenland, S. (2000). Principles of multilevel modelling. Int. J. Epidemiol. 29, 158–167. Austin, P. C. (2008b). The performance of different propensity-score methods doi: 10.1093/ije/29.1.158 for estimating relative risks. J. Clin. Epidemiol. 61, 537–545. doi: 10.1016/j. Greenland, S., and Morgenstern, H. (2001). Confounding in health research. Annu. jclinepi.2007.07.011 Rev. Public Health 22, 189–212. doi: 10.1146/annurev.publhealth.22.1.189 Austin, P. C. (2009). Balance diagnostics for comparing the distribution of baseline Greenland, S., Robins, J. M., and Pearl, J. (1999). Confounding and collapsibility in covariates between treatment groups in propensity-score matched samples. causal inference. Stat. Sci. 14, 29–46. doi: 10.1214/ss/1009211805 Stat. Med. 28, 3083–3107. doi: 10.1002/sim.3697 Hansen, B. B. (2004). Full matching in an observational study of coaching for the Austin, P. C. (2011). An introduction to propensity score methods for reducing sat. J. Am. Stat. Assoc. 99, 609–618. doi: 10.1198/016214504000000647 the effects of confounding in observational studies. Multivariate Behav. Res. Hansen, B. B. (2008). The prognostic analogue of the propensity score. Biometrika 46, 399–424. doi: 10.1080/00273171.2011.568786 95, 481–488. doi: 10.1093/biomet/asn004 Austin, P. C. (2019). Assessing covariate balance when using the generalized Hernán, M. Á., Brumback, B., and Robins, J. M. (2000). Marginal structural models propensity score with quantitative or continuous exposures. Stat. Methods Med. to estimate the causal effect of zidovudine on the survival of hiv-positive men. Res 28, 1365–1377. doi: 10.1177/0962280218756159 Epidemiology 11(5), 561–570. doi: 10.1097/00001648-200009000-00012 Austin, P. C., and Stuart, E. A. (2015). Moving towards best practice when using Hirano, K., and Imbens, G. W. (2001). Estimation of causal effects using propensity inverse probability of treatment weighting (iptw) using the propensity score score weighting: An application to data on right heart catheterization. Health to estimate causal treatment effects in observational studies. Stat. Med. 34, Serv. Outcomes Res. Methodol. 2, 259–278. doi: 10.1023/A:1020371312283 3661–3679. doi: 10.1002/sim.6607 Ho, D., Imai, K., King, G., and Stuart, E. (2007). Matching as nonparametric Bärnighausen, T., Tugwell, P., Røttingen, J.-A., Shemilt, I., Rockers, P., Geldsetzer, P., preprocessing for reducing model dependence in parametric causal inference. et al. (2017). Quasi-experimental study designs series—paper 4: uses and value. Polit. Anal. 15 (3), 199–236. doi:10.1093/pan/mpl013 J. Clin. Epidemiol. 89, 21–29. doi: 10.1016/j.jclinepi.2017.03.012 Imai, K., and Ratkovic, M. (2014). Covariate balancing propensity score. J. R. Stat. Belitser, S. V., Martens, E. P., Pestman, W. R., Groenwold, R. H., De Boer, A., and Soc. Series. B. Stat. Methodol. 76, 243–263. doi: 10.1111/rssb.12027 Klungel, O. H. (2011). Measuring balance and model selection in propensity score Imai, K., and Van Dyk, D. A. (2004). Causal inference with general treatment methods. Pharmacoepidemiol. Drug Saf. 20, 1115–1129. doi: 10.1002/pds.2188 regimes: Generalizing the propensity score. J. Am. Stat. Assoc. 99, 854–866. Blake, K. V., deVries, C. S., Arlett, P., Kurz, X., and Fitt, H. of Centres for doi: 10.1198/016214504000001187 Pharmacoepidemiology Pharmacovigilance, E. N. (2012). Increasing scientific Imbens, G. W. (2000). The role of the propensity score in estimating dose-response standards, independence and transparency in post-authorisation studies: functions. Biometrika 87, 706–710. doi: 10.1093/biomet/87.3.706 the role of the european network of centres for pharmacoepidemiology and Imbens, G. W. (2004). Nonparametric estimation of average treatment effects pharmacovigilance. Pharmacoepidemiol. Drug Saf. 21, 690–696. doi: 10.1002/ under exogeneity: A review. Rev. Econ. Stat. 86, 4–29. doi: 10.1162/​ pds.3281 003465304323023651 Brookhart, M. A., Schneeweiss, S., Rothman, K. J., Glynn, R. J., Avorn, J., and Imbens, G. W., and Rubin, D. B. (2015). Causal inference in statistics, social, and Stürmer, T. (2006). Variable selection for propensity score models. Am. J. biomedical sciences an introduction. Cambridge: Cambridge University Press. doi: Epidemiol. 163, 1149–1156. doi: 10.1093/aje/kwj149 10.1017/​CBO9781139025751 Cepeda, M. S., Boston, R., Farrar, J. T., and Strom, B. L. (2003). Comparison of Joffe, M. M., and Rosenbaum, P. R. (1999). Invited commentary: propensity scores. logistic regression versus propensity score when the number of events is Am. J. Epidemiol. 150, 327–333. doi: 10.1093/oxfordjournals.aje.a010011 low and there are multiple confounders. Am. J. Epidemiol. 158, 280–287. Johnson, S. R., Tomlinson, G. A., Hawker, G. A., Granton, J. T., and Feldman, B. M. doi: 10.1093/aje/kwg115 (2018). Propensity score methods for bias reduction in observational studies of Cole, S. R., and Hernán, M. A. (2008). Constructing inverse probability weights treatment effect. Rheum. Dis. Clin. 44, 203–213. doi: 10.1016/j.rdc.2018.01.002 for marginal structural models. Am. J. Epidemiol. 168, 656–664. doi: 10.1093/ King, G., and Nielsen, R. (2016). Why propensity scores should not be used for aje/kwn164 matching. Copy at. http://j.mp/1sexgVw Download Citation BibTex Tagged Cole, S. R., Platt, R. W., Schisterman, E. F., Chu, H., Westreich, D., Richardson, D., XML Download Paper 378. doi: 10.1017/pan.2019.11 et al. (2009). Illustrating bias due to conditioning on a collider. Int. J. Epidemiol. Leacy, F. P., and Stuart, E. A. (2014). On the joint use of propensity and prognostic 39, 417–420. doi: 10.1093/ije/dyp334 scores in estimation of the average treatment effect on the treated: a simulation Concato, J., Shah, N., and Horwitz, R. I. (2000). Randomized, controlled trials, study. Stat. Med. 33, 3488–3508. doi: 10.1002/sim.6030 observational studies, and the hierarchy of research designs. N. Engl. J. Med. Lee, B. K., Lessler, J., and Stuart, E. A. (2010). Improving propensity score weighting 342, 1887–1892. doi: 10.1056/NEJM200006223422507 using machine learning. Stat. Med. 29, 337–346. doi: 10.1002/sim.3782 Crowe, B. J., Lipkovich, I. A., and Wang, O. (2010). Comparison of several Lee, B. K., Lessler, J., and Stuart, E. A. (2011). Weight trimming and propensity imputation methods for missing baseline data in propensity scores analysis of score weighting. PloS One 6, e18174. doi: 10.1371/journal.pone.0018174 binary outcome. Pharm. Stat. 9, 269–279. doi: 10.1002/pst.389 Leyrat, C., Seaman, S. R., White, I. R., Douglas, I., Smeeth, L., Kim, J., et al. Cummings, P. (2013). Missing data and multiple imputation. JAMA Pediatr. (2019). Propensity score analysis with partially observed covariates: How 167, 656–661. doi: 10.1001/jamapediatrics.2013.1329 should multiple imputation be used? Stat. Methods Med. Res 28, 3–19. Dehejia, R. H., and Wahba, S. (2002). Propensity score-matching methods for doi: 10.1177/0962280217713032 nonexperimental causal studies. Rev. Econ. Stat. 84, 151–161. doi: 10.1162/​ Li, F., Morgan, K. L., and Zaslavsky, A. M. (2018). Balancing covariates via 003465302317331982 propensity score weighting. J. Am. Stat. Assoc. 113, 390–400. doi: 10.1080/​ D’Agostino, Jr, R. B. (2007). Propensity scores in cardiovascular research. 01621459.2016.1260466 Circulation 115, 2340–2343. doi: 10.1161/CIRCULATIONAHA.105.594952 Li, F., Zaslavsky, A. M., and Landrum, M. B. (2013). Propensity score weighting Eichler, H.-G., Abadie, E., Breckenridge, A., Flamion, B., Gustafsson, L. L., with multilevel data. Stat. Med. 32, 3373–3387. doi: 10.1002/sim.5786 Leufkens, H., et al. (2011). Bridging the efficacy–effectiveness gap: a regulator’s Linden, A., and Samuels, S. J. (2013). Using balance statistics to determine the perspective on addressing variability of drug response. Nat. Rev. Drug. Discov. optimal number of controls in matching studies. J. Eval. Clin. Pract. 19, 968–975. 10, 495. doi: 10.1038/nrd3501 doi: 10.1111/jep.12072 Funk, M. J., Westreich, D., Wiesen, C., Stürmer, T., Brookhart, M. A., and Davidian, M. Lunt, M. (2013). Selecting an appropriate caliper can be essential for achieving (2011). Doubly robust estimation of causal effects. Am. J. Epidemiol. 173, 761– good balance with propensity score matching. Am. J. Epidemiol. 179, 226–235. 767. doi: 10.1093/aje/kwq439 doi: 10.1093/aje/kwt212 Frontiers in Pharmacology | www.frontiersin.org 17 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA Lunt, M., Solomon, D., Rothman, K., Glynn, R., Hyrich, K., Symmons, D. P., et al. Rubin, D. B. (1997). Estimating causal effects from large data sets using (2009). Different methods of balancing covariates leading to different effect propensity scores. Ann. Intern. Med. 127, 757–763. doi: 10.7326/​ estimates in the presence of effect modification. Am. J. Epidemiol. 169, 909–917. 0003-4819-127-8_Part_2-199710151-00064 doi: 10.1093/aje/kwn391 Rubin, D. B. (2001). Using propensity scores to help design observational studies: Makady, A., de Boer, A., Hillege, H., Klungel, O., Goettsch, W., and on behalf application to the tobacco litigation. Health Serv. Outcomes Res. Methodol. of GetReal Work Package 1. (2017). What is real-world data? a review of 2, 169–188. rubin2001using. doi: 10.1023/A:1020363010465 definitions based on literature and stakeholder interviews. Value Health 20, Rubin, D. B. (2004a). Multiple imputation for nonresponse in surveys Vol. 81. New 858–865. doi: 10.1016/j.jval.2017.03.008 York: John Wiley & Sons. McCaffrey, D. F., Ridgeway, G., and Morral, A. R. (2004). Propensity score Rubin, D. B. (2004b). On principles for modeling propensity scores in medical estimation with boosted regression for evaluating causal effects in observational research. Pharmacoepidemiol. Drug Saf. 13, 855–857. doi: 10.1002/pds.968 studies. Psychol. Methods 9, 403. doi: 10.1037/1082-989X.9.4.403 Rubin, D. B. (2005). Causal inference using potential outcomes: Design, Myers, J. A., Rassen, J. A., Gagne, J. J., Huybrechts, K. F., Schneeweiss, S., Rothman, modeling, decisions. J. Am. Stat. Assoc. 100, 322–331. doi: 10.1198/​ K. J., et al. (2011a). Effects of adjusting for instrumental variables on bias and 016214504000001880 precision of effect estimates. Am. J. Epidemiol. 174, 1213–1222. doi: 10.1093/ Rubin, D. B. (2007). The design versus the analysis of observational studies for aje/kwr364 causal effects: parallels with the design of randomized trials. Stat. Med. Myers, J. A., Rassen, J. A., Gagne, J. J., Huybrechts, K. F., Schneeweiss, S., Rothman, K. J., 26, 20–36. doi: 10.1002/sim.2739 et al. (2011b). Myers et al. respond to “understanding bias amplification”. Am. J. Rubin, D. B., and Thomas, N. (1996). Matching using estimated propensity scores: Epidemiol. 174, 1228–1229. doi: 10.1093/aje/kwr353 relating theory to practice. Biometrics 52, 249–264. doi: 10.2307/2533160 Nguyen, T.-L., Collins, G. S., Spence, J., Daurès, J.-P., Devereaux, P., Landais, P., Rubin, D. B., and Thomas, N. (2000). Combining propensity score matching with et al. (2017). Double-adjustment in propensity score matching analysis: additional adjustments for prognostic covariates. J. Am. Stat. Assoc. 95, 573–585. choosing a threshold for considering residual imbalance. BMC Med. Res. doi: 10.1080/01621459.2000.10474233 Methodol. 17, 78. doi: 10.1186/s12874-017-0338-0 Schafer, J. L., and Kang, J. (2008). Average causal effects from nonrandomized Patrick, A. R., Schneeweiss, S., Brookhart, M. A., Glynn, R. J., Rothman, K. J., Avorn, J., studies: a practical guide and simulated example. Psychol. Methods 13, 279. et al. (2011). The implications of propensity score variable selection strategies doi: 10.1037/a0014268 in pharmacoepidemiology: an empirical illustration. Pharmacoepidemiol. Drug Schneeweiss, S. (2006). Sensitivity analysis and external adjustment for Saf. 20, 551–559. doi: 10.1002/pds.2098 unmeasured confounders in epidemiologic database studies of therapeutics. Pearl, J. (2011). Invited commentary: understanding bias amplification. Am. J. Pharmacoepidemiol. Drug Saf. 15, 291–303. doi: 10.1002/pds.1200 Epidemiol. 174, 1223–1227. doi: 10.1093/aje/kwr352 Schneeweiss, S., Gagne, J., Glynn, R., Ruhl, M., and Rassen, J. (2011). Assessing Pearl, J. (2012). On a class of bias-amplifying variables that endanger effect the comparative effectiveness of newly marketed medications: methodological estimates. arXiv e-prints e1203.3503. https://ui.adsabs.harvard.edu/abs/2012​ challenges and implications for drug development. Clin. Pharmacol. Ther. arXiv1203.3503P. 90, 777–790. doi: 10.1038/clpt.2011.235 Peduzzi, P., Concato, J., Feinstein, A. R., and Holford, T. R. (1995). Importance of Schneeweiss, S., Rassen, J. A., Glynn, R. J., Avorn, J., Mogun, H., and Brookhart, events per independent variable in proportional hazards regression analysis ii. M. A. (2009). High-dimensional propensity score adjustment in studies accuracy and precision of regression estimates. J. Clin. Epidemiol. 48, 1503– of treatment effects using health care claims data. Epidemiology 20, 512. 1510. doi: 10.1016/0895-4356(95)00048-8 doi: 10.1097/EDE.0b013e3181a663cc Peduzzi, P., Concato, J., Kemper, E., Holford, T. R., and Feinstein, A. R. (1996). Setoguchi, S., Schneeweiss, S., Brookhart, M. A., Glynn, R. J., and Cook, E. F. (2008). A simulation study of the number of events per variable in logistic regression Evaluating uses of data mining techniques in propensity score estimation: a analysis. J. Clin. Epidemiol. 49, 1373–1379. doi: 10.1016/S0895-4356(96)00236-3 simulation study. Pharmacoepidemiol. Drug Saf. 17, 546–555. doi: 10.1002/ Platt, R. W., Schisterman, E. F., and Cole, S. R. (2009). Time-modified confounding. pds.1555 Am. J. Epidemiol. 170, 687–694. doi: 10.1093/aje/kwp175 Shah, B. R., Laupacis, A., Hux, J. E., and Austin, P. C. (2005). Propensity score Qu, Y., and Lipkovich, I. (2009). Propensity score estimation with missing values methods gave similar results to traditional regression modeling in observational using a multiple imputation missingness pattern (mimp) approach. Stat. Med. studies: a systematic review. J. Clin. Epidemiol. 58, 550–559. doi: 10.1016/j. 28, 1402–1414. doi: 10.1002/sim.3549 jclinepi.2004.10.016 Rassen, J. A., Glynn, R. J., Brookhart, M. A., and Schneeweiss, S. (2011). Sibbald, B., and Roland, M. (1998). Understanding controlled trials. why are randomised Covariate selection in high-dimensional propensity score analyses of controlled trials important? BMJ 316, 201. doi: 10.1136/bmj.316.7126.201 treatment effects in small samples. Am. J. Epidemiol. 173, 1404–1413. doi: Sterne, J. A., White, I. R., Carlin, J. B., Spratt, M., Royston, P., Kenward, M. G., et al. 10.1093/aje/kwr001 (2009). Multiple imputation for missing data in epidemiological and clinical Rassen, J. A., Glynn, R. J., Rothman, K. J., Setoguchi, S., and Schneeweiss, S. (2012). research: potential and pitfalls. BMJ 338, b2393. doi: 10.1136/bmj.b2393 Applying propensity scores estimated in a full cohort to adjust for confounding Stuart, E. A. (2008). Developing practical recommendations for the use of in subgroup analyses. Pharmacoepidemiol. Drug Saf. 21, 697–709. doi: 10.1002/ propensity scores: Discussion of ‘a critical appraisal of propensity score pds.2256 matching in the medical literature between 1996 and 2003’by peter austin, Robins, J. M., Hernan, M. A., and Brumback, B. (2000). Marginal structural statistics in medicine. Stat. Med. 27, 2062–2065. doi: 10.1002/sim.3207 models and causal inference in epidemiology. Epidemiology 11 (5), 550–560 Stuart, E. A. (2010). Matching methods for causal inference: A review and a look doi: 10.1097/00001648-200009000-00011 forward. Stat. Sci. 25, 1. stuart2010matching. doi: 10.1214/09-STS313 Rosenbaum, P. R. (2005). Sensitivity analysis in observational studies. In: Stürmer, T., Joshi, M., Glynn, R. J., Avorn, J., Rothman, K. J., and Schneeweiss, S. Encyclopedia of statistics in behavioral science. Eds. B. S. Everitt and D. C. Howell (2006a). A review of the application of propensity score methods yielded (John Wiley & Sons, Ltd.) 4, 1809–1814. doi: 10.1002/0470013192.bsa606 increasing use, advantages in specific settings, but not substantially different Rosenbaum, P. R., and Rubin, D. B. (1983). The central role of the propensity score estimates compared with conventional multivariable methods. J. Clin. in observational studies for causal effects. Biometrika 70, 41–55. doi: 10.1093/ Epidemiol. 59, 437–4e1. doi: 10.1016/j.jclinepi.2005.07.004 biomet/70.1.41 Stürmer, T., Rothman, K. J., and Glynn, R. J. (2006b). Insights into different results Rosenbaum, P. R., and Rubin, D. B. (1984). Reducing bias in observational studies from different causal contrasts in the presence of effect-measure modification. using subclassification on the propensity score. J. Am. Stat. Assoc. 79, 516–524. Pharmacoepidemiol. Drug Saf. 15, 698–709. doi: 10.1002/pds.1231 doi: 10.1080/01621459.1984.10478078 Stürmer, T., Schneeweiss, S., Avorn, J., and Glynn, R. J. (2005). Adjusting effect Rosenbaum, P. R., and Rubin, D. B. (1985). Constructing a control group using estimates for unmeasured confounding with validation data using propensity multivariate matched sampling methods that incorporate the propensity score. score calibration. Am. J. Epidemiol. 162, 279–289. doi: 10.1093/aje/kwi192 Am. Stat. 39, 33–38. doi: 10.1080/00031305.1985.10479383 Stürmer, T., Schneeweiss, S., Rothman, K. J., Avorn, J., and Glynn, R. J. (2007). Rubin, D. B. (1996). Multiple imputation after 18+ years. J. Am. Stat. Assoc. Performance of propensity score calibration—a simulation study. Am. J. 91, 473–489. doi: 10.1080/01621459.1996.10476908 Epidemiol. 165, 1110–1118. doi: 10.1093/aje/kwm074 Frontiers in Pharmacology | www.frontiersin.org 18 September 2019 | Volume 10 | Article 973 Ali et al. Propensity Score Methods in HTA Tarricone, R., Boscolo, P. R., and Armeni, P. (2016). What type of clinical evidence Westreich, D., Lessler, J., and Funk, M. J. (2010). Propensity score estimation: is needed to assess medical devices? Eur. Respir. Rev. 25, 259–265. doi: 10.1183/​ neural networks, support vector machines, decision trees (cart), and meta- 16000617.0016-2016 classifiers as alternatives to logistic regression. J. Clin. Epidemiol. 63, 826–833. Uddin, M. J., Groenwold, R. H., Ali, M. S., de Boer, A., Roes, K. C., Chowdhury, doi: 10.1016/j.jclinepi.2009.11.020 M. A., et al. (2016). Methods to control for unmeasured confounding in Yuan, H., Ali, M. S., Brouwer, E. S., Girman, C. J., Guo, J. J., Lund, J. L., et al. pharmacoepidemiology: an overview. Int. J. Clin. Pharmacol. Res. 38, 714–723. (2018). Real-world evidence: What it is and what it can tell us according to doi: 10.1007/s11096-016-0299-0 the international society for pharmacoepidemiology (ispe) comparative Vandenbroucke, J. P., Broadbent, A., and Pearce, N. (2016). Causality and causal effectiveness research (cer) special interest group (sig). Clin. Pharmacol. Ther. inference in epidemiology: the need for a pluralistic approach. Int. J. Epidemiol. 104, 239–241. doi: 10.1002/cpt.1086 45, 1776–1786. doi: 10.1093/ije/dyv341 Zhang, X., Faries, D. E., Li, H., Stamey, J. D., and Imbens, G. W. (2018). Von Elm, E., Altman, D. G., Egger, M., Pocock, S. J., Gøtzsche, P. C., Vandenbroucke, Addressing unmeasured confounding in comparative observational research. J. P., et al. (2007). The strengthening the reporting of observational studies Pharmacoepidemiol. Drug Saf. 27, 373–382. doi: 10.1002/pds.4394 in epidemiology (strobe) statement: guidelines for reporting observational studies. PLoS Med. 4, e296. doi: 10.1371/journal.pmed.0040296 Conflict of Interest Statement: LS holds research grants from GSK, Wellcome, Wang, S. V., He, M., Jin, Y., Wyss, R., Shin, H., Ma, Y., et al. (2017). A review of MRC, NIHR, BHF, and Diabetes UK and is a Trustee of the British Heart the performance of different methods for propensity score matched subgroup Foundation. DP-A holds research grants from NIHR, AMGEN, and UCB. analyses and a summary of their application in peer-reviewed research studies. Pharmacoepidemiol. Drug Saf. 26, 1507–1512. doi: 10.1002/pds.4328 The remaining authors declare that the research was conducted in the absence of Wang, S. V., Jin, Y., Fireman, B., Gruber, S., He, M., Wyss, R., et al. (2018). Relative any commercial or financial relationships that could be construed as a potential performance of propensity score matching strategies for subgroup analyses. conflict of interest. Am. J. Epidemiol. 187 (8): 1799–1807. doi: 10.1093/aje/kwy049 Weitzen, S., Lapane, K. L., Toledano, A. Y., Hume, A. L., and Mor, V. (2005). Copyright © 2019 Ali, Prieto-Alhambra, Lopes, Ramos, Bispo, Ichihara, Pescarini, Weaknesses of goodness-of-fit tests for evaluating propensity score models: the Williamson, Fiaccone, Barreto and Smeeth. This is an open-access article distributed case of the omitted confounder. Pharmacoepidemiol. Drug Saf. 14, 227–238. under the terms of the Creative Commons Attribution License (CC BY). The use, doi: 10.1002/pds.986 distribution or reproduction in other forums is permitted, provided the original Westreich, D., Cole, S. R., Funk, M. J., Brookhart, M. A., and Stürmer, T. (2011). author(s) and the copyright owner(s) are credited and that the original publication The role of the c-statistic in variable selection for propensity score models. in this journal is cited, in accordance with accepted academic practice. No use, Pharmacoepidemiol. Drug Saf. 20, 317–320. doi: 10.1002/pds.2074 distribution or reproduction is permitted which does not comply with these terms. Frontiers in Pharmacology | www.frontiersin.org 19 September 2019 | Volume 10 | Article 973