3 Basic functionality

This page illustrates the basic functionality of the pstratreg function. This function conducts parametric principal stratification analysis to estimate average causal effect among the always-valid subgroup whose outcome would exist in either treatment condition.

Jargon? Start with the first page on the goal!

The package automates the process to

estimate a mediator model
estimate an outcome model
- allowing heteroskedasticity if needed
calculate the conditional probability of being always-valid
implement monotonicity assumptions
bound estimates using the conditional outcome distribution and the proportion in the always-valid subgroup
return estimates, conditional on population subgroups if requested

3.1 Simulate data

We first simulate some data for illustration.

The data has four variables

continuous confounder x
binary treatment a
binary mediator m
continuous outcome y
- y is NA when m = FALSE

library(tidyverse)
library(pstratreg)
data <- pstratreg_sim(n = 100)

#>             x     a     s          y
#> 1 -0.28163841 FALSE  TRUE  1.8702244
#> 2  0.09644999 FALSE  TRUE -2.0035928
#> 3 -0.80410616  TRUE  TRUE  0.8978557
#> 4 -2.50169791 FALSE FALSE         NA
#> 5 -0.39149331 FALSE FALSE         NA

3.2 The `pstratreg` function

The call below runs a principal stratification regression analysis with default options, returning estimates of the average treatment effect among the latent stratum who would have valid outcomes regardless of treatment.

result <- pstratreg(
  formula_y = formula(y ~ x*a),
  formula_s = formula(s ~ x*a),
  data = data,
  treatment_name = "a"
)

#> Effect on survival, where S = 1 indicates the outcome exists
#> # A tibble: 1 × 3
#>      s0    s1 effect_s
#>   <dbl> <dbl>    <dbl>
#> 1 0.768 0.846   0.0787
#> 
#> Effect on outcome among those who would have a valid outcome regardless of treatment
#> # A tibble: 1 × 2
#>   effect_y_lower effect_y_upper
#>            <dbl>          <dbl>
#> 1          0.109           1.73

3.3 Positive monotonicity

If you believe that the TRUE value of the treatment never causes the outcome to be undefined, then you might add the monotonicity_positive = TRUE option.

Sometimes, the monotonicity assumption disagrees with the empirical estimates in at least some cases. For example, we assume that treatment never prevents a valid outcome but empirically we estimate that the treatment increases the probability that the treatment increases the value of the mediator in some subgroups. When this happens, the package issues a warning.

Empirical monotonicity violations may be non-troubling; they can occur in estimates due to random chance from sampling variability. Because the user has assumed monotonicity, the package assumes that any violations arise purely from estimation errors. The predicted values of the mediator under each treatment condition in these cases are forced to be equal, at the midpoint of the two predicted values.

Generally, if the warning tells you that monotonicity is violated in only a small percent of cases, it may be warranted to proceed. If monotonicity is empirically violated in many cases, then you may need to rethink this assumption.

result <- pstratreg(
  formula_y = formula(y ~ x*a),
  formula_s = formula(s ~ x*a),
  data = data,
  treatment_name = "a",
  monotonicity_positive = TRUE
)
#> Warning in pstratreg(formula_y = formula(y ~ x * a), formula_s = formula(s ~ : Monotonicity violated in 42 % of cases
#> Forcing s1_trunc = s0_trunc at midpoint of estimates for those

#> Effect on survival, where S = 1 indicates the outcome exists
#> # A tibble: 1 × 3
#>      s0    s1 effect_s
#>   <dbl> <dbl>    <dbl>
#> 1 0.768 0.846   0.0787
#> 
#> Effect on outcome among those who would have a valid outcome regardless of treatment
#> # A tibble: 1 × 2
#>   effect_y_lower effect_y_upper
#>            <dbl>          <dbl>
#> 1          0.804           1.05

3.4 Negative monotonicity

Conversely, you can assume negative monotonicity with monotonicity_negative = TRUE. If you are not sure what negative monotonicity is, see the previous page on the big idea!

In this particular simulation, negative monotonicity does not hold and you see below that the warning appropriately alerts us that monotonicity is frequently empirically violated.

