*! version 1.1.0 28Oct2016

*! version 1.0.0 21Aug2016
help pstrata

Title

pstrata -- Optimal propensity score stratification

Syntax

pstrata treat [if] [in]
,
pscore(varlist)
[ smin(#)
smax(#)
plevel(#)
common
replace
prefix(string)
]

treat must contain integer values representing the treatment levels

pscore(varlist) required. One or more propensity scores must be provided, depending on the number of treatment levels

smin(#) minimum number of quantiles to start with; default is 5

smax(#) maximum number of quantiles to be tested; default is 50

plevel(#) significance level of the test for assessing balance on the propensity score; default is 0.05

common use only those observations within the region of common support when generating quantiles

replace replace the strata variables created by pstrata if they already exist

prefix adds a prefix to the names of the strata variables created by pstrata; default name is strata1, strata2, etc.

Description

pstrata stratifies the propensity score into an optimal number of quantiles (meaning, strictly, quantile-based bins), where optimal means that no statistical
differences are found between treatment groups within any quantile of the propensity score. In general, stratification (also referred to as subclassification) allows the
investigator to analyze outcomes between treatment groups within each quantile as well as to observe overall differences between groups across all quantiles (Linden & Adams 2008).
It has been shown that stratification of the propensity score into five quantiles can remove over 90% of the initial bias due to the covariates used to create
the propensity score (Cochran 1968; Rosenbaum & Rubin 1984). While the strata generated by pstrata can be used in any subsequent analysis, it was specifically intended
for use in conjunction with mmws; which can be downloaded from SSC.

Before using pstrata, the propensity score for a given treatment must be estimated. A logit or probit regression model can be used for estimating the propensity score
for a binary treatment and a multinomial logistic or probit regression can be used for estimating propensity scores for multiple treatment levels. pstrata then generates
one or more strata variables (named strata1, strata2, etc.), depending on the number of treatment groups in the study.

Remarks

pstrata implements an iterative process to generate the optimal number of quantiles of the propensity score. xtile initially stratifies the propensity
score into 5 quantiles (or any other user-specified minimum number, using smin()). The propensity score is then regressed on the treatment variable
within each quantile to assess whether the propensity score is balanced between the treatment groups under study (i.e. treatment groups are not statistically different,
based on plevel()). If the propensity score is imbalanced within any quantile, xtile is reissued adding one additional quantile, and balance is again assessed.
This procedure continues until one of the following conditions are met:

(1) the propensity score is balanced within all quantiles or,

(2) there is at least one treatment group within a given stratum with zero observations or,

(3) the maximum number of quantiles has been reached (i.e. smax()) but balance is not achieved or,

(4) there is an insufficient number of observations within any given quantile to statistically test balance.

If the first condition is met, then the program ends with the optimal number of quantiles being determined (that is, the minimum number of quantiles has been found with
no statistically significant differences on the propensity score between treatment groups). If any of conditions 2-4 are met, the program terminates and an error
message is displayed indicating in which strata the error was identified. When the study involves multiple treatment levels (groups), this process continues iteratively
across each respective propensity score. With multiple treatments, it is imperative that balance be achieved across all propensity scores, not only a subset. Thus pstrata
will error out if conditions 2-4 are met within any of the propensity scores.

When pstrata errors out due to conditions 2-4, the user should consider re-estimating the propensity score(s), collapsing one or more treatment arms, and/or experimenting
with the various command options, until a viable solution is achieved.

It is important to note that balance on the propensity score does not ensure that balance has been achieved on the underlying covariates. Thus, after implementing pstrata,
balance should be tested on the covariates underlying each propensity score, and within each strata using covbal; downloadable from SSC. If
covariate balance has not been achieved, or if balance on the propensity score itself has not achieved, the propensity score(s) should be re-estimated using additional
covariates (if available) as well as interactions between covariates and polynomials of continuous variables. At least in theory, a well estimated propensity score should accurately
represent the underlying covariates.

Options

pscore(varlist) required. For binary or ordinal level treatments, one propensity score must
be provided. For multiple nominal level treatments, one propensity score must be provided for each treatment level and they must be
ordered correspondingly (e.g. ps1 ps2 ps3 correspond to treatment levels 0 1 2, respectively).

smin(#) indicates the number of quantiles that pstrata will start with; the default is 5.

smax(#) indicates the maximum number of quantiles that pstrata will test; the default is 50.

plevel(#) indicates the significance level of the ANOVA model used for testing balance on the propensity score; the default is 0.05.

common restricts generation of quantiles to those observations within the region of common support.

replace replaces the strata variables created by pstrata if they already exist. If prefix() is specified, only the strata variables created by pstrata with
the same prefix will be replaced.

prefix(string) adds a prefix to the names of the strata variables created by pstrata. Short prefixes are recommended.

Examples

Load example data

. webuse cattaneo2, clear

(1) Binary treatment

Estimate propensity score for the binary treatment mbsmoke and use pstrata to generate an optimal number of quantiles with the default settings

. logit mbsmoke mmarried c.mage##c.mage fbaby medu

. predict pscore, pr

. pstrata mbsmoke, ps(pscore)

Review the p-values estimated for each stratum

. matrix list r(pval1)

Use the strata variable generated by pstrata to generate mmws weights

. mmws mbsmoke, pscore(pscore) strata(strata1)

Examine covariate balance, within strata, with covbal using the weights generated in the previous run

. bys strata1: covbal mbsmoke pscore mmarried mage fbaby medu, abs for(%9.3f) wt( _mmws)

Estimate average treatment effects of mbsmoke on bweight using MMWS weights generated in the previous run

. regress bweight mbsmoke [pw=_mmws], robust

(3) Multiple treatments

Load example data

. webuse cattaneo2, clear

Estimate propensity scores for the four treatment levels of msmoke and generate four strata variables corresponding to each treatment, respectively.

. mlogit msmoke mage medu fage fedu c.mage#c.mage c.mage#c.medu c.mage#c.fage c.mage#c.fedu c.medu#c.medu, base(0)

. predict double (ps1 ps2 ps3 ps4), pr

. pstrata msmoke, ps(ps1 ps2 ps3 ps4)

Generate weights for nominal treatments using the strata variable created in the previous step

. mmws msmoke, pscore(ps1 ps2 ps3 ps4) nominal strata(strata1 strata2 strata3 strata4)

Estimate average treatment effect of each level of msmoke on bweight using MMWS weights generated in the previous run

. regress bweight i.msmoke [pw=_mmws], robust

Run margins to get reweighted mean bweight for each level of msmoke

. margins msmoke

Run marginsplot to visually display reweighted mean bweight for each level of msmoke

. marginsplot, plotopts(connect(i))

Run pairwise comparisons between reweighted treatment levels, with Bonferroni adjustment for multiple tests

. margins msmoke, pwcompare(effects) mcompare(bonferroni)

Saved results

By default, pstrata returns the following results, which
can be displayed by typing return list after
pstrata is finished (see return).

Scalars

r(suppmin) the minimum value of common support (will have a suffix for multiple treatments)

r(suppmax) the maximum value of common support (will have a suffix for multiple treatments)

Matrices

r(pval) the p-values for each quantile (will have a suffix for each respective propensity score, e.g. pval1)

References

Cochran, W. G. 1968. The effectiveness of adjustment by subclassification
in removing bias in observational studies. Biometrics 24: 205-213.

Linden, A. 2014. Combining propensity score-based stratification and weighting to improve
causal inference in the evaluation of health care interventions. Journal of Evaluation in Clinical Practice 20: 1065-1071.

Linden, A. & Adams, J. L. 2008. Improving participant selection in
disease management programs: insights gained from propensity score
stratification. Journal of Evaluation in Clinical Practice 14, 914-918.

Linden, A., Uysal, S. D., Ryan, A., & Adams, J. L. 2016. Estimating causal effects for multivalued treatments:
A comparison of approaches. Statistics in Medicine 35: 534-552.

Rosenbaum, P. R. & Rubin, D. B. 1984. Reducing bias in observational
studies using subclassification on the propensity score. Journal of the American Statistical Association 79, 516-524.

Citation of pstrata

pstrata is not an official Stata command. It is a free contribution
to the research community, like a paper. Please cite it as such:

Linden, A. 2016. pstrata: Stata module for implementing optimal propensity score stratification. http://ideas.repec.org/c/boc/bocode/s458232.html

Author

Ariel Linden

President, Linden Consulting Group, LLC

Ann Arbor, MI, USA

alinden@lindenconsulting.org

http://www.lindenconsulting.org

Acknowledgments

I wish to thank Nicholas J. Cox for his support while developing pstrata.

Also see

Online: xtile, logit, mlogit, margins, marginsplot, teffects,
mmws (if installed), covbal (if installed)