Vignette updates

vignettes/WeightIt.Rmd

@@ -148,33 +148,28 @@ summary(Wmsm.out)
 Displayed are summaries of how the weights perform at each time point with respect to variability. Next, we'll examine how well they perform with respect to covariate balance.
 
 ```{r}
-bal.tab(Wmsm.out, stats = c("m", "ks"), which.time = .none)
+bal.tab(Wmsm.out, stats = c("m", "ks")
+        which.time = .none)
 ```
 
 By setting `which.time = .none` in `bal.tab()`, we can focus on the overall balance assessment, which displays the greatest imbalance for each covariate across time points. We can see that our estimated weights balance all covariates all time points with respect to means and KS statistics. Now we can estimate our treatment effects.
 
 First, we fit a marginal structural model for the outcome using `glm()` with the weights included:
 
 ```{r}
-# Attach weights to dataset
-msmdata$weights <- Wmsm.out$weights
-
 # Fit outcome model
-fit <- glm(Y_B ~ A_1 * A_2 * A_3 * (X1_0 + X2_0),
-           data = msmdata, weights = weights,
-           family = quasibinomial)
-
+fit <- glm_weightit(Y_B ~ A_1 * A_2 * A_3 * (X1_0 + X2_0),
+                    data = msmdata,
+                    weightit = Wmsm.out,
+                    family = binomial)
 ```
 
 Then, we compute the average expected potential outcomes under each treatment regime using `marginaleffects::avg_predictions()`:
 
 ```{r, eval = me_ok}
 library("marginaleffects")
 (p <- avg_predictions(fit,
-                      vcov = "HC3",
-                      newdata = datagridcf(A_1 = 0:1, A_2 = 0:1, A_3 = 0:1),
-                      by = c("A_1", "A_2", "A_3"),
-                      type = "response"))
+                      variables = c("A_1", "A_2", "A_3")))
 ```
 
 We can compare the expected potential outcomes under each regime using `marginaleffects::hypotheses()`. To get all pairwise comparisons, supply the `avg_predictions()` output to `hypotheses(., "pairwise")`. To compare individual regimes, we can use `hypotheses()`, identifying the rows of the `avg_predictions()` output. For example, to compare the regimes with no treatment for all three time points vs. the regime with treatment for all three time points, we would run

vignettes/estimating-effects.Rmd

@@ -131,7 +131,9 @@ Our primary focus will be on marginal effects, which are appropriate for all eff
 
 To estimate marginal effects, we use a method known as g-computation [@snowdenImplementationGComputationSimulated2011] or regression estimation [@schaferAverageCausalEffects2008]. This involves first specifying a model for the outcome as a function of the treatment and covariates. Then, for each unit, we compute their predicted values of the outcome setting their treatment status to treated, and then again for control, leaving us with two predicted outcome values for each unit, which are estimates of the potential outcomes under each treatment level. We compute the mean of each of the estimated potential outcomes across the entire sample, which leaves us with two average estimated potential outcomes. Finally, the contrast of these average estimated potential outcomes (e.g., their difference or ratio, depending on the effect measure desired) is the estimate of the treatment effect.
 
-When doing g-computation after weighting, a few additional considerations are required. First, the outcome model should be fit incorporating the estimated weights (e.g., using weighted least squares or weighted maximum likelihood estimation) [@vansteelandtInvitedCommentaryGComputation2011]. Second, when we take the average of the estimated potential outcomes under each treatment level, if an estimand other than the ATT, ATC, or ATE is used or sampling weights are included, this must be a weighted average that incorporates the weights\^[Note: Previously, this guide recommended always including the weights. We believe best practice is to omit them, though this does not mean previous analyses using weighted averages are invalid.]. Third, if we want to target the ATT or ATC, we only estimate potential outcomes for the treated or control group, respectively (though we still generate predicted values under both treatment and control) [@wangGcomputationAverageTreatment2017].
+When doing g-computation after weighting, a few additional considerations are required. First, the outcome model should be fit incorporating the estimated weights (e.g., using weighted least squares or weighted maximum likelihood estimation) [@vansteelandtInvitedCommentaryGComputation2011]. Second, when we take the average of the estimated potential outcomes under each treatment level, if an estimand other than the ATT, ATC, or ATE is used or sampling weights are included, this must be a weighted average that incorporates the weights[^1]. Third, if we want to target the ATT or ATC, we only estimate potential outcomes for the treated or control group, respectively (though we still generate predicted values under both treatment and control) [@wangGcomputationAverageTreatment2017].
+
+[^1]: Note: Previously, this guide recommended always including the weights. We believe best practice is to omit them, though this does not mean previous analyses using weighted averages are invalid.
 
 G-computation as a framework for estimating effects after weighting has a number of advantages over other approaches. It works the same regardless of the form of the outcome model or type of outcome (e.g., whether a linear model is used for a continuous outcome or a logistic model is used for a binary outcome); the only difference might be how the average expected potential outcomes are contrasted in the final step. In simple cases, the estimated effect is numerically identical to effects estimated using other methods; for example, if no covariates are included in the outcome model, the g-computation estimate is equal to the difference in means from a t-test or coefficient of the treatment in a linear model for the outcome. There are analytic and bootstrap approximations to the SEs of the g-computation estimate. The analytic approximation is computed using the delta method, a technique for computing the SE of a quantity derived from the regression model parameters, such as the g-computation estimate.
 
