rpietro/twang_examples.r

## twang_examples.r
#this documented script is based on the article "Toolkit for Weighting and Analysis of Nonequivalent Groups: A tutorial for the twang package" http://goo.gl/2xYBX as well as the package documentation at http://goo.gl/BWaLH

library(twang)
library(lattice)
set.seed(1)
data(lalonde)
ps.lalonde <- ps(treat ~ age + educ + black + hispan + nodegree + married + re74 + re75, data = lalonde, n.trees=5000, #n.trees is the maximum number of iterations that gbm will run
                 interaction.depth=2, #The maximum depth of variable interactions. 1 implies an additive model, 2 implies a model with up to 2-way interactions, etc.
                 shrinkage=0.01, #ps() will issue a warning if the estimated optimal number of iterations is too close to the bound selected in this argument because it indicates that balance may improve if more complex models (i.e., those with more trees) are considered. The user should increase n.trees or decrease shrinkage if this warning appears.
                 perm.test.iters=0, #specifies whether p-values for KS statistics should be calculated using Monte Carlo methods, which is slow but can be accurate, or estimated using an analytic approximation that is fast, but produces poor estimates in the presence of many ties. If perm.test.iters=0 is called, then analytic approximations are used. If perm.test.iters=500 is called, then 500 Monte Carlo trials are run to establish the reference distribution of KS statis- tics for each covariate. Higher numbers of trials will produce more precise p-values.
                 stop.method=c("es.mean","ks.max"), #The stop.method argument specifies a set (or sets) of rules and measures for assessing the balance, or equivalence, established on the pretreatment covariates of the treatment and weighted control group. The ps function selects the optimal number of GBM iterations to min- imize the differences between the treatment and control groups as measured by the rules of the given stop.method object. The package includes four built-in stop.method objects. 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 default stopping rules in twang use two balance metrics: absolute standardized bias (also referred to as the absolute standardized mean difference of the Effect Size) and the Kolmogorov-Smirnov (KS) statistic. The stopping rule use two different rules for summarizing across covariates: the mean of the covariate balance metrics (“mean”) or the maximum of the balance metrics (“max”). The first piece of the stopping rule name identifies the balance metric (ES or KS) and the second piece specifies the method for summarizing across balance metrics. 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 re- maining two combinations of balance metrics and summary statistics. The variable distributions used in the balance metrics depend on whether we are interested in estimating the ATT or ATE, and correct specification of these distributions is set automatically by the specification of the estimand in the ps() function.
                 estimand = "ATT", #The estimand argument 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 estimand argument was added to the 2012 revision of the package which integrated ATE weighting into the package and the ps function estimate of the propensity score.
                 verbose=FALSE) #if TRUE, lots of information will be printed to monitor the the progress of the fitting

#------------------------------------
#Checking whether the specified value of n.trees allowed GBM to explore sufficiently complicated models
#------------------------------------

plot(ps.lalonde) #If it appears that additional iterations would be likely to result in lower values of the balance statistic, n.trees should be increased. However, after a point, additional complexity typically makes the balance worse, as in the example below. 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.
plot(ps.lalonde, subset = 2) #If we wish to focus on only one stopping rule, the plotting commands also take a subset argument.
summary(ps.lalonde$gbm.obj, n.trees=ps.lalonde$desc$ks.max.ATT$n.trees, plot=FALSE) #computes the relative influence of each variable for estimating the probability of treatment assignment. optimal number for minimizing the largest of the KS statistics. This value can be found in the ps.lalonde$desc$ks.max.ATT$n.trees.
summary(ps.lalonde$gbm.obj, n.trees=ps.lalonde$desc$ks.max.ATT$n.trees, plot=TRUE) #barchart of the relative influence

#----------------------------------------
#Assessing “balance” using balance tables
#----------------------------------------
lalonde.balance <- bal.table(ps.lalonde)
lalonde.balance #returns information on the pretreatment covariates before and after weighting. The object is a list with named components, one for an unweighted analysis (named unw) and one for each stop.method specied, here es.mean and ks.max. McCafrey et al (2004) essentially used es.mean for the analyses, but our more recent work has sometimes used ks.max. tx.mn, ct.mn The treatment means and the control means for each of the variables. Theunweighted table (unw) shows the unweighted means. For each stopping rule the means are weighted using weights corresponding to the gbm model selected by ps() 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 tx.mn. tx.sd, ct.sd The propensity score weighted treatment and control groups' standard deviations for each of the variables. The unweighted table (unw) shows the unweighted standard deviations std.eff.sz The standardized efect 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 \standard- ized 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, ks.pval 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 specifcations made in the perm.test.iters argument of the ps function.  For categorical variables this is just the  2 test p-value

library(xtable)
pretty.tab <- lalonde.balance$ks.max.ATT[,c("tx.mn","ct.mn","ks")]
pretty.tab <- cbind(pretty.tab, lalonde.balance$unw[,"ct.mn"])
names(pretty.tab) <- c("E(Y1|t=1)","E(Y0|t=1)","KS","E(Y0|t=0)")
xtable(pretty.tab, caption = "Balance of the treatment and comparison groups", label = "tab:balance", digits = c(0, 2, 2, 2, 2), align=c("l","r","r","r","r"))

summary(ps.lalonde)

