Toolkit for Weighting and Analysis of Nonequivalent Groups A Tutorial on the Twang Commands for Stata Users

This document and trademark(s) contained herein are protected by law. This representation of RAND intellectual property is provided for noncommercial use only. Unauthorized posting of this publication online is prohibited. Permission is given to duplicate this document for personal use only, as long as it is unaltered and complete. Permission is required from RAND to reproduce, or reuse in another form, any of its research documents for commercial use. For information on reprint and linking permissions, please visit The RAND Corporation is a research organization that develops solutions to public policy challenges to help make communities throughout the world safer and more secure, healthier and more prosperous. RAND is nonprofit, nonpartisan, and committed to the public interest. RAND's publications do not necessarily reflect the opinions of its research clients and sponsors.


Introduction
The Toolkit for Weighting and Analysis of Nonequivalent Groups, TWANG, contains a set of macros to support causal modeling of observational data through the estimation and evaluation of propensity scores and associated weights (Ridgeway et al., 2013). The commands call functions from the twang package in the R environment for statistical computing and graphics (R Core Team, 2013). The twang package was developed in 2004, and after extensive use, it received a major update in 2012. The Stata twang macros were developed in 2015 to support the use of the twang tools without requiring analysts to learn R. This tutorial provides an introduction to twang and demonstrates its use through illustrative examples.
The foundation to the methods supported by twang is the propensity score. The propensity score is the probability that a particular case would be assigned or exposed to a treatment condition. Rosenbaum & Rubin (1983) showed that knowing the propensity score is sufficient to separate the effect of a treatment on an outcome from observed confounding factors that influence both treatment assignment and outcomes, provided the necessary conditions hold. The propensity score has the balancing property that given the propensity score the distribution of features for the treatment cases is the same as that for the control cases. While the treatment selection probabilities are generally not known, good estimates can be effective at diminishing or eliminating confounds between pretreatment group differences and treatment outcomes in the estimation of treatment effects.
There are now numerous propensity score methods in the literature. They differ in how they estimate the propensity score (e.g. logistic regression, CART), the target estimand (e.g. treatment effect on the treated, population treatment effect), and how they utilize the resulting estimated propensity scores (e.g. stratification, matching, weighting, doubly robust estimators). We originally developed the twang package with a particular process in mind, namely, generalized boosted regression to estimate the propensity scores and weighting of the comparison cases to estimate the average treatment effect on the treated (ATT). However, we have updated the package to also handle the case where interest lies in using the population weights (e.g., weighting of comparison and treatment cases to estimate the population average treatment effect, ATE). The main workhorse of twang is the ps command, which implements generalized boosted regression modeling (GBM) to estimate the propensity scores. Other tools are included that allow users to assess the success of the resulting weights at obtaining equivalence (or "balance") in the pretreatment covariate distributions of treatment and control groups. However, the framework and functions of the package are flexible enough to allow the user to use propensity score estimates from other methods and to assess the usefulness of those estimates for ensuring balance between the treatment and control groups using tools from the package. The same set of macros is also useful for other tasks, such as non-response weighting, as discussed in Section 4.
The twang Stata package aims to (i) compute from the data estimates of the propensity scores which yield accurate causal effect estimates, (ii) check the quality of the resulting propensity score weights by assessing whether or not they have the balancing properties that we expect in theory, and (iii) use them in computing treatment effect estimates.
2 An ATT example to start