@@ -141,7 +143,7 @@ There are other methods to incorporate the outcome model into estimation of the
 
 ### Modeling the Outcome
 
-The goal of the outcome model is to generate good predictions for use in the g-computation procedure described above. The type and form of the outcome model should depend on the outcome type. For continuous outcomes, one can use a linear model regressing the outcome on the treatment; for binary outcomes, one can use a generalized linear model with, e.g., a logistic link; for time-to-event outcomes, one can use a Cox proportional hazards model. Note that the purpose of including the outcome model is not to arrive at a doubly robust estimator (i.e., one that is consistent if either the outcome or propensity score model is correct); rather, it is simply to increase the precision of the weighted estimate essentially for free. To take advantage of this feature, it is important to use a canonical link (i.e., the default link for a given family), as recommended by @gabrielInverseProbabilityTreatment2023.
+The goal of the outcome model is to generate good predictions for use in the g-computation procedure described above. The type and form of the outcome model should depend on the outcome type. For continuous outcomes, one can use a linear model regressing the outcome on the treatment; for binary outcomes, one can use a generalized linear model with, e.g., a logistic link; for time-to-event outcomes, one can use a Cox proportional hazards model. Note that the purpose of including the outcome model is not to arrive at a doubly robust estimator (i.e., one that is consistent if either the outcome or propensity score model is correct); rather, it is simply to increase the precision of the weighted estimate essentially for free. To take advantage of this feature, it is important to use a canonical link (i.e., the default link for a given family), as recommended by @gabrielInverseProbabilityTreatment2024.
 
 An additional decision to make is whether (and how) to include covariates in the outcome model. One may ask, why use weighting at all if you are going to model the outcome with covariates anyway? Weighting reduces the dependence of the effect estimate on correct specification of the outcome model; this is the central thesis of @hoMatchingNonparametricPreprocessing2007 (though applied to matching in their case). Including covariates in the outcome model after weighting has several functions: it can increase precision in the effect estimate, reduce the bias due to residual imbalance, and make the effect estimate "doubly robust", which means it is consistent if either the weighting reduces sufficient imbalance in the covariates or if the outcome model is correct. For these reasons, we recommend covariate adjustment after weighting when possible. There is some evidence that covariate adjustment is most helpful for covariates with standardized mean differences greater than .1 [@nguyenDoubleadjustmentPropensityScore2017], so these covariates and covariates thought to be highly predictive of the outcome should be prioritized in treatment effect models if not all can be included due to sample size constraints.
 
@@ -157,7 +159,7 @@ Uncertainty estimation (i.e., of SEs, confidence intervals, and p-values) may co
 
 M-estimation involves specifying a stack of estimating equations, one for each parameter to be estimated, that are each a function of the parameters. The final parameter estimates are those that yield a vector of 0s for the solutions of the estimating equations. These estimating equations include the parameters of the weighting model (e.g., the coefficients in the logistic regression model for the propensity score) and the parameters of the outcome model (i.e., the coefficients on treatment and the covariates). It is possible to compute the joint covariance matrix of all the estimated parameters using functions of the estimating equations and the estimated parameters. We refer curious readers to @stefanskiCalculusMEstimation2002 for an introduction.
 
-Most theoretical developments on the estimation of asymptotically correct standard errors for propensity score-weighted treatment effect estimates rely on M-estimation theory [@luncefordStratificationWeightingPropensity2004; @reifeisVarianceTreatmentEffect2020; @gabrielInverseProbabilityTreatment2023]. This method can only be used when the model used to estimate the weights involves estimating equations. Currently, only weights estimated using a generalized linear model propensity score, entropy balancing, inverse probability tilting, or the just-identified covariate balancing propensity score can be accommodated with this method. `glm_weightit()` can be used to fit (generalized) linear models for the outcome that account for estimation of the weights. When estimating equations are not available for a given method, the weights are treated as fixed and known, and the M-estimation standard errors are equal to the robust standard errors described below.
+Most theoretical developments on the estimation of asymptotically correct standard errors for propensity score-weighted treatment effect estimates rely on M-estimation theory [@luncefordStratificationWeightingPropensity2004; @reifeisVarianceTreatmentEffect2020; @gabrielInverseProbabilityTreatment2024]. This method can only be used when the model used to estimate the weights involves estimating equations. Currently, only weights estimated using a generalized linear model propensity score, entropy balancing, inverse probability tilting, or the just-identified covariate balancing propensity score can be accommodated with this method. `glm_weightit()` can be used to fit (generalized) linear models for the outcome that account for estimation of the weights. When estimating equations are not available for a given method, the weights are treated as fixed and known, and the M-estimation standard errors are equal to the robust standard errors described below.
 
 #### Robust Standard Errors
 
