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.05982853 FALSE  TRUE -0.8101817
#> 2 -1.33520583 FALSE  TRUE -2.2969309
#> 3 -0.82352375  TRUE  TRUE -1.8947736
#> 4  0.58452730  TRUE FALSE         NA
#> 5  0.55000337  TRUE  TRUE  1.6149997

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.710 0.870    0.160
#> 
#> 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.107           1.66

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
)

#> Effect on survival, where S = 1 indicates the outcome exists
#> # A tibble: 1 × 3
#>      s0    s1 effect_s
#>   <dbl> <dbl>    <dbl>
#> 1 0.710 0.870    0.160
#> 
#> 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.625           1.14

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 100 % 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.710 0.870    0.160
#> 
#> 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.867          0.867

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.05982853 FALSE  TRUE -0.8101817  FALSE  FALSE
#> 2 -1.33520583 FALSE  TRUE -2.2969309   TRUE  FALSE
#> 3 -0.82352375  TRUE  TRUE -1.8947736   TRUE  FALSE
#> 4  0.58452730  TRUE FALSE         NA  FALSE   TRUE
#> 5  0.55000337  TRUE  TRUE  1.6149997  FALSE   TRUE
#> 6  0.44628076  TRUE  TRUE  2.8372122  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.749 0.887   0.138 
#> 2 FALSE  TRUE   0.880 0.934   0.0539
#> 3 TRUE   FALSE  0.476 0.780   0.304 
#> 
#> 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.170           1.53
#> 2 TRUE   FALSE          -1.36            2.45
#> 3 FALSE  TRUE            0.731           1.44

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.700 0.833    0.133
#> 
#> 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.0558           2.00

2 Why regression

4 Relaxing homoskedasticity