3 Continuous Outcome: Average Treatment Effect (ATE)
3.1 Computing Weights
The OVtool can either take a vector of weights estimated using any method or a propensity score (ps) object produced by TWANG (Ridgeway et al., 2021). The chunk of code below demonstrates how to specify your propensity score model and generate propensity score weights using the ps function from TWANG. If you already have a column of weights in your dataset, you can skip to step 2.
## Create Formula
my_formula = as.formula(treat ~ eps7p_0 + sfs8p_0 + sati_0 + ada_0 + recov_0 +
tss_0 + dss9_0)
## Get weights
# library(twang)
ps.twang <- twang::ps(my_formula, data = sud, estimand = 'ATE', booster = "gbm",
stop.method = "ks.max", verbose=F, ks.exact = T)The following line of code produces a balance table that demonstrates TWANG does a reasonable job of balancing. There are additional diagnostics we could check to ensure we have sufficient balance, but we move on without diving in further because the purpose of this tutorial is to showcase OVtool. See Ridgeway et al. for further information on balance diagnostics.
# Check Balance
twang::bal.table(ps.twang)#> $unw
#> tx.mn tx.sd ct.mn ct.sd std.eff.sz stat p ks ks.pval
#> eps7p_0 0.256 0.196 0.219 0.192 0.187 5.937 0.000 0.103 0.000
#> sfs8p_0 11.253 13.134 10.571 12.162 0.054 1.703 0.089 0.045 0.032
#> sati_0 8.233 22.128 2.145 10.658 0.345 11.088 0.000 0.121 0.000
#> ada_0 48.748 33.400 54.236 32.454 -0.166 -5.271 0.000 0.081 0.000
#> recov_0 0.246 0.431 0.240 0.427 0.015 0.479 0.632 0.006 1.000
#> tss_0 2.277 3.525 1.924 3.115 0.106 3.365 0.001 0.043 0.050
#> dss9_0 2.750 2.604 2.638 2.492 0.044 1.390 0.165 0.023 0.666
#>
#> $ks.max.ATE
#> tx.mn tx.sd ct.mn ct.sd std.eff.sz stat p ks ks.pval
#> eps7p_0 0.239 0.194 0.233 0.192 0.028 0.880 0.379 0.021 0.765
#> sfs8p_0 10.885 12.655 10.727 12.372 0.012 0.392 0.695 0.012 0.998
#> sati_0 5.345 17.937 4.198 15.691 0.065 1.797 0.072 0.023 0.662
#> ada_0 51.455 32.894 52.377 32.898 -0.028 -0.858 0.391 0.020 0.823
#> recov_0 0.248 0.432 0.238 0.426 0.023 0.716 0.474 0.010 1.000
#> tss_0 2.099 3.358 2.065 3.279 0.010 0.313 0.754 0.012 0.999
#> dss9_0 2.675 2.556 2.706 2.535 -0.012 -0.381 0.703 0.012 0.999The last step needed prior to moving onto the specification of your outcome model is to attach the weights as a column in your dataset. We extract the weights produced from ps in the line below.
# Get weights (not needed if user inserts a ps object in OVtool::outcome_model)
sud$w_twang = ps.twang$w$ks.max.ATE3.2 Running Outcome Model
The next step is to estimate the observed treatment effect. The outcome_model function in OVtool automatically uses the svydesign and svyglm functions from the survey package in R (survey) and formats and returns a list of components that are needed throughout step 3.
There are two options the researcher can take to input the relevant information to get their outcome results using outcome_model().
- Input a
ps.objectfromTWANGand astop.method(e.g."ks.max") or - Input a data frame containing the
dataused, the column name representing the weights,weights, and the column name representing the treatment indicator,treatment.
The analyst must also provide:
- a column name representing the outcome,
outcome - a vector of covariates to be included in the final outcome model,
model_covariates - the estimand of interest,
estimand
# Run Models -- first standardize outcome
sud$eps7p_3_std = sud$eps7p_3/sd(sud$eps7p_3)
# Run outcome model (function in OVtool that calls svyglm from survey)
results = outcome_model(ps_object = NULL,
stop.method = NULL,
data = sud,
weights = "w_twang",
treatment = "treat",
outcome = "eps7p_3_std",
model_covariates = c("eps7p_0", "sfs8p_0",
"sati_0", "ada_0",
"recov_0", "tss_0",
"dss9_0"),
estimand = "ATE")
summary(results$mod_results)#>
#> Call:
#> svyglm(formula = formula, design = design_u)
#>
#> Survey design:
#> survey::svydesign(ids = ~1, weights = ~w_orig, data = data)
#>
#> Coefficients:
#> Estimate Std. Error t value Pr(>|t|)
#> (Intercept) 0.4235011 0.0657531 6.441 1.33e-10 ***
#> treat 0.0712001 0.0278454 2.557 0.0106 *
#> eps7p_0 1.7397954 0.1166266 14.918 < 2e-16 ***
#> sfs8p_0 0.0015028 0.0020843 0.721 0.4709
#> sati_0 0.0028393 0.0012348 2.299 0.0215 *
#> ada_0 -0.0005808 0.0007725 -0.752 0.4522
#> recov_0 -0.0704673 0.0317912 -2.217 0.0267 *
#> tss_0 0.0329637 0.0069454 4.746 2.15e-06 ***
#> dss9_0 0.0540816 0.0080153 6.747 1.72e-11 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> (Dispersion parameter for gaussian family taken to be 0.7074412)
#>
#> Number of Fisher Scoring iterations: 2The outcome model results show an adjusted treatment effect estimate that accounts for confounding from observed covariates between youth in the two treatment programs (A = 1 and B = 0). From the results, we observe that the estimated treatment effect is 0.07 (p = 0.011), whereby youth receiving treatment A have slightly higher emotional problems at the 3-month follow-up than youth in treatment program B.
3.3 Sensitivity Analysis
At this stage, researchers may ask themselves if the observed effect in Step 2 is real and how sensitive to unobserved confounders it is. Our tool is used to help answer this logical next step question.
3.3.1 Run the sensitivity analysis
The following snippet of code presents the main function in OVtool: ov_sim(). This function requires results from outcome_model() plus additional parameters including:
plot_covariates: a vector of column names potentially representing the covariates used to produce the analysts’ propensity score weights (these may or may not be the same as the list of covariates used for the outcome model)es_grid: a vector on an effect size scale representing the range of association between the unobserved confounder (omitted variable) and the treatment indicator over which to run the simulations.rho_grid: a vector of absolute correlations to simulate over. These correlations represent the range of absolute correlations between the omitted variable and the outcomen_reps: the number of repetitions represents the number of times an unobserved confounder is simulated at each effect size and rho combination. The package defaults to 50. We find that fifty repetitions is sufficient in most cases, but the analyst may need to increase the number of repetitions if the plotted contours are not sufficiently smooth.
Optional arguments include:
progress: whether or not the function progress should print to screen. The default value is TRUE. If the user does not want the output to print to the console, they should set this parameter to FALSE.add: The default value is FALSE. This is only set to TRUE if the user is running additional repetitions after the first call toov_sim. It should be used if the original value ofn_repsdoes not yield smooth contours.sim_archive: The default value is NULL and is only used ifadd_repsis called. This parameter should never be defined by the user.
The ov_sim function does not require the user to specify vectors for es_grid and rho_grid. If the user specifies NULL for both, the tool will automatically suggest appropriate values. We define an effect size value in es_grid to show the strength of the relationship between the simulated unobserved covariate (U) and the treatment group indicator; it is defined as the standardized mean difference in U for the treatment A and treatment B groups. Typical rules of thumb for effect sizes (Cohen’s D) follow such that effect sizes greater than 0.2 would be considered small, 0.5 would be moderate and 0.8 would be large (Cohen, J., 1988). We define a particular value of rho specified in rho_grid as the absolute correlation the unobserved covariate (U) has with the outcome of interest, with larger values indicating stronger relationships between U and the outcome.
Please see Burgette et al. (2021) for additional details on the theory behind OVtool.
# Run OVtool (with weights (not a ps object))
ovtool_results_twang = ov_sim(model_results=results,
plot_covariates=c("eps7p_0", "sfs8p_0",
"sati_0", "ada_0",
"recov_0", "tss_0",
"dss9_0"),
es_grid = NULL,
rho_grid = seq(0, 0.40, by = 0.05),
n_reps = 50,
progress = TRUE,
add = FALSE,
sim_archive = NULL)The grid values are not required; if es_grid and/or rho_grid are set to NULL, the tool will calculate reasonable values to simulate over. For this tutorial we set es_grid to NULL. The tool then iterated over all observed covariate associations with the treatment indicator on an effect size scale to derive a reasonable range.
The function, ov_sim, produced a warning saying “You specified a rho grid whose maximum value is less than the maximum absolute correlation at least one observed covariate has with the outcome. The rho grid was automatically expanded to include all plot_covariates specified in the relevant graphics. If you want the rho grid range to remain from 0 to 0.4 then you must exclude the following variables from the plot_covariates argument: eps7p_0, tss_0, dss9_0.” The grid was expanded to ensure all plot_covariates could be seen on the contour plot. If the user does not want the grid expanded, they can leave out the observed covariates used in the propensity score model that have an absolute correlation with the outcome that is greater than the maximum rho value the user specifies (in this example: eps_7p_0, dss9_0, and tss_0).
There are a few methodological assumptions that are important for an analyst to understand.
First, when generating the omitted variable (U), the empirical cumulative distribution function (CDF) for the residuals within each treatment is used. The residuals are generated from the observed outcome model. When multiple ties are present, the tool handles this by randomly ranking multiple ties. For example, imagine a scenario where you have ten observations (shown in rank order) and the first two observations are 0s. The tool will randomly assign the rank of the two zeroes and assign the first observation with 0.1 and the second observation with 0.2 as opposed to assigning the two zeroes a value of 0.2. This process is repeated for each effect size and rho combination n_reps times.
Another key assumption this method draws upon is that the omitted variable is independent from all observed covariates included in the propensity score model.
3.3.2 Visualize the sensitivity analysis
To visualize our results, the plot.ov function can be called to produce three graphics. The main arguments to plot.ov are as follows:
x: the object returned from the call toov_simcol: if specified as “color” a color graphic will be produced and if specified as “bw” a black and white graphic will be produced. If and only ifprint_graphicis set to “1”, the “col” option will produce a heat mapp_contours: a vector of p-value contours plotted. The default is 0.01, 0.05, and 0.1. We only recommend changing from the default if the p-value for the raw effect is close to one of these values. This parameter is only used whenprint_graphicequals “2” or “3”.print_graphic: takes values “1”, “2”, or “3”, representing Figures 1, 2, and 3 described below.The first graphic (Figure 1) plots the treatment effect contours without covariate labels. This graphic can be produced by specifying
print_graphicto “1”. If the user specifies the parametercolas"color", the contours will overlay a colored heat map.The second graphic (Figure 2) plots the treatment effect contours in black and the p-value contours in red. The p-value contours can be specified with
p_contours. The default values are 0.01, 0.05, and 0.1.The third graphic (Figure 3) plots the treatment effect contours with the p-value contour overlaid and covariate labels (i.e. the column names submitted to
plot_covariates) plotted by their raw rho and effect size.
The following code snippet will produce Figure 1.
plot.ov(ovtool_results_twang, print_graphic = "1", col = "bw")
The grid, as shown by the x-axis and y-axis in Figure 1 presents the effect size and rho, respectively.
The y-axis in Figure 1 represents the unobserved confounder’s absolute correlation with the outcome and the x-axis is the association between the unobserved confounder and the treatment indicator on an effect size scale. The black lines represent treatment effect contours that run along the grid. The PS weighted treatment effect of Treatment A versus Treatment B equals 0.071 and is significant with a p-value equal to 0.011. However, looking at this graphic alone will not give us an idea of how sensitive statistical significance is.
plot.ov(ovtool_results_twang, print_graphic = "2", col = "color")
#> [1] "NOTE: Covariates with absolute correlation with
outcome greater than 0.4: eps7p_0 (Actual: 0.509), tss_0
(Actual: 0.423), dss9_0 (Actual: 0.420)"Figure 2 is a different variation of Figure 1, but now adds p-value contours. This graphic will allow the user to see what treatment effect will switch the significance level at critical p-values (e.g., 0.05). This graphic will now give the user an idea of how sensitive both the treatment effect estimate and the statistical significance are. The analyst can specify what red p-value thresholds they want plotted with the p_contours argument.
plot.ov(ovtool_results_twang, print_graphic = "3", col = "color")
#> [1] "NOTE: Covariates with absolute correlation with
outcome greater than 0.4: eps7p_0 (Actual: 0.509), tss_0
(Actual: 0.423), dss9_0 (Actual: 0.420)"Figure 3 adds a final dimension (observed covariate labels) to Figure 2. Similar to Figures 1 and 2, plotted at the top of the figure margin is the ps weighted treatment effect (0.071) and associated p-value of 0.011. The solid black contours represent the treatment effect contour lines and the red lines (sometimes dashed) represent the p-value thresholds. The key on the right side of the graphic shows where various p-value cutoff lines are, including p = 0.05. The blue points on the plot represent the observed covariate correlations with the outcome (y-axis) and effect size associations with the treatment indicator (x-axis). For instance, ada_0 and the outcome have approximately a 0.18 absolute correlation with the emotional problem scale at three months and an absolute association of approximately 0.17 effect size difference between the two treatment groups (magnitude of its relationship with the treatment indicator). In this case, if an unobserved confounder had a relationship with the outcome (eps7p_3) and the treatment indicator (treat) that was equivalent to ada_0, tss_0, eps7p_0, sati_0, or dss9_0, then the researcher would conclude that their results would likely be sensitive to an unobserved confounder at a p-value threshold of 0.05. All of these observed relationships are sensitive at the 0.10 p-value threshold other than dss9_0. If the blue points all existed in contours greater than the 0.05 p-value contour, then unobserved confounders with similar associations would retain the significant effect and allow the user to conclude that the results are reasonably robust. Overall, the graphic above showcases an example where findings could be highly sensitive to the effects of an unobserved confounder.
Note: When the outcome model shows a significant effect for all observed covariates, regardless treatment effect direction, we force the sign of the magnitude to go with the direction of the significant effect. The blue points are meant to give the analyst an idea of possible associations (using observed covariates as an indicator) that would change the interpretation of their results.
3.3.3 Produce additional simulations
These results were produced with 50 simulations (n_reps) of the unobserved confounder. If the contours are not smooth or the user wants to add simulations, they can call add_reps and specify the number (more_reps) of additional simulations. Once this is completed, a user can recreate the plots. Example code to add repetitions:
# If you want to add repetitions, run the following line.
# ovtool_results_twang = add_reps(OVtool_results = ovtool_results_twang,
# model_results = results,
# more_reps = 30)
#
# Recreate Graphic
# plot.ov(ovtool_results_twang, print_graphic = "1", col = "bw")3.3.4 Recommendations for reporting the sensitivity analysis
Finally, we can interpret this graphic by running the summary command on the ov object:
summary.ov(object = ovtool_results_twang, model_results = results)#> [1] "Recommendation for reporting the sensitivity
analyses"
#> [1] "The sign of the estimated effect is expected to
remain consistent when simulated unobserved confounders
have the same strength of associations with the treatment
indicator and outcome that are seen in 6 of the 7 observed
confounders. In the most extreme observed case in which the
sign changes, the estimated treatment effect shifts from
0.071 to -0.025. The sign of the estimate would not be
expected to be preserved for unobserved confounders that
have the same strength of association with the treatment
indicator and outcome as eps7p_0."
#> [1] "Statistical significance at the 0.05 level is
expected to be robust to unobserved confounders with
strengths of associations with the treatment indicator and
outcome that are seen in 2 of the 7 observed confounders.
In the most extreme observed case, the p-value would be
expected to increase from 0.011 to 0.372. Significance at
the 0.05 level would not be expected to be preserved for
unobserved confounders that have the same strength of
association with the treatment indicator and outcome as
eps7p_0, sati_0, ada_0, tss_0, dss9_0."The OVtool gives a recommendation on how to report findings regarding the direction of the treatment effect and statistical significance. An analyst could take the results produced by summary.ov() and adapt the desciption to be included a manuscript. The tool produces these results by iterating through each observed confounder the user specified in the plot_covariates argument to ov_sim and runs the simulation for those exact effect size and rho combinations to produce the treatment effect and p-value. It then determines which observed relationships would produce a treatment effect that changes direction and p-value that crosses over the 0.05 threshold when a significant raw effect is observed.
In the next section, we will show how our method works with the average treatment effect on the treated (ATT) using a continuous outcome.
