Calculate risk ratio in weighted population

Question

I have a propensity score weighted population (using IPTW) and I want to compute risk ratios on my weighted population. For this, I am using a weighted Poisson regression.

Let's suppose that "married" is the outcome and I want to see the risk ratio of treated vs non treated considering married as outcome.

library(WeightIt)
library(survey)
W.out <- weightit(treat ~ age + educ + race + nodegree + re74 + re75,
                  data = lalonde, estimand = "ATE", method = "ps")

d.w <- svydesign(~1, weights = W.out$weights, data = lalonde)
fit <- svyglm(married ~ treat , design = d.w, family="poisson") 
exp(coefficients(fit))
exp(confint.default(fit))

For instance, in this case the weighted rate of "married" for each group is:

married_treated <- with(lalonde, weighted.mean(married[treat == 1], W.out$weights[treat == 1]))
married_nontreated <- with(lalonde, weighted.mean(married[treat == 0], W.out$weights[treat == 0]))

The result of the Poisson regression is basically equal to the ratio between married_treated / married_nontreated. This should be correct because RR is risk of outcome in exposed people / risk of outcome in unexposed.

The risk of outcome in exposed patients is equal to the frequency of married people in the total exposed group divided by the total n° of people in the exposed group (married/total n° exposed). Subsequently, the weighted rate of married in the treated group should be the weighted risk in the exposed while the weighted rate of married in the non treated group should be the weighted risk in the non exposed. Therefore, RR = weighted risk treated / weighted risk non treated. However, I would like a confirmation.

Is this the correct way of computing risk ratios in a population balanced through IPTW?

Answer 1

That's one way to do it. But if you include covariates in the outcome model or if you want to generalize the effect to a subset of your population (e.g., for the average treatment effect in the treated), you should use weighted g-computation instead. This is explained in the WeightIt vignette on estimating effects.

Weighted g-computation involves fitting an outcome model with the IPW weights applied, then generating predicted values of the outcome under treatment and under control for each unit. Then, you computed the IPW-weighted mean of the predicted outcomes under treatment and under control, which are the marginal risks. Finally, you can take the ratio of the marginal risks and that is your IPW- and covariate-adjusted estimate of the marginal risk ratio. A benefit of this method is that you can use whatever model you want to fit the outcome model because the risk ratio does not correspond to a coefficient in the outcome model but rather is a function of the covariates, the model parameters, and the weights. This is important when including covariates in the outcome model because a Poisson model is rarely if ever the right model for a binary outcome. It is useful as a convenience because the coefficient on treatment is equal to the log risk ratio, but that is only true when no covariates are present in the model. In contrast, a logistic regression model is more likely to fit the data well and generate valid predictions.

Below is how to do weighted g-computation manually. Getting standard errors is hard unless you bootstrap.

library(WeightIt)

data("lalonde", package = "cobalt")

W.out <- weightit(treat ~ age + educ + race + nodegree + re74 + re75,
                  data = lalonde, estimand = "ATE", method = "ps")

library(survey)
d.w <- svydesign(~1, weights = W.out$weights, data = lalonde)
fit <- svyglm(married ~ treat * (age + educ + race + nodegree + re74 + re75),
              design = d.w,  family = "quasibinomial")

# Predicted outcomes under treatment
pred_1 <- predict(fit, newdata = transform(lalonde, treat = 1),
                  type = "response")

# Predicted outcomes under control
pred_0 <- predict(fit, newdata = transform(lalonde, treat = 0),
                  type = "response")

# Marignal risks udner treatment and control
Epred_1 <- weighted.mean(pred_1, W.out$weights)
Epred_0 <- weighted.mean(pred_0, W.out$weights)

# Marginal risk ratio
(RR <- Epred_1 / Epred_0)
#> [1] 0.6167727

Below is how you do it using marginaleffects as explained in the WeightIt vignette.

# Using marginaleffects
library(marginaleffects)
avg_comparisons(fit, variables = "treat",
                comparison = "lnratioavg",
                transform = "exp",
                wts = "(weights)")
#> 
#>   Term              Contrast Estimate Pr(>|z|) 2.5 % 97.5 %
#>  treat ln(mean(1) / mean(0))    0.617    0.048 0.382  0.996
#> 
#> Columns: term, contrast, estimate, p.value, conf.low, conf.high, predicted, predicted_hi, predicted_lo

You should always use the marginaleffects strategy because it produces the right answer no matter what outcome model you use or whether you have covariates in the outcome model or not. When there are no covariates in the outcome model, it produces the same estimate as the exponentiated coefficient on treatment in a Poisson regression model.