Is the emmeans R package performing causal inference G-computation? So I am trying to get an understanding of causal inference and how it differs from the usual contrasts. I regularly use the emmeans package in R, and I am wondering when the function emmeans() mentions it has averaged over the covariates is this essentially performing G-computation? At least for regular OLS or identity link GLM models it seems like that
One possible difference I see is that G-computation takes place on the response scale, so I wonder if you use the transform argument in emmeans when working with GLMs with a non-identity link, then would it be performing the G-computation? It seems like it does the delta method to convert the SEs to the response scale from the link scale.
And then if I modeled the treatment probability, inverted it, and used it as a covariate in a model, and then used emmeans -- would this be doubly robust estimation?
As my reference I am using the book here https://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/ pages 163-167.
Edit 09/07/22: Anyone wondering about the connection --it is written about here https://vincentarelbundock.github.io/marginaleffects/articles/gformula.html. The marginaleffects package can do this
 A: Re-reading your question, my understanding is that you are asking if emmeans() does G-computation as part of what it ordinarily does. And based on my very limited understanding of causal models and G-computation, I would say the answer is NO. That is simply because we don't treat covariates in any special way. For a numerical covariate, the default action is to compute its mean and use that as a reference value for all subsequent estimates, regardless of whether it is regarded as a mediator or not. We just treat it as a direct effect.
There may be some options in emmeans() that do allow the user to treat covariates in a different way. For example, we can fit a model y ~ treat + M where treat is a treatment and M is a mediator. Then suppose we subsequently do
emmeans(model, "treat", cov.reduce = M ~ treat)

This instructs emmeans to not use the average value of M, but rather to use lm() to fit the model M ~ treat (with the same dataset) and use its predictions for the value of M. In that way, the reference value of M is different for each treatment level. This is equivalent to creating a covariate C that is equal to the residuals of the M ~ treat model, fitting the model y ~ treat + C, and using emmeans() in the ordinary way by using C's mean (which is zero) as the reference value. Perhaps this is similar to what G-computation does -- I am not sure, but perhaps someone else can shed some light on this. But at least it does something special with covariates thought to be mediators, and that seems more akin to what is needed in causal inference.
Addendum
A comment to this answer suggests doing something like emmeans(model, "treat", cov.reduce = FALSE, weights = "prop") but that this is very inefficient as it creates a huge reference grid. I believe that the following may do the same thing:
emmeans(model, "treat", submodel = ~ treat)

The above puts a linear constraint on the estimates whereby all the effects other than those of treat are replaced by predictions of those effects from the given submodel. See vignette("xplanations", "emmeans") for the gory details. But in words, what happens is that we are trying to obtain the predictions we would have obtained from the submodel, while still accounting for the reduction in error variance achieved by including the covariate in the model. I think this in fact does relate to some causal-inference methods, but I lack the depth of knowledge in that are to be sure.
In the case of mixed models and generalized linear models, the submodel constraint will not be quite the same as would be obtained by fitting the submodel with the same method. To accomplish this (or at least get closer), one can use a new feature in version 1.6.0 of emmeans to bring in covariate predictions from an external model. Suppose model was fitted using something like model <- lmer(y ~ M + treat + (1|SUBJ), ...)
Mmod <- lmer(M ~ treat + (1|SUBJ, ...)
Mpred <- function(grid)
             list(M = predict(Mmod, newdata = grid, re.form = ~ 0))

emmeans(model, "treat", cov.reduce = extern ~ Mpred)

This is like the original cov.reduce = M ~ treat, except it uses Mmod instead of lm(M ~ treat) to do the predictions of M.