@@ -246,9 +248,9 @@ If, in addition to the effect estimate, we want the average estimated potential
 avg_predictions(fit, variables = "A")
 ```
 
-We can see that the difference in potential outcome means is equal to the average treatment effect computed previously[^1]. The arguments to `avg_predictions()` are the same as those to `avg_comparisons()`.
+We can see that the difference in potential outcome means is equal to the average treatment effect computed previously[^2]. The arguments to `avg_predictions()` are the same as those to `avg_comparisons()`.
 
-[^1]: To verify that they are equal, supply the output of `avg_predictions()` to `hypotheses()`, e.g., `avg_predictions(…) |> hypotheses("revpairwise")`; this explicitly compares the average potential outcomes and should yield identical estimates to the `avg_comparisons()` call.
+[^2]: To verify that they are equal, supply the output of `avg_predictions()` to `hypotheses()`, e.g., `avg_predictions(…) |> hypotheses("revpairwise")`; this explicitly compares the average potential outcomes and should yield identical estimates to the `avg_comparisons()` call.
 
 ### Adjustments to the Standard Case
 
@@ -276,9 +278,9 @@ wts = W$weights
 
 Estimating effects on binary outcomes is essentially the same as for continuous outcomes. The main difference is that there are several measures of the effect one can consider, which include the odds ratio (OR), risk ratio/relative risk (RR), and risk difference (RD), and the syntax to `avg_comparisons()` depends on which one is desired. The outcome model should be one appropriate for binary outcomes (e.g., logistic regression) but is unrelated to the desired effect measure because we can compute any of the above effect measures using `avg_comparisons()` after the logistic regression.
 
-To fit a logistic regression model, change `lm_weightit()` to `glm_weightit()` and set `family = binomial`. To compute the marginal RD, we can use exactly the same syntax as in the Standard Case; nothing needs to change[^2].
+To fit a logistic regression model, change `lm_weightit()` to `glm_weightit()` and set `family = binomial`. To compute the marginal RD, we can use exactly the same syntax as in the Standard Case; nothing needs to change[^3].
 
-[^2]: Note that for low or high average expected risks computed with `predictions()`, the confidence intervals may go below 0 or above 1; this is because an approximation is used. To avoid this problem, bootstrapping or simulation-based inference can be used instead.
+[^3]: Note that for low or high average expected risks computed with `predictions()`, the confidence intervals may go below 0 or above 1; this is because an approximation is used. To avoid this problem, bootstrapping or simulation-based inference can be used instead.
 
 To compute the marginal log RR, we need to add `comparison = "lnratioavg"` to `avg_comparisons()`; this computes the marginal log RR. To get the marginal RR, we need to add `transform = "exp"` to `avg_comparisons()`, which exponentiates the marginal log RR and its confidence interval. The code below computes the effects and displays the statistics of interest:
 
@@ -457,7 +459,9 @@ ggplot(p, aes(x = Ac)) +
 
 We can see that the ADRF has an elbow-shape, with a zone of no evidence of treatment effect followed by a zone of increasing values of the outcome as the treatment increases.
 
-Another way to characterize the effect of continuous treatments is to examine the average marginal effect function (AMEF), which is a function that relates the value of treatment to the derivative of the ADRF. When this derivative is different from zero, there is evidence of a treatment effect at the corresponding value of treatment. Below, we use `avg_slopes()` to compute the pointwise derivatives of the ADRF across levels of `Ac` and then plot it\^[You can also use `plot_slopes()`.].
+Another way to characterize the effect of continuous treatments is to examine the average marginal effect function (AMEF), which is a function that relates the value of treatment to the derivative of the ADRF. When this derivative is different from zero, there is evidence of a treatment effect at the corresponding value of treatment. Below, we use `avg_slopes()` to compute the pointwise derivatives of the ADRF across levels of `Ac` and then plot it[^4].
+
+[^4]: You can also use `plot_slopes()`
 
 ```{r, fig.height=3.5, fig.width=7}
 # Estimate the pointwise derivatives at representative

vignettes/references.bib

@@ -478,3 +478,17 @@ @article{xuApplicationsFractionalRandomWeightBootstrap2020a
 publisher: Taylor & Francis},
 	langid = {en}
 }
+
+@article{
+  gabrielInverseProbabilityTreatment2024,
+  title={Inverse probability of treatment weighting with generalized linear outcome models for doubly robust estimation},
+  volume={43},
+  ISSN={1097-0258},
+  url={https://onlinelibrary.wiley.com/doi/abs/10.1002/sim.9969},
+  DOI={10.1002/sim.9969},
+  note={_eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/sim.9969},
+  number={3},
+  journal={Statistics in Medicine},
+  author={Gabriel, Erin E. and Sachs, Michael C. and Martinussen, Torben and Waernbaum, Ingeborg and Goetghebeur, Els and Vansteelandt, Stijn and Sjölander, Arvid},
+  year={2024},
+pages={534–547} }