Using models to compute the law of total probabilities

[lorenzoFabbri]

I am going through an example of the parametric G-formula here (in the PDF it starts on page 16). We simply want to use the G-formula for time-invariant exposures, to compute the marginal mean of the outcome.

# arrange into long data
C <- c(0, 0, 1, 1)
A <- c(0, 1, 0, 1)
Y <- c(94.3, 119.2, 130.6, 155.7)
N <- c(344052, 154568, 154560, 346820)
D <- NULL
for (i in 1:4) {
  d <- data.frame(cbind(rep(C[i], N[i]), rep(A[i], N[i]), rep(Y[i], N[i])))
  D <- rbind(D, d)
}
names(D) <- c("C", "A", "Y")
D <- tibble::as_tibble(D)

# fit models
mC <- glm(C ~ 1, data = D, family = binomial("logit"))
mA <- glm(A ~ C, data = D, family = binomial("logit"))
mY <- glm(Y ~ A + C, data = D, family = gaussian("identity"))

# obtain C predictions
pC <- predict(mC, type = "response")
# use predicted C to obtain predicted A
pA <- predict(mA, newdata = data.frame(C = pC), type = "response")
# use predicted A and C to obtain predicted Y
pY <- predict(mY, newdata = data.frame(A = pA, C = pC), type = "response")

# compute marginal mean of predicted Y
mean(pY)

In the PDF the author explains why we can use models and predictions to implement the law of total probability. The key is modeling each of exposure, confounder, and outcome, using predictions from the previous steps.

I do not understand then, how I can simply use marginaleffects::avg_predictions(mY)$estimate to obtain the same result (which is not really the same, since the 3rd decimal is different). I did not provide a model for $C$ and $A$, just for $Y$, and yet I believe this is how to implement the parametric G-formula with marginaleffects. I feel ashamed for the simplicity of this issue but I really do not understand.

@vincentarelbundock @ngreifer @ainaimi

[ngreifer]

This isn't really a marginaleffects question; marginaleffects::avg_predictions(mY)$estimate is just the same as mean(mY$fitted) which is always equal to mean(D$Y) when a linear outcome model is used (as you did here). So I think you need to ask why mean(pY) is not equal to mean(D$Y).

G-computation is done by running

avg_predictions(mY, variables = "A")

This computes the predictions of Y under each value of A and averages each over the empirical distribution of C, which is exactly what Ashley describes as g-computation on p15 in the formula starting with $E(Y^a) =$. I honestly have never seen a code or numerical implementation like the one you posted and I don't understand what it is supposed to compute or illustrate, so I can't speak to it.