plot(ps.lalonde, plots=2)
plot(ps.lalonde, plots=3) #generates an error
plot(ps.lalonde, plots = 4)
	#this documented script is based on the article "Toolkit for Weighting and Analysis of Nonequivalent Groups: A tutorial for the twang package" http://goo.gl/2xYBX as well as the package documentation at http://goo.gl/BWaLH

	library(twang)
	library(lattice)
	set.seed(1)
	data(lalonde)
	ps.lalonde <- ps(treat ~ age + educ + black + hispan + nodegree + married + re74 + re75, data = lalonde, n.trees=5000, #n.trees is the maximum number of iterations that gbm will run
	interaction.depth=2, #The maximum depth of variable interactions. 1 implies an additive model, 2 implies a model with up to 2-way interactions, etc.
	shrinkage=0.01, #ps() will issue a warning if the estimated optimal number of iterations is too close to the bound selected in this argument because it indicates that balance may improve if more complex models (i.e., those with more trees) are considered. The user should increase n.trees or decrease shrinkage if this warning appears.
	perm.test.iters=0, #specifies whether p-values for KS statistics should be calculated using Monte Carlo methods, which is slow but can be accurate, or estimated using an analytic approximation that is fast, but produces poor estimates in the presence of many ties. If perm.test.iters=0 is called, then analytic approximations are used. If perm.test.iters=500 is called, then 500 Monte Carlo trials are run to establish the reference distribution of KS statis- tics for each covariate. Higher numbers of trials will produce more precise p-values.
	stop.method=c("es.mean","ks.max"), #The stop.method argument specifies a set (or sets) of rules and measures for assessing the balance, or equivalence, established on the pretreatment covariates of the treatment and weighted control group. The ps function selects the optimal number of GBM iterations to min- imize the differences between the treatment and control groups as measured by the rules of the given stop.method object. The package includes four built-in stop.method objects. 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 default stopping rules in twang use two balance metrics: absolute standardized bias (also referred to as the absolute standardized mean difference of the Effect Size) and the Kolmogorov-Smirnov (KS) statistic. The stopping rule use two different rules for summarizing across covariates: the mean of the covariate balance metrics (“mean”) or the maximum of the balance metrics (“max”). The first piece of the stopping rule name identifies the balance metric (ES or KS) and the second piece specifies the method for summarizing across balance metrics. 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 re- maining two combinations of balance metrics and summary statistics. The variable distributions used in the balance metrics depend on whether we are interested in estimating the ATT or ATE, and correct specification of these distributions is set automatically by the specification of the estimand in the ps() function.
	estimand = "ATT", #The estimand argument 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 estimand argument was added to the 2012 revision of the package which integrated ATE weighting into the package and the ps function estimate of the propensity score.
	verbose=FALSE) #if TRUE, lots of information will be printed to monitor the the progress of the fitting

	#------------------------------------
	#Checking whether the specified value of n.trees allowed GBM to explore sufficiently complicated models
	#------------------------------------

	plot(ps.lalonde) #If it appears that additional iterations would be likely to result in lower values of the balance statistic, n.trees should be increased. However, after a point, additional complexity typically makes the balance worse, as in the example below. 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.
	plot(ps.lalonde, subset = 2) #If we wish to focus on only one stopping rule, the plotting commands also take a subset argument.
	summary(ps.lalonde$gbm.obj, n.trees=ps.lalonde$desc$ks.max.ATT$n.trees, plot=FALSE) #computes the relative influence of each variable for estimating the probability of treatment assignment. optimal number for minimizing the largest of the KS statistics. This value can be found in the ps.lalonde$desc$ks.max.ATT$n.trees.
	summary(ps.lalonde$gbm.obj, n.trees=ps.lalonde$desc$ks.max.ATT$n.trees, plot=TRUE) #barchart of the relative influence

	#----------------------------------------
	#Assessing “balance” using balance tables
	#----------------------------------------
	lalonde.balance <- bal.table(ps.lalonde)
	lalonde.balance #returns information on the pretreatment covariates before and after weighting. The object is a list with named components, one for an unweighted analysis (named unw) and one for each stop.method specied, here es.mean and ks.max. McCafrey et al (2004) essentially used es.mean for the analyses, but our more recent work has sometimes used ks.max. tx.mn, ct.mn The treatment means and the control means for each of the variables. Theunweighted table (unw) shows the unweighted means. For each stopping rule the means are weighted using weights corresponding to the gbm model selected by ps() 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 tx.mn. tx.sd, ct.sd The propensity score weighted treatment and control groups' standard deviations for each of the variables. The unweighted table (unw) shows the unweighted standard deviations std.eff.sz The standardized efect 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 \standard- ized 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, ks.pval 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 specifcations made in the perm.test.iters argument of the ps function. For categorical variables this is just the 2 test p-value

	library(xtable)
	pretty.tab <- lalonde.balance$ks.max.ATT[,c("tx.mn","ct.mn","ks")]
	pretty.tab <- cbind(pretty.tab, lalonde.balance$unw[,"ct.mn"])
	names(pretty.tab) <- c("E(Y1\|t=1)","E(Y0\|t=1)","KS","E(Y0\|t=0)")
	xtable(pretty.tab, caption = "Balance of the treatment and comparison groups", label = "tab:balance", digits = c(0, 2, 2, 2, 2), align=c("l","r","r","r","r"))

	summary(ps.lalonde)

	plot(ps.lalonde, plots=2)
	plot(ps.lalonde, plots=3) #generates an error
	plot(ps.lalonde, plots = 4)