Manual (Version 0.0.3) Updated: 8/28/2023

The package PStrata fits Bayesian principal stratification models using Stan. A wide variety of models, priors and assumptions are allowed to provide flexibility.


Let \(Z\) denote the assigned treatment, \(D\) denote the post-randomization intermediate variable, \(X\) denote the covariates and \(Y\) denote the outcome. Define \(S = (D(0), D(1))\) as the principal stratum.

PStrata requires specification of two models.



Consider fitting a Bayesian principal stratification model on dataset sim_data_normal, with both S-model and Y-model involving only intercepts and no covariates. The Y-model is a linear model (gaussian distribution with link function being the identity function). Assume monotonicity and exclusion restriction on both \(S = (0, 0)\) and \(S = (1, 1)\).

We specify standard normal priors for the intercepts included in the S-model and Y-model. Also, the Y-model introduces a standard deviation parameter \(\sigma\) on which we specify a standard inverse gamma prior.

The following script runs the Bayesian analysis (may take around 2 minutes).

obj <- PStrata(
    S.formula = Z + D ~ 1, # treatment + intermediate ~ covariates
    Y.formula = Y ~ 1,     # outcome ~ covariates = gaussian(link = "identity"),
    data = sim_data_normal,
    strata = c(n = "00*", c = "01", a = "11*"), # * for ER, names are optional
    prior_intercept = prior_normal(0, 1),
    prior_sigma = prior_inv_gamma(1),
    chains = 4, warmup = 500, iter = 1000 # optional parameters to pass to rstan()

Raw output of obj provides an overview of the estimated proportion for each strata and summary(obj) provides quantiles and confidence intervals.

To obtain the estimated mean effects of each principal stratum and treatment arm, run the following script.

res <- PSOutcome(obj)
summary(res, "data.frame") # returns a data.frame object. Other options: "array", "matrix"

To obtain contrasts, use function PSContrast. For example, to compare the stratum-specific treatment effects (defined by \(Y(1) - Y(0)\)), use the following script.

cont <- PSContrast(res, Z = TRUE)
summary(cont, "data.frame")