The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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.


Overview

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.


Assumption


Example

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
    Y.family = 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"
plot(res)

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")
plot(cont)

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.