Set-up
If you have not already done so, you will need to download twang ado files and supporting documents to a folder on your computer. The files are available at http://www.rand.org/statistics/twang/downloads.html. They include: • twang Stata package -files containing the Stata program and help files with details on implementing the commands • Stata_Start-up.pdf -step by step details on installing R • lalonde.dta -Example dataset from the Lalonde Study • lindner.dta -Example dataset from the Lindner Study The datasets and example code will be useful for trying the code presented in this tutorial. Those files are not necessary for you to run your own applications.
To use the ado files, you could place the files in the PERSONAL ado-path directory, which can be identified using adopath. Alternately, you can place the files in any directory, and add the directory to ado-path using the command adopath + "C:\Users\uname\adofile" The help files should also be placed in the same folder. 1 After adding the directory to ado-path, the help files can be opened using help command together with the command that you need help with (e.g., help ps). Note that adopath temporarily adds the directory to the ado-path -you must rerun the command each time Stata is opened.
The ado files will run code in R and import the results into your Stata session. The ado files will export the users' data to a .csv file that can be read into R. They will also create an R script file that is run in R batch mode. The script file exports weights and diagnostic information in .csv files that are then ported back into Stata. All files created by the macro are stored in the directory specified by the user in the objpath option as seen in subsequent sections. Any files in this directory created from previous calls of ps will be overwritten. 1 Users might be tempted to copy and paste code from this PDF document into an editor to run this example code. We advise against this.
Text from the PDF file may not appear the same in a text editor as it does in the PDF file; symbols or spaces may be added. The file "tutorial_code.do" file, that is available with the twang ado files contains all the code from this tutorial in text file. Analysts can use that file to copy the code and run the examples.
To manage file transfer between Stata and R and ensure all the necessary R functions are working, we require the operating system to be Windows Vista or later, Mac OS-X or later, or UNIX.
The ado files rely on the twang package in R. You will need to install R from The Comprehensive R Archive Network (http://cran.us.r-project.org/). The software can be installed by clicking on the link for the users computer platform (e.g., Windows users would click on "Download R for Windows" and then click on the "base" link to download the standard R software). For assistance on installing R please see the Start-up file "Toolkit for Weighting and Analysis of Nonequivalent Groups Stata Commands Start-Up" (Stata_Start-up.pdf) or you can view tutorial videos on Youtube such as https://www.youtube.com/watch?v=PwfVCaMCO8U.
Users will need to note the directory where the R software is installed and the name of the executable file. For Windows users with a 64-bit processor the directory information for the standard installment is C:\Program Files\R\R-3.0.2\bin\x64 where 3.0.2 is replaced by the current version of R at the time of installation. For users with a 32-bit processor the directory information for the standard installment is C:\Program Files\R\R-3.0.2\bin\i386 Again 3.0.2 is replaced by the current version of R at the time of installation. For both 64 or 32 bit processors, the executable is Rscript.exe for batch implementation. Although not necessary, users with some familiarity with R might also want to install the twang package in R. The ado files will install the package if the users do not.
The ado files also create a "TWANG" folder in a standard location based on the operating system. In Windows, the folder is in the user's AppData\Local folder (C:\Users\username\AppData\Local\TWANG would be the default for a user with the "username" as his or her username). In Mac OS-X, the folder is in the user's Library folder (/Users/username/Library/Twang). This folder will remain on the user's hard drive until it is removed. Users can remove the folder using any method they would normally use for removing a folder when they no longer plan to use twang.

Estimating propensity scores with the ps command
To demonstrate the package we utilize data from Lalonde's National Supported Work Demonstration (NSWD) analysis (Lalonde, 1986, Dehejia & Wahba, 1999, http://users.nber.org/ ~rdehejia/nswdata2.html). The NSWD was a temporary employment program that gave work experience and counseling service to disadvantaged workers to help them move into the labor market. A comparison group that did not participate in the NSWD was constructed using the Current Population Survey (CPS) for the same years as the NSWD.
In this example, we will estimate the causal effect of the NSWD on earnings among those who participated in the program, or, in other words, the average treatment effect on the treated (ATT). Pretreatment covariates include age, education, race, ethnicity, education level, marital status, earnings in 1974 (pretreatment), and earnings in 1975 (pretreatment). As we will show, the challenge in this analysis is that the distribution of these pretreatment covariates differs between the individuals who participated in the NSWD and those who did not. The dataset is provided with the twang Stata package (lalonde.dta).
In the lalonde dataset, the variable treat is the 0/1 treatment indicator, 1 indicates "treatment" by being part of the NSWD and 0 indicates "comparison" cases drawn from the CPS. In order to estimate a treatment effect for this demonstration program that is unbiased by pretreatment group differences on other observed covariates, we include the pretreatment covariates listed above in a propensity score model of treatment assignment. The ps command is the primary method in twang for estimating propensity scores. This step is computationally intensive and can take a few minutes.
All the options of the command go after the comma. ntrees, intdepth, and shrinkage are parameters for the GBM that ps computes. The option ntrees is the maximum number of iterations that the GBM will run. There will be a warning if the estimated optimal number of iterations is too close to the bound selected in this option because it indicates that balance may improve if more complex models (i.e., those with more trees or a larger value for ntrees) are considered. The user should increase ntrees or decrease shrinkage and rerun it if this warning appears. The option intdepth controls the level of interactions allowed in the GBM, with larger values specifying more complex models. We specified a value of 2, indicating that the algorithm will consider all two-way interactions between covariates. The GBM estimation algorithm uses shrinkage to enhance the smoothness of resulting model. The shrinkage option controls the amount of shrinkage. Small values such as 0.005 or 0.001 yield smooth fits but require greater values of ntrees to achieve adequate fits. Computational time increases inversely with shrinkage. Additional details on ntrees, intdepth, and shrinkage can be found in McCaffrey, Ridgeway, and Morral (2004).
permtestiters specifies whether p-values for KS statistics are calculated using Monte Carlo methods, which is slow but accurate, or estimated using an analytic approximation that is fast, but produces poor estimates in the presence of many ties. If permtestiters=0 (default) is called, then analytic approximations are used. If permtestiters =500 is called, then 500 Monte Carlo trials are run to establish the reference distribution of KS statistics for each covariate. Higher numbers of trials will produce more precise p-values for the test of the KS statistic. Specifying permtestiters greater than zero can greatly slow down the twang computations. We tend to rely on the approximations (permtestiters =0) when using twang in practice.
The stopmethod option specifies a set (or sets) of rules and measures for assessing the balance, or equivalence, established on the pretreatment covariates of the weighted treatment and control groups. The iterations used in the GBM minimize the differences between the treatment and control groups as measured by the balance statistics specified by values given to stopmethod option. The package includes four built-in stop methods. They are "es.mean", "es.max", "ks.mean", and "ks.max". The four stopping rules are defined by two components: a balance metric for covariates and rule for summarizing across covariates. The balance metric summarizes the difference between two univariate distributions of a single pre-treatment variable (e.g., age). The stopping rules in twang use two balance metrics: absolute standardized bias (also referred to as the absolute standardized mean difference or the Effect Size) and the Kolmogorov-Smirnov (KS) statistic. The stopping rules use two different rules for summarizing across covariates: the mean of the covariate balance metrics ("mean") or the maximum of the individual covariate balance metrics ("max"). The first piece of the stopping rule name identifies the balance metric ("es" for the effect size or standardized bias or "ks" for the KS statistic) and the second piece specifies the method for summarizing across balance metrics ("mean" or "max"). For instance, "es.mean" uses the effect size or the absolute standardized bias and summarizes across variables with the mean and the "ks.max" uses the KS statistics to assess balances and summarizes using the maximum across variables and the other two stopping rules use the remaining two combinations of balance metrics and summary statistics. The balance metrics depend on the estimand and correct specification of the metrics is set automatically by the specification of the estimand option.
The sampw option is the name of the variable that contains sampling weights if they exist. If there are no sampling weights, the parameter can be left unspecified as it is in this example.
The estimand option is used to indicate whether the analyst is interested in estimating the average treatment effect (ATE) or the average treatment effect on the treated (ATT), as we do above. ATE addresses the question of how outcomes would differ if everyone in the sample were given the treatment versus everyone being given the control (Wooldridge, 2002). ATT, on the other hand, estimates the analogous quantity averaging only over the subjects who were actually treated.
The primary results of the ps command are the weights which can be used for estimating effects. It also produces checks of the balance of covariates in the form of balance tables and an overall summary table. It will also produce diagnostic plots to help assess the balance and the GBM fit, if the user requests them by specifying the plotname option.
The propensity score based weights created by the ps command could be found in the two new variables at the end of the original dataset. In this example, we save this new dataset as "lalonde_att_wgts" in the output folder specified by the cd command (save lalonde_att_wgts, replace). There is one weight variable for each stopping rule specified in stopmethod. The weight variables are named according to the stopping rule and estimand so that in this example there is a weight variable 'esmeanatt' with the weights from a GBM with the iterations chosen to minimize the mean standardized bias (effect size) and a second weight variable 'ksmaxatt' with the weights from a GBM with the iterations chosen to minimize the maximum KS statistic. Because the estimand is ATT, the weight equals 1 for every individual in the treatment group. This dataset can be saved and used to estimate the treatment effects.
The options rcmd specifies the R program executable file for running R. The location of the file is determined through the installation of R. The specification of rcmd is not necessary for Mac OS-X, but the default location of the executable is "/usr/bin/Rscript". The default setup of R Version 3.0.2 on Windows 7 resulted in the executable being "C:\Program Files\R\R 3.0.2\bin\x64\Rscript.exe". For other versions of R "3.0.2" will be replaced by the version number. If the analyst has added R to the path environmental variable then the path does not need to be included in the rcmd option, "rcmd(Rscript)" will work. 2 (Similarly specifying the ".exe" extension is not necessary.) The option plotname gives the name for a pdf file of default diagnostic plots that twang creates. Creation of the plots is optional. If plotname is not given, then no plots are created. If the option contains a path, then the pdf file with plots will be stored in the folder specified by it. Otherwise the file will be stored in the folder specified in the option objpath. In our example, we specify plotname as "C:\Users\uname\twang\output\lalonde_twang.pdf" so the plots will be stored in that folder. Users will need to specify an appropriate folder where they can write file to use if they specify it.
The final option objpath specifies a folder where files created by the command to run the twang functions in R and return the results to Stata are stored. Namely, an "R object" ("ps.RData") with the GBM fit information and a log of the R session ("ps.Rout"). The objpath option is required for running the ps command and must reference an existing directory.
Having fit the GBM, the analyst should perform several diagnostic checks before estimating the causal effect in question. The first of these diagnostic checks makes sure that the specified value of ntrees allowed the GBM to explore sufficiently complicated models. We can do this quickly using the "convergence" or "optimization" plot created by the R functions. There are two ways to generate this plot. The first way to obtain the plot is to specify plotname option in ps command. This will create all the default diagnostics plots available in twang. They will be stored in the single pdf file specified by the option. In this example the file is: "C:\Users\uname\twang\output\lalonde_twang.pdf". The default file created by specifying plotname contains one page for each type of diagnostic plot and each page contains a multi-panel plot with one panel for each stopping rule specified in the stopmethod option. If only one stopping rule is specified each page contains a single panel. Figure  1 presents the plots for the first page of the lalonde_twang.pdf file. 2 To add the R directory to the PATH the user can open a command (cmd) window and type: setx PATH "%PATH%;C:\Program Files\R\R--3.0.2\bin\x64 where "C:\Program Files\R\R--3.0.2\bin\x64" should be replaced by actual directory where R is stored. Alternatively, the R directory can be added to the PATH by 1 8. Click OK until no longer prompted and close out windows that were opened. Again C:\Program FilesR\R--3.0.2\bin\x64" should be replaced by actual directory where R is stored. Figure 1: Example of an optimization plot for two stopping rules ("es.mean," and "ks.max") for estimating ATT weights for the Lalonde dataset. It was generated by setting the plotname option in the ps command and appears as the first page of a multiple page document of diagnostics plots.
The second way to obtain the plot is to run the psplot command which can create specific plots and store them using pdf or other file formats. The inputobj option specifies the R object created by the ps command. plotname gives the file name for the plot and plotformat specifies the file format. Allowable file formats are: jpg -JPEG, pdf -PDF, png -PNG, wmf -Windows enhanced metafile, and ps -postscript. The plots option specifies the diagnostic plot to be created. Only one type of diagnostic plot can be created by a psplot command. 3 The plot can be specified by number or name and the names should not be in quotation marks. The convergence plot is specified by "1" or "optimize". The results of the following code produce the same plot as Figure 1, except there are no titles on the plots produced by the psplot command. The psplot command also allows for restricting the plot to results for a single stopping rule by specifying it by number in the subset option. Stopping rules are numbered by alphabetical order; not the order in which they are specified. See Figure 1. The convergence plot plots the balance measures as a function of the number of iterations in the GBM algorithm, with higher iterations corresponding to more complicated fitted models. In this example, 2127 iterations minimized the average effect size difference and 1756 iterations minimized the largest of the eight Kolmogorov-Smirnov (KS) statistics computed for the covariates. This can be observed in Figure 1. The maximum of KS statistics starts large, decreases and then increases somewhere between 1000 and 2000 iterations. The plot suggest ntrees=5000 was sufficient. However, if it had appeared that additional iterations would be likely to result in lower values of the balance statistic -for instance, if the maximum was still declining without appearing to have attained a minimum by the maximum number of iterations, ntrees should be increased and ps rerun. As shown in the plot, after a point, additional complexity typically makes the balance worse. This figure also gives information on how compatible two or more stopping rules are: if the minima for multiple stopping rules under consideration are near one another, the results should not be sensitive to which stopping rule one uses for the final analysis. See Section 5.3 for a discussion of these and other balance measures. Figure 2: Example of an optimization plot for a single stopping rule (ks.max) for estimating ATT weights for the Lalonde dataset.

Assessing "balance" using balance tables
The ps command generates a "balance table" which provides a tabular summary of the balance between the covariate distributions for the treatment and control groups. The table created by the ps command could be found in a csv file "baltab", in the folder specified by objpath, or can be printed by the use of the post-estimation balance command. The weighted tables show how well the resulting weights succeed in manipulating the control group so that its weighted pretreatment characteristics match, or balance, those of the unweighted treatment group if estimand = "ATT" or in adjusting both the control and treatment groups so that their weighted pretreatment characteristics match, or balance, with one another if estimand = "ATE". The unweighted table provides the same information without weighting to give a sense of the general imbalance in the groups. For both the weighted and unweighted tables, balance is assessed separately for each covariate by each stop method. Statistics for each of the specified values to stopmethod are identified by the stop method label with the specified estimand appended, here "es.mean.ATT" and "ks.max.ATT." There is no research on which stopping rule is best and the choice is likely to depend on the application. If there are missing values in the covariates, twang will attempt to construct weights that also balance rates of missingness in the treatment and control groups. In this case, the balance table will have an extra row for each variable that has missing entries. The columns of the table consist of the following items: txmn, ctmn The treatment means and the control means for each of the variables.
The unweighted table shows the unweighted means. For each stopping rule the means are weighted using weights corresponding to the gbm model selected by ps command using the stopping rule. When estimand= "ATT" the weights for the treatment group always equal 1 for all cases and there is no difference between unweighted and propensity score weighted txmn txsd, ctsd The propensity score weighted treatment and control groups' standard deviations for each of the variables. The unweighted table shows the unweighted standard deviations stdeffsz The standardized effect size, defined as the treatment group mean minus the control group mean divided by the treatment group standard deviation if estimand = "ATT" or divided by the pooled sample (treatment and control) standard deviation if estimand = "ATE". (In discussions of propensity scores this value is sometimes referred to as "standardized bias".) Occasionally, lack of treatment group or pooled sample variance on a covariate results in very large (or infinite) standardized effect sizes. For purposes of analyzing mean effect sizes across multiple covariates, we set all standardized effect sizes larger than 500 to NA (missing values) stat, p Depending on whether the variable is continuous or categorical, stat is a t-statistic or a χ 2 statistic. p is the associated p-value ks, kspval The Kolmogorov-Smirnov test statistic and its associated p-value. P-values for the KS statistics are either derived from Monte Carlo simulations or analytic approximations, depending on the specifications made in the permtestiters option of the ps command. For categorical variables this is just the χ 2 test p-value.
The balance command with the summary option returns a compact summary of the sample sizes of the groups and the balance measures in the results window as well as in a csv file "summary" in the same folder specified by objpath. The summary table includes one row for the results using the weights produced by each stopping rule specified by the stopmethod option and one row for the unweighted data. Each row contains: the row name which specifies the stopping rule used for the weights or that the data are unweighted ("unw"), the estimand is also appended to the label; ntreat and nctrl, the treatment and control group sample sizes, respectively; esstreat and essctrl, the effective sample sizes of the treatment and control groups, esstreat equals the treatment group sample size for ATT because the weights are one (additional details on the effective sample size follow); maxes and meanes, and maxks and meanks, the maximum and mean or average of the standardized effect sizes or KS statistics for the covariates, respectively; maxksp, the p-value for testing the maximum of the KS statistics is greater than zero; and iter, the number of iterations or trees in the GBM that minimizes the stopping rule, missing for unw. The maxksp is only produced when permtestiters> 0; otherwise it is missing for all rows of the summary table, as it is for our example. If permtestiters > 0 was used in the call to ps, then Monte Carlo simulation is used to estimate p-values for the maximum KS statistic that would be expected across the covariates, had individuals with the same covariate values been assigned to groups randomly. Thus, a p-value of 0.04 for maxksp indicates that the largest KS statistic found across the covariates is larger than would be expected in 96% of trials in which the same cases were randomly assigned to groups.
In general, weighted means can have greater sampling variance than unweighted means from a sample of equal size. The effective sample size (ESS) of the weighted comparison group captures this increase in variance as where summation is over cases in the control group. The ESS is approximately the number of observations from a simple random sample that yields an estimate with sampling variation equal to the sampling variation obtained with the weighted comparison observations. Therefore, the ESS will give an estimate of the number of comparison participants that are comparable to the treatment group when estimand= "ATT". When the estimand of interest is "ATE", there is an analogous ESS for the treatment group because the weights are no longer equal to one for that group. The ESS is an accurate measure of the relative size of the variance of means when the weights are fixed or they are uncorrelated with outcomes. Otherwise the ESS underestimates the effective sample size (Little & Vartivarian, 2004). With propensity score weights, it is rare that weights are uncorrelated with outcomes. Hence the ESS typically gives a lower bound on the effective sample size, but it still serves as a useful measure for choosing among alternative models and assessing the overall quality of a model, even if it provides a possibly conservative picture of the loss in precision due to weighting.

Graphical assessments of balance
The psplot command can generate useful diagnostic plots to evaluate the propensity scores. The full set of plots available in twang and the option value of plot to produce each one are given in Table 1. The convergence or optimization plot was discussed above. Other diagnostic plots are specified by the value of the plots option 4 . For example, specifying plots(2) or plots(boxplot) produces boxplots illustrating the spread of the estimated propensity scores in the treatment and comparison groups (Figure 3). Whereas propensity score stratification requires considerable overlap in these spreads, excellent covariate balance can often be achieved with weights, even when the propensity scores estimated for the treatment and control groups show little overlap.  The effect size plot (see Figure 4) illustrates the effect of weights on the magnitude of differences between groups on each pretreatment covariate. These magnitudes are standardized using the standardized effect size described earlier. In these plots, substantial reductions in effect sizes are observed for most variables (blue lines), with only one variable showing an increase in effect size (red lines), but only a seemingly trivial increase. Closed red circles indicate a statistically significant difference, many of which occur before weighting, none after. In some analyses variables can have very little variance in the treatment group sample or the entire sample and group differences can be very large relative to the standard deviations. In these situations, the user is warned that some effect sizes are too large to plot. When many of the p-values testing individual covariates for balance are very small, the groups are clearly imbalanced and inconsistent with what we would expect had the groups been formed by random assignment. After weighting we would expect the p-values to be larger if balance had been achieved. We use a QQ plot comparing the quantiles of the observed p-values to the quantiles of the uniform distribution (45 degree line) to conduct this check of balance. Ideally, the p-values from independent tests in which the null hypothesis is true will have a uniform distribution. Although the ideal is unlikely to hold even if we had random assignment (Bland, 2013), severe deviation of the pvalues below the diagonal suggests lack of balance and p-values running at or above the diagonal suggests balance might have been achieved. The p-value plot (plots=4) allows users to visually inspect the p-values of the t-tests for group differences in the covariate means.  Figure 5, presents the t-test p-value plot for the Lalonde example. Before weighting (closed circles), the groups have statistically significant differences on many variables (i.e., p-values are near zero). After weighting (open circles) the p-values are generally above the 45-degree line, which represents the cumulative distribution of a uniform variable on [0,1]. This indicates that the p-values are even larger than would be expected in an ideal randomized study, so that balance is generally good. One can inspect similar plots for the KS statistic with the option plots = 5 or "ks" (see Figure  6).

Analysis of outcomes
The aim of the National Supported Work Demonstration analysis is to determine whether the program was effective at increasing earnings in 1978. We will estimate this effect as the difference in the treatment and weighted control group means and test that it is not zero using a Wald test. The propensity score adjusted test can be computed using regress along with Stata's built in weighting features (svyset followed by the svy: prefix). We start with an analysis using the weights derived from the GBM selected to minimize the mean standardized bias ("es.mean" stopping rule).We use svyset to declare that our weighting variable is esmeanatt. The propensity score adjusted results are then estimated using the svy prefix with regress.  The analysis estimates an increase in earnings of $733 for those that participated in the NSWD compared with similarly situated people observed in the CPS. The effect, however, does not appear to be statistically significant.
Some authors have recommended utilizing both propensity score adjustment and additional covariate adjustment to minimize mean square error or to obtain "doubly robust" estimates of the treatment effect (Huppler-Hullsiek & Louis 2002, Bang & Robins 2005). These estimators are consistent if either the propensity scores are estimated correctly or the regression model is specified correctly. For example, note that the balance table for ks.max.att made the two groups more similar on nodegree, but still some differences remained, 70.8% of the treatment group had no degree while 60.1% of the comparison group had no degree. While linear regression is sensitive to model misspecification when the treatment and comparison groups are dissimilar, the propensity score weighting has made them more similar, perhaps enough so that additional modeling with covariates can adjust for any remaining differences. In addition to potential bias reduction, the inclusion of additional covariates can reduce the standard error of the treatment effect if some of the covariates are strongly related to the outcome. Adjusting for the remaining group difference in the nodegree variable slightly increased the estimate of the program's effect to $920, but the difference is still not statistically significant. We can further adjust for the other covariates, but that too in this case has little effect on the estimated program effect.

Estimating the program effect using linear regression
Users may be wondering whether using twang and weighting to adjust for differences between groups yields different results than the more familiar regression approaches to adjusting for group differences on observed covariates. We now compare our weighted estimates of the program effect to results from a more traditional analysis in which the program effect is estimated by a linear model with a treatment indicator and linear terms for each of the covariates. regress is the standard procedure for fitting such models in Stata.

Propensity scores estimated from logistic regression
Propensity score analysis is intended to avoid problems associated with the misspecification of covariate adjusted models of outcomes, but the quality of the balance and the treatment effect estimates can be sensitive to the method used to estimate the propensity scores. For instance, we consider estimating the propensity scores using logistic regression instead of the ps command and compare the results to the weights from twang. ⋅ ⋅ logit treat age educ black hispan nodegree married re74 re75 ⋅ ⋅ predict phat logit fits the logistic regression model. The default is to model the probability that treat=1. The predict command used after running logit generates the predicted probabilities, which are then saved in a variable named phat. The dataset includes all the variables and appends the predicted probabilities. We can create the ATT weights as shown below with weights for the treatment group equal to one and weights for the control group equal to the odds of treatment.

⋅ ⋅ gen w_logit_att=treat+(1-treat)*phat/(1-phat)
The dxwts command provides the balance assessment tools of twang for weights generated using any method, not just by ps. The options are similar to those in ps except weight variables are now specified. Multiple weights can be assessed but they must all be set for a common estimand. The command produces summary and balance tables that are just like those produced by the ps command except there is no iter variable in the summary because the weights might not come from GBM model.

An ATE example
In the analysis of Section 2, we focused on estimating ATT for the lalonde dataset. In that example, the ATE is not of great substantive interest because not all people who are offered entrance into the program could be expected to take advantage of the opportunity. Further, there is some evidence that the treated subjects were drawn from a subset of the covariate space. In particular, in an ATE analysis, we see that we are unable to achieve balance, especially for the black indicator.
We now turn to an ATE analysis that is feasible and meaningful. We focus on the lindner dataset, which was included in the USPS package in R (Obenchain 2011), and is included with the twang Stata package as lindner.dta. A tutorial by Helmreich and Pruzek (2009;HP) for the PSAgraphics package also uses propensity scores to analyze a portion of these data. HP describe the data as follows on p. 3 with our minor recodings in square braces: The lindner data contain data on 996 patients treated at the Lindner Center, Christ Hospital, Cincinnati in 1997. Patients received a Percutaneous Coronary Intervention (PCI). The data consists of 10 variables. Two are outcomes: The question is whether this association is causal. If health care policies were to be made on the basis of these data, we would wish to elicit expert opinion as to whether there are likely to be other confounding pretreatment variables. For this tutorial, we simply follow HP in choosing the pretreatment covariates. The twang model is fit as follows ⋅ ⋅ ps abcix stent height female diabetic acutemi ejecfrac ves1proc, /// stopmethod(es.mean ks.mean) estimand(ATE) /// rcmd(C:\Program Files\R\R-3.0.3\bin\Rscript.exe) /// objpath(C:\Users\uname\twang\output) /// plotname(C:\Users\uname\twang\output\abcix_twang.pdf) We set estimand = "ATE" because we are interested in the effects of abciximab on everyone in the population. We specify the stopping rules to be es.mean and ks.mean. We then inspect preand post-weighting balance using the balance table.      This balance table shows that stent, acutemi, and ves1proc were all significantly imbalanced before weighting. After weighting (using either stop.method considered) we do not see problems in this regard. Examining the diagnostic plots created by the specification of the plotname does not reveal problems, either. In regard to the optimize plot, we note that the scales of the KS and ES statistics presented in the optimize plots are not necessarily comparable. The fact that the KS values are lower than the ES values in the optimize plot does not suggest that the KS stopping rule is finding superior models. Each panel of the optimize plot indicates the gbm that minimizes each stopping rule. The panels should not be compared other than to compare the number of iterations selected by each rule.   From the summary table generated by ps, we see that the "es_mean_ATE" stopping rule results in a slightly higher ESS with comparable balance measures, so we proceed with those weights. Also, we note that esstreat is no longer equal to ntreat since we are focusing on ATE rather than ATT.   The reweighting does not diminish the association between the treatment and the outcome. Indeed, it is marginally more significant after the reweighting. Alternatively, one could use logistic regression to assess the relationship between the dichotomous outcome and dichotomous treatment.

Non-response weights
The twang commands were designed to estimate propensity score weights for the evaluation of treatment effects in observational or quasi-experimental studies. However, we find that the package includes functions and diagnostic tools that are highly valuable for other applications, such as for generating and diagnosing nonresponse weights for survey nonresponse or study attrition. We now present an example that uses the tools in twang. This example uses the subset of the US Sustaining Effects Study data distributed with the HLM software (Bryk, Raudenbush, Congdon, 1996), also available in the R package mlmRev, and included with the twang Stata package as egsingle.dta. The data include mathematics test scores for 1,721 students in kindergarten to fourth grade. They also include student race (black, Hispanic, or other), gender, an indicator for whether or not the student had been retained in grade, the percent low income students at the school, the school size, the percent of mobile students, the students' grade-levels, student and school IDs, and grades converted to year by centering. The study analysis plans to analyze growth in math achievement from grade 1 to grade 4 using only students with complete data. However, the students with complete data differ from other students. To reduce bias that could potentially result from excluding incomplete cases, our analysis plan is to weight complete cases with nonresponse weights.
The goal of nonresponse weighting is to develop weights for the respondents that make them look like the entire sample --both the respondents and nonrespondents. Since the respondents already look like themselves, the hard part is to figure out how well each respondent represents the nonrespondents. Nonresponse weights equal the reciprocal of the probability of response and are applied only to respondents.
Note that the probability of response is equivalent to the propensity score if we consider subjects with an observed outcome to be the "treated" group, and those with an unobserved outcome to be the "controls". We wish to reweight the sample to make it equivalent to the population from which the sample was drawn, so ATE weights are appropriate in this case. Further, recall that the weights for the treated subjects are 1/p in an ATE analysis. Therefore we can reweight the sample of respondents using the weights returned from the ps command.
Before we can generate nonresponse weights, we need to prepare the data using the following commands. The data contain zero, one or two observations for students from grade "0" (kindergarten) and zero or one observation for each of grades 1 to 4. Only students with data from each of grades 1 to 4 will be included so we need to identify those students. There are 1,721 children in the study and 843 (49%) have the necessary four years of outcome data. As discussed above, to use ps to estimate nonresponse, we let respondents be the treatment group by modeling an indicator of response.             By default the balance table generated by ps and balance compares the weighted treatment group (respondents) to the weighted comparison group (nonresponders) --both groups weighted to equal the overall population. However, the goal is to weight the respondents to match the population not to compare the weighted respondents and nonrespondents. The default balance table may be useful for evaluating the propensity scores, but it does not directly assess the quality of the weights for balancing the weighted respondents with the overall population.
We can "trick" the dxwts command in twang into making the desired comparison. We want to compare the weighted respondents to the unweighted full sample. When evaluating ATT weights, we compare the weighted comparison group with the unweighted treatment group. If we apply dxwts to a data set where the "treatment" group is the entire esingle sample and the "control" group is the esingle respondents and the weights equal one for every student in the pseudo-treatment group and equal the weights from ps for every student in the pseudo-control group, we can obtain the balance statistics we want. We begin by setting up the data with the pseudo-treatment and control groups. We add ATE weights from the "ks.max" stopping rule as our nonresponse weights. We now stack the full sample and the respondents. The variable "nr2" is the pseudo-treatment indicator. We set it equal to one for the full sample and 0 for the respondents. Similarly, "wgt2" is the pseudo-ATT weight which is set equal to one for the full sample and equal to the nonresponse weights for the respondents.

⋅ ⋅ use _inputds, clear
⋅ ⋅ use egsingle_nrespwt, clear ⋅ ⋅ append using egsingle, gen(nr2) ⋅ ⋅ gen wgt2=1 if nr2==1 ⋅ ⋅ replace wgt2=wgt if nr2==0 We now run dxwts to obtain the balance statistics.     The resulting balance table includes a table for an unweighted comparison of the respondents with the overall sample and a weighted comparison. We reproduce only the weighted comparison here. In the table, the columns for the treatment group mean and standard deviation ("tx.mn" and "tx.sd") contain the sample statistics for the full sample (egsingle) and the columns for the comparison group ("ct.mn" and "ct.sd") contain the weighted respondent. The following codes prepares an analysis file with all of the data from the respondents with the nonresponse weights included.

Propensity scores and weighting
Propensity scores can be used to reweight comparison cases so that the distribution of their features match the distribution of features of the treatment cases, for ATT, or cases from both treatment and control groups to match each other, for ATE (Rosenbaum 1987, Wooldridge 2002, Hirano and Imbens 2001, McCaffrey et al. 2004) Let f(x|t = 1) be the distribution of features for the treatment cases and f(x|t = 0) be the distribution of features for the comparison cases. If treatments were randomized then we would expect these two distributions to be similar. When they differ for ATT we will construct a weight, w(x), so that For example, if f(age=65, sex=F|t = 1) = 0.10 and f(age=65; sex=F|t = 0) = 0.05 (i.e. 10% of the treatment cases and 5% of the comparison cases are 65 year old females) then we need to give a weight of 2.0 to every 65 year old female in the comparison group so that they have the same representation as in the treatment group. More generally, we can solve (2) for w(x) and apply Bayes Theorem to the numerator and the denominator to give an expression for the propensity score weight for comparison cases, where K is a normalization constant that will cancel out in the outcomes analysis. Equation (3) indicates that if we assign a weight to comparison case i equal to the odds that a case with features x i would be exposed to the treatment, then the distribution of their features would balance. Note that for comparison cases with features that are atypical of treatment cases, the propensity score P(t = 1|x) would be near 0 and would produce a weight near 0. On the other hand, comparison cases with features typical of the treatment cases would receive larger weights. For ATE, each group is weighted to match the population. The weight must satisfy: f(x|t = 1) = w(x)f(x); and (4) f(x|t = 0) = w(x)f(x); and Again using Bayes Theorem we obtain w(x) ∝1/f(t = 1|x) for the treatment group and w(x) ∝ 1/f(t =0|x) for the control group.

Estimating the propensity score
In randomized studies P(t = 1|x) is known and fixed in the study design. In observational studies the propensity score is unknown and must be estimated, but poor estimation of the propensity scores can cause just as much of a problem for estimating treatment effects as poor regression modeling of the outcome. Linear logistic regression is the common method for estimating propensity scores, and can suffice for many problems. Linear logistic regression for propensity scores estimates the log-odds of a case being in the treatment given x as Usually, β is selected to maximize the logistic log-likelihood Maximizing (7) provides the maximum likelihood estimates of β. However, in an attempt to remove as much confounding as possible, observational studies often record data on a large number of potential confounders, many of which can be correlated with one another. Standard methods for fitting logistic regression models to such data with the iteratively reweighted least squares algorithm can be statistically and numerically unstable. To improve the propensity score estimates we might also wish to include non-linear effects and interactions in x. The inclusion of such terms only increases the instability of the models. One increasingly popular method for fitting models with numerous correlated variables is the lasso (least absolute subset selection and shrinkage operator) introduced in statistics in Tibshirani (1996). For logistic regression, lasso estimation replaces (7) with a version that penalizes the absolute magnitude of the coefficients The second term on the right-hand side of the equation is the penalty term since it decreases the overall of when there are coefficients that are large in absolute value. Setting = 0 returns the standard (and potentially unstable) logistic regression estimates of . Setting to be very large essentially forces all of the ! to be equal to 0 (the penalty excludes ! ). For a fixed value of the estimated can have many coefficients exactly equal to 0, not just extremely small but precisely 0, and only the most powerful predictors of t will be non-zero. As a result the absolute penalty operates as a variable selection penalty. In practice, if we have several predictors of t that are highly correlated with each other, the lasso tends to include all of them in the model, shrink their coefficients toward 0, and produce a predictive model that utilizes all of the information in the covariates, producing a model with greater out-of-sample predictive performance than models fit using variable subset selection methods.
Our aim is to include as covariates all piecewise constant functions of the potential confounders and their interactions. That is, in x we will include indicator functions for continuous variables like I(age <15); I(age < 16),…, I(age < 90), etc., for categorical variables like I(sex = male), I(prior MI = TRUE), and interactions among them like I(age < 16)I(sex = male)I(prior MI = TRUE). This collection of basis functions spans a plausible set of propensity score functions, are computationally efficient, and are at the extremes of x reducing the likelihood of propensity score estimates near 0 and 1 that can occur with linear basis functions of x. Theoretically with the lasso we can estimate the model in (8), selecting a small enough so that it will eliminate most of the irrelevant terms and yield a sparse model with only the most important main effects and interactions. Boosting (Friedman 2001, 2003, Ridgeway 1999) effectively implements this strategy using a computationally efficient method that Efron et al. (2004) showed is equivalent to optimizing (8). With boosting it is possible to maximize (8) for a range of values of with no additional computational effort than for a specific value of . We use boosted logistic regression as implemented in the generalized boosted modeling (gbm) package in R (Ridgeway 2005).

Evaluating the weights
As with regression analyses, propensity score methods cannot adjust for unmeasured covariates that are uncorrelated with the observed covariates. Nonetheless, the quality of the adjustment for the observed covariates achieved by propensity score weighting is easy to evaluate. The estimated propensity score weights should equalize the distributions of the cases' features as in (2). This implies that weighted statistics of the covariates of the comparison group should equal the same statistics for the treatment group. For example, the weighted average of the age of comparison cases should equal the average age of the treatment cases. To assess the quality of the propensity score weights one could compare a variety of statistics such as means, medians, variances, and Kolmogorov-Smirnov statistics for each covariate as well as interactions. The twang commands provide both the standardized effect sizes and KS statistics and p-values testing for differences in the means and distributions of the covariates for analysts to use in assessing balance.

Analysis of outcomes
With propensity score analyses the final outcomes analysis is generally straightforward, while the propensity score estimation may require complex modeling. Once we have weights that equalize the distribution of features of treatment and control cases by reweighting. For ATT, we give each treatment case a weight of 1 and each comparison case a weight w i = p(x i )/(1 -p(x i )). To estimate the ATE, we give control cases weight w i = 1/p(x i ) and we give the treatment cases w i = 1/(1-p(x i )). We then estimate the treatment effect estimate with a weighted regression model that contains only a treatment indicator. No additional covariates are needed if the weights account for differences in x.
A combination of propensity score weighting and covariate adjustment can be useful for several reasons. First, the propensity scores may not have been able to completely balance all of the covariates. The inclusion of these covariates in addition to the treatment indicator in a weighted regression model may correct this if the imbalance is relatively small. Second, in addition to exposure, the relationship between some of the covariates and the outcome may also be of interest. Their inclusion can provide coefficients that can estimate the direction and magnitude of the relationship. Third, as with randomized trials, stratifying on covariates that are highly correlated with the outcome can improve the precision of estimates. Lastly, the some treatment effect estimators that utilize an outcomes regression model and propensity scores are "doubly robust" in the sense that if either the propensity score model is correct or the regression model is correct then the treatment effect estimator will be unbiased (Bang & Robins 2005).