result <- pstratreg(
  formula_y = formula(y ~ x*a),
  formula_s = formula(s ~ x*a),
  data = data,
  treatment_name = "a",
  monotonicity_negative = TRUE
)
#> Warning in pstratreg(formula_y = formula(y ~ x * a), formula_s = formula(s ~ : Monotonicity violated in 58 % of cases
#> Forcing s1_trunc = s0_trunc at midpoint of estimates for those

#> Effect on survival, where S = 1 indicates the outcome exists
#> # A tibble: 1 × 3
#>      s0    s1 effect_s
#>   <dbl> <dbl>    <dbl>
#> 1 0.768 0.846   0.0787
#> 
#> Effect on outcome among those who would have a valid outcome regardless of treatment
#> # A tibble: 1 × 2
#>   effect_y_lower effect_y_upper
#>            <dbl>          <dbl>
#> 1          0.856          0.919

3.5 Aggregate in subgroups

Instead of sample average effect estimates, you might want the estimate within subgroups defined by a grouping variable in the data. The group_vars argument allows you to specify a vector of variable names from data within which to aggregate results.

First we create some groups for illustration

data_with_groups <- data %>%
  mutate(group1 = x < -.5,
         group2 = x > .5)

#>             x     a     s          y group1 group2
#> 1 -0.28163841 FALSE  TRUE  1.8702244  FALSE  FALSE
#> 2  0.09644999 FALSE  TRUE -2.0035928  FALSE  FALSE
#> 3 -0.80410616  TRUE  TRUE  0.8978557   TRUE  FALSE
#> 4 -2.50169791 FALSE FALSE         NA   TRUE  FALSE
#> 5 -0.39149331 FALSE FALSE         NA  FALSE  FALSE
#> 6  0.24832642  TRUE FALSE         NA  FALSE  FALSE

and then we apply the function to estimate within those groups.

result <- pstratreg(
  formula_y = formula(y ~ x*a),
  formula_s = formula(s ~ x*a),
  data = data_with_groups,
  treatment_name = "a",
  group_vars = c("group1","group2")
)

#> Effect on survival, where S = 1 indicates the outcome exists
#> # A tibble: 3 × 5
#>   group1 group2    s0    s1 effect_s
#>   <lgl>  <lgl>  <dbl> <dbl>    <dbl>
#> 1 FALSE  FALSE  0.806 0.852   0.0468
#> 2 FALSE  TRUE   0.941 0.913  -0.0273
#> 3 TRUE   FALSE  0.446 0.733   0.287 
#> 
#> Effect on outcome among those who would have a valid outcome regardless of treatment
#> # A tibble: 3 × 4
#> # Groups:   group1, group2 [3]
#>   group1 group2 effect_y_lower effect_y_upper
#>   <lgl>  <lgl>           <dbl>          <dbl>
#> 1 FALSE  FALSE          0.0792           1.51
#> 2 TRUE   FALSE         -2.56             3.44
#> 3 FALSE  TRUE           0.829            1.46

3.6 Sample weights

If you have case weights from sampling with unequal probabilities, they can be provided in a vector of length nrow(data) using the weights argument.

Here we generate some hypothetical weights

sim_weights <- runif(nrow(data))

and then call the function. Note that the glm() used internally to estimate logistic regression will create a warning for non-integer #successes in a binomial glm! which is to be expected when weights are used in this function.

result <- pstratreg(
  formula_y = formula(y ~ x*a),
  formula_s = formula(s ~ x*a),
  data = data,
  treatment_name = "a",
  weights = sim_weights
)
#> Warning in eval(family$initialize): non-integer #successes
#> in a binomial glm!

#> Effect on survival, where S = 1 indicates the outcome exists
#> # A tibble: 1 × 3
#>      s0    s1 effect_s
#>   <dbl> <dbl>    <dbl>
#> 1 0.774 0.801   0.0277
#> 
#> Effect on outcome among those who would have a valid outcome regardless of treatment
#> # A tibble: 1 × 2
#>   effect_y_lower effect_y_upper
#>            <dbl>          <dbl>
#> 1         0.0985           1.84

2 Why regression

4 Relaxing homoskedasticity