The g-formula (g-computation algorithm) is a method to estimate the causal effect of a time varying treatment or covariates, developed by Robins (1986). Note this is distinct from g-estimation of structural nested models (Robins and Rotnitzky, 1992), a different g-methods estimator that uses instrumental-variable-like estimating equations and does not require modeling the joint distribution of all time-varying confounders. This chapter — and the gfoRmula/gmethods packages used below — implements the g-formula, not g-estimation.
We can imagine such a scenario. Suppose \(A0\) and \(A1\) are time 0 and time 1 treatment, can be binary or continuous. \(L1\) is a time 1 covariate. It can be some intermediate outcome. Say HbA1c. HbA1c can be a result of treatment \(A0\). Then it can be a cause for the next period treatment \(A1\). \(Y\) is the outcome. \(U\) is the unmeasured confounder. The causal diagram is shown below. How do we estimate the effect of treatment?
library(dagitty)
library(ggdag)
library(ggplot2)
library(gfoRmula)
library(gmethods)
library(tidyverse)
dag <- dagify(
Y ~ A1 + A0 + L1 ,
A1 ~ A0 + L1,
L1 ~ A0 ,
coords = time_ordered_coords()
)
dag %>%
ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
geom_dag_point() +
geom_dag_edges_arc() +
geom_dag_text() +
theme_dag()
There are four potential outcomes, with \(A_0\) and \(A_1\) can be 0 or 1. The causal estimand is \[ \tau_{a_0,a_1,a_0',a_1'} = E[Y^{a_0,a_1} - Y^{a_0',a_1'}] \tag{30.1}\]
To identify the effects, we need assumptions.
In general, treatment at time \(t\) is randomized with probabilities depending on the observed past (covariates and intermediate outcomes): at any time \(t\), for every regime \(\bar a\), \[ Y^{\bar a} \perp A_t \mid H_t \tag{30.4}\] where \(H_t\) is the observed history of treatment and covariates up to time \(t\), \[H_t=(A_0, A_1, ..., A_{t-1}; L_1,..., L_{t}). \tag{30.5}\] Read this in terms of the observed treatment \(A_t\) and the observed history \(H_t\): conditional on what has actually happened up to time \(t\), the next treatment is independent of the future potential outcome. The superscripted objects in the two-period statement above (e.g. \(A_1^{a_0}\), \(L_1^{a_0}\)) denote the values that would arise under setting \(A_0=a_0\); by consistency these equal the observed \(A_1\), \(L_1\) on the subset of units with \(A_0=a_0\), which is what lets the conditional-independence statement be checked on observed data.
\[ 0 < P(A_0=1) < 1 \tag{30.6}\] \[ 0 < P(A_1=1|A_0=1, L_1) < 1 \tag{30.7}\]
In general, \[ 0 < P(A_t=1|H_t) < 1 \tag{30.8}\]
\[ Y^{a_0,a_1} = Y | (A_0=a_0,A_1=a_1) \tag{30.9}\]
\[ \begin{aligned} E[Y^{a_0,a_1}] &= E[Y^{a_0,a_1} | A_0=a_0] \\ &= E[Y^{A_0,a_1} | A_0=a_0] \\ &= \sum_{l=0,1}E[Y^{A_0,a_1} | A_0=a_0, L=l] P(L=l | A_0=a_0) \\ &= \sum_{l=0,1}E[Y^{A_0,a_1} | A_0=a_0, L=l, A_1=a_1] P(L=l | A_0=a_0) \\ &= \sum_{l=0,1}E[Y^{A_0,A_1} | A_0=a_0, L=l, A_1=a_1] P(L=l | A_0=a_0) \end{aligned} \tag{30.10}\]
This can be expanded to the general case. Basically start from the first period, then the second, based on the first, and so on.
There are two types of components in the g-estimation. The first is the outcome model, then there are models for time-varying covariates, or intermediate outcomes (confounders).
In this simple case, if \(Y\) is continuous, we can use a linear model. If \(Y\) is binary, we can use a logistic model.
\[ E[Y | A_0, A_1, L_1] = \beta_0 + \beta_1 A_0 + \beta_2 L_1 + \beta_3 A_1 \tag{30.11}\]
Suppose \(L_1\) is binary,
\[P(L_1=1 | A_0) = logit^{-1}(\gamma_0 + \gamma_1 A_0) \tag{30.12}\] \[ L_1 | A_0 \sim Bernoulli\big(logit^{-1}(\gamma_0 + \gamma_1 A_0)\big) \tag{30.13}\]
We may need to simulate the distribution of the time-varying covariates for the expected value of the potential outcomes.
I am following Christopher Boyer’s example in his teaching page: https://christopherbboyer.com/about.
The first example is from his lab3 teaching page.
First we set up a simulated data set.
# setting up the data -----------------------------------------------------
# The code below creates the frequency table, and then turns it into long format
# data. Read through the code carefully and run it.
# re-create the frequency table for homework 3
hw3_freq <- tribble(
~A_0, ~L_1, ~A_1, ~N, ~Y_1,
0, 0, 0, 6000, 60,
0, 0, 1, 2000, 60,
0, 1, 0, 2000, 210,
0, 1, 1, 6000, 210,
1, 0, 0, 3000, 240,
1, 0, 1, 1000, 240,
1, 1, 0, 3000, 120,
1, 1, 1, 9000, 120
)
# expand the frequency table to 1 row per person
# (first just assign everyone the average Y value)
dat <- uncount(hw3_freq, N) %>%
# give those rows an id number
rowid_to_column(var = "id") %>%
# for each of the variables except for id, split it into two
# rows, one for each time point
pivot_longer(-id,
names_to = c(".value", "time"),
names_sep = "_"
) %>%
# make sure that time is read as a number
mutate(time = parse_number(time),
# add some random error to Y
# (nrow(.) means a different value for each row of the dataset in use
Y = Y + rnorm(nrow(.), 0, 10))
# create lagged variables
# group_by(id) makes this robust to row order -- lag() would otherwise
# silently lag across a different person's rows if dat were ever re-sorted
# (e.g. by time first). Currently harmless here since dat is in id-major
# order with exactly one row per id per time, but safer to be explicit.
t1 <- dat %>%
group_by(id) %>%
mutate(lag_A = lag(A)) %>%
ungroup() %>%
filter(time == 1)Then calculate the g-formula by hand.
The recipe is the same in both cases: fit a model for the time-1 covariate \(L\) given past treatment, fit a model for \(Y\) given treatment, past treatment and \(L\), then simulate 10,000 people forward with the treatments fixed at the regime’s values, drawing \(L\) from its fitted conditional distribution rather than conditioning on the observed one.
First for \(E[Y^{1,1}]\).
# fit models for covariate and outcome given past on t1
L_mod <- glm(L ~ lag_A, data = t1, family = binomial())
Y_mod <- glm(Y ~ A * lag_A * L, data = t1)
# create new data set with space for 10,000 simulated entries
nsim <- 10000
baseline_pop <- tibble(id = 1:nsim)
# fix intervention values in the new data set to be consistent with the regime
new_dat_11 <- baseline_pop %>%
mutate(lag_A = 1, A = 1)
# predict covariate means at time 1
new_dat_11$pL <- predict(L_mod, newdata = new_dat_11, type = "response")
# simulate covariate values at time 1
new_dat_11$L <- rbinom(nsim, 1, new_dat_11$pL)
# predict outcome at time 1 based on simulated covariates and fixed treatments
new_dat_11$Ey <- predict(Y_mod, newdata = new_dat_11, type = "response")
# take mean across all simulations to get g-formula estimate!
mean(new_dat_11$Ey)[1] 150.1166That gives \(E[Y^{1,1}] = 150.12\).
Then for \(E[Y^{0,0}]\).
# for 0, 0
new_dat_00 <- baseline_pop %>%
mutate(lag_A = 0, A = 0)
# predict covariate means at time 1
new_dat_00$pL <- predict(L_mod, newdata = new_dat_00, type = "response")
# simulate covariate values at time 1
new_dat_00$L <- rbinom(nsim, 1, new_dat_00$pL)
# predict outcome at time 1 based on simulated covariates and fixed treatments
new_dat_00$Ey <- predict(Y_mod, newdata = new_dat_00, type = "response")
# take mean across all simulations to get g-formula estimate!
mean(new_dat_00$Ey)[1] 134.7198And \(E[Y^{0,0}] = 134.72\), so always-treat versus never-treat is worth \(150.12 - 134.72 = 15.40\) by hand.
Now with the gfoRmula package.
# running the g-formula ---------------------------------------------------
# add your comments to the lines below
id <- "id"
time_name <- "time"
covnames <- c("L", "A")
outcome_name <- "Y"
#
outcome_type <- "continuous_eof"
#
histories <- c(lagged)
#
histvars <- list(c("A"))
#
covtypes <- c("binary", "binary")
covparams <- list(covmodels = c(
L ~ lag1_A,
A ~ lag1_A * L
))
ymodel <- Y ~ A * L * lag1_A
#
intvars <- list("A", "A", "A", "A")
#
interventions <- list(
list(c(static, c(0, 0))),
list(c(static, c(0, 1))),
list(c(static, c(1, 0))),
list(c(static, c(1, 1)))
)
int_descript <- c(
"0, 0",
"0, 1",
"1, 0",
"1, 1"
)
# put it all together
# (the warning "obs_data was coerced to a data table" is fine!)
gform_res <- gfoRmula::gformula(
obs_data = dat,
id = id,
time_name = time_name,
covnames = covnames,
outcome_name = outcome_name,
outcome_type = outcome_type,
covtypes = covtypes,
covparams = covparams,
ymodel = ymodel,
intvars = intvars,
interventions = interventions,
int_descript = int_descript,
ref_int = 1,
histories = histories,
histvars = histvars,
sim_data_b = TRUE,
seed = 345636
)
# Use the following code to explore the results. (You can ignore Intervention 0,
# the natural course intervention, for now.) How many models were fit? How can
# you tell which of the data in the `sim_data` object was simulated or predicted
# from a model? Does these results match your answers to the earlier questions,
gform_resPREDICTED RISK UNDER MULTIPLE INTERVENTIONS
Intervention Description
0 Natural course
1 0, 0
2 0, 1
3 1, 0
4 1, 1
Sample size = 32000, Monte Carlo sample size = 32000
Number of bootstrap samples = 0
Reference intervention = 1
k Interv. NP mean g-form mean Mean ratio Mean difference % Intervened On
<num> <num> <num> <num> <num> <num> <num>
1 0 142.4572 141.5895 1.0523929 7.0489627 0.00000
1 1 NA 134.5405 1.0000000 0.0000000 74.87812
1 2 NA 134.5276 0.9999039 -0.0129268 75.12188
1 3 NA 150.3494 1.1175029 15.8088974 80.85625
1 4 NA 150.0516 1.1152895 15.5111169 69.14375
Aver % Intervened On
<num>
0.00000
49.74531
50.25469
55.99375
44.00625gform_res$coeffs$L
(Intercept) lag1_A
8.985178e-14 1.098612e+00
$A
(Intercept) lag1_A L lag1_A:L
-1.098612e+00 -2.273084e-13 2.197225e+00 2.983558e-14
$Y
(Intercept) A L lag1_A A:L
59.97804299 0.07788309 149.84611288 180.24000339 -0.18249806
A:lag1_A L:lag1_A A:L:lag1_A
-0.37075696 -269.83091494 0.17594707 gform_res$sim_data$`Natural course`
id time L A eligible_pt intervened lag1_A Ey
<int> <num> <num> <num> <lgcl> <num> <num> <num>
1: 1 0 NA 0 TRUE 0 0 NA
2: 1 1 1 1 TRUE 0 0 209.7195
3: 2 0 NA 0 TRUE 0 0 NA
4: 2 1 1 1 TRUE 0 0 209.7195
5: 3 0 NA 0 TRUE 0 0 NA
---
63996: 31998 1 1 1 TRUE 0 1 119.9338
63997: 31999 0 NA 1 TRUE 0 0 NA
63998: 31999 1 1 1 TRUE 0 1 119.9338
63999: 32000 0 NA 1 TRUE 0 0 NA
64000: 32000 1 1 1 TRUE 0 1 119.9338
$`0, 0`
id time L A eligible_pt intervened lag1_A Ey
<int> <num> <num> <num> <lgcl> <lgcl> <num> <num>
1: 1 0 NA 0 TRUE FALSE 0 NA
2: 1 1 1 0 TRUE TRUE 0 209.82416
3: 2 0 NA 0 TRUE FALSE 0 NA
4: 2 1 1 0 TRUE TRUE 0 209.82416
5: 3 0 NA 0 TRUE FALSE 0 NA
---
63996: 31998 1 0 0 TRUE FALSE 0 59.97804
63997: 31999 0 NA 0 TRUE TRUE 0 NA
63998: 31999 1 0 0 TRUE FALSE 0 59.97804
63999: 32000 0 NA 0 TRUE TRUE 0 NA
64000: 32000 1 1 0 TRUE TRUE 0 209.82416
$`0, 1`
id time L A eligible_pt intervened lag1_A Ey
<int> <num> <num> <num> <lgcl> <lgcl> <num> <num>
1: 1 0 NA 0 TRUE FALSE 0 NA
2: 1 1 1 1 TRUE FALSE 0 209.71954
3: 2 0 NA 0 TRUE FALSE 0 NA
4: 2 1 1 1 TRUE FALSE 0 209.71954
5: 3 0 NA 0 TRUE FALSE 0 NA
---
63996: 31998 1 0 1 TRUE TRUE 0 60.05593
63997: 31999 0 NA 0 TRUE TRUE 0 NA
63998: 31999 1 0 1 TRUE TRUE 0 60.05593
63999: 32000 0 NA 0 TRUE TRUE 0 NA
64000: 32000 1 1 1 TRUE FALSE 0 209.71954
$`1, 0`
id time L A eligible_pt intervened lag1_A Ey
<int> <num> <num> <num> <lgcl> <lgcl> <num> <num>
1: 1 0 NA 1 TRUE TRUE 0 NA
2: 1 1 1 0 TRUE TRUE 1 120.2332
3: 2 0 NA 1 TRUE TRUE 0 NA
4: 2 1 1 0 TRUE TRUE 1 120.2332
5: 3 0 NA 1 TRUE TRUE 0 NA
---
63996: 31998 1 1 0 TRUE TRUE 1 120.2332
63997: 31999 0 NA 1 TRUE FALSE 0 NA
63998: 31999 1 1 0 TRUE TRUE 1 120.2332
63999: 32000 0 NA 1 TRUE FALSE 0 NA
64000: 32000 1 1 0 TRUE TRUE 1 120.2332
$`1, 1`
id time L A eligible_pt intervened lag1_A Ey
<int> <num> <num> <num> <lgcl> <lgcl> <num> <num>
1: 1 0 NA 1 TRUE TRUE 0 NA
2: 1 1 1 1 TRUE FALSE 1 119.9338
3: 2 0 NA 1 TRUE TRUE 0 NA
4: 2 1 1 1 TRUE FALSE 1 119.9338
5: 3 0 NA 1 TRUE TRUE 0 NA
---
63996: 31998 1 1 1 TRUE FALSE 1 119.9338
63997: 31999 0 NA 1 TRUE FALSE 0 NA
63998: 31999 1 1 1 TRUE FALSE 1 119.9338
63999: 32000 0 NA 1 TRUE FALSE 0 NA
64000: 32000 1 1 1 TRUE FALSE 1 119.9338The package agrees with the hand calculation. Intervention 4 (“1, 1”, always treat) gives a g-formula mean of 150.0516 against the 150.12 computed by hand, and intervention 1 (“0, 0”, never treat) gives 134.5405 against 134.72 — the small gaps are Monte Carlo error, since both routes simulate rather than integrate. The table’s own contrast of the two, 15.51, matches the 15.40 from the hand calculation.
The natural course is 141.59, sitting between the two regimes, and the non-parametric mean of the observed data is 142.46, which is the check that the model-based natural course is not badly misspecified.
The other two regimes are worth a glance. Treating only at time 0 gives 150.35 and treating only at time 1 gives 134.53 — so essentially all of the effect comes from the first period’s treatment, and the second period adds nothing.
Reading gform_res$coeffs shows why this needs the g-formula rather than a regression. The \(L\) model has a coefficient of 1.0986 on lagged treatment, so past treatment moves the time-1 covariate; and the \(Y\) model has 149.85 on \(L\) and \(-269.83\) on the \(L \times \text{lag}_1 A\) interaction, so that covariate in turn moves the outcome. \(L\) is a time-varying confounder affected by past treatment, which is exactly the configuration where conditioning on it blocks part of the effect and omitting it leaves confounding.
Chistopher Boyer has a package “gmethods” that can be used to estimate the g-formula. To me it’s a bit more intuitive than the gfoRmula package. Here is an example in the “gmethods” package, with the same model in the “gfoRmula” package.
The example dataset continuous_eofdata again consists of 7,500 observations on 2,500 individuals with a maximum of 7 follow-up times, where the outcome corresponds to a characteristic only in the last interval (e.g., systolic blood pressure in interval 7). The variables in the dataset are: t0: The time index id: A unique identifier for each individual L1: A time-varying covariate; categorical L2: A time-varying covariate; continuous L3: A baseline covariate; continuous A: The treatment variable; binary Y: The outcome of interest; continuous
The goal is to estimate the mean outcome at time \(K + 1 = 7\) under ‘‘never treat’’ versus ‘‘always treat.’’
The DAG would be something like:
dag <- dagify(
Y ~ A + L1 + L2 + L3 + lag1_A + lag1_L1 + lag1_L2 ,
L1 ~ lag1_A + lag1_L1 + lag1_L2 + L3 ,
L2 ~ lag1_A + L1 + lag1_L1 + lag1_L2 + L3,
A ~ lag1_A + L1 + L2 + lag1_L1 + lag1_L2 + L3,
coords = time_ordered_coords()
)
dag %>%
ggplot(aes(x = x, y = y, xend = xend, yend = yend)) +
geom_dag_point() +
geom_dag_edges_arc() +
geom_dag_text() +
theme_dag()
First use the gFormula package.
id <- 'id'
time_name <- 't0'
covnames <- c('L1', 'L2', 'A')
outcome_name <- 'Y'
covtypes <- c('categorical', 'normal', 'binary')
histories <- c(lagged)
histvars <- list(c('A', 'L1', 'L2'))
covparams <- list(covmodels = c(L1 ~ lag1_A + lag1_L1 + L3 + t0 +
lag1_L2,
L2 ~ lag1_A + L1 + lag1_L1 + lag1_L2 + L3 + t0,
A ~ lag1_A + L1 + L2 + lag1_L1 + lag1_L2 + L3 + t0))
ymodel <- Y ~ A + L1 + L2 + lag1_A + lag1_L1 + lag1_L2 + L3
intvars <- list('A', 'A')
interventions <- list(list(c(static, rep(0, 7))),
list(c(static, rep(1, 7))))
int_descript <- c('Never treat', 'Always treat')
nsimul <- 50000
gform_cont_eof <- gformula_continuous_eof(obs_data = continuous_eofdata,
id = id,
time_name = time_name,
covnames = covnames,
outcome_name = outcome_name,
covtypes = covtypes,
covparams = covparams, ymodel = ymodel,
intvars = intvars,
interventions = interventions,
int_descript = int_descript,
histories = histories,
histvars = histvars,
model_fits = TRUE,
basecovs = c("L3"),
nsimul = nsimul, seed = 1234)
gform_cont_eofPREDICTED RISK UNDER MULTIPLE INTERVENTIONS
Intervention Description
0 Natural course
1 Never treat
2 Always treat
Sample size = 2500, Monte Carlo sample size = 50000
Number of bootstrap samples = 0
Reference intervention = natural course (0)
k Interv. NP mean g-form mean Mean ratio Mean difference % Intervened On
<num> <num> <num> <num> <num> <num> <num>
6 0 -4.414543 -4.363027 1.0000000 0.0000000 0.000
6 1 NA -3.133081 0.7180981 1.2299454 100.000
6 2 NA -4.602952 1.0549904 -0.2399248 74.038
Aver % Intervened On
<num>
0.00000
82.35200
17.49571Now use the gmethods package.
gf <- gmethods::gformula(
outcome_model = list(
formula = Y ~ A + L1 + L2 + lag1_A + lag1_L1 + lag1_L2 + L3,
link = "identity",
family = "normal"
),
covariate_model = list(
"L1" = list(
formula = L1 ~ lag1_A + lag1_L1 + L3 + t0 + lag1_L2,
family = "categorical"
),
"L2" = list(
formula = L2 ~ lag1_A + L1 + lag1_L1 + lag1_L2 + L3 + t0,
link = "identity",
family = "normal"
),
"A" = list(
formula = A ~ lag1_A + L1 + L2 + lag1_L1 + lag1_L2 + L3 + t0,
link = "logit",
family = "binomial"
)
),
data = continuous_eofdata,
survival = FALSE,
id = 'id',
time = 't0'
)
s <- simulate(
gf,
interventions = list(
"Never treat" = function(data, time) set(data, j = "A", value = 0),
"Always treat" = function(data, time) set(data, j = "A", value = 1)
),
n_samples = 50000
)
s$results$means
intervention estimate conf.low conf.high
1 0 -4.360358 NA NA
2 1 -3.139260 NA NA
3 2 -4.602994 NA NA
$ratios
intervention estimate conf.low conf.high
1 0 NA NA NA
2 1 0.7199547 NA NA
3 2 1.0556459 NA NA
$diffs
intervention estimate conf.low conf.high
1 0 NA NA NA
2 1 1.221098 NA NA
3 2 -0.242636 NA NABoth packages agree closely. gfoRmula gives a natural course of \(-4.363\), never-treat \(-3.133\) and a mean difference of 1.230; gmethods gives \(-4.360\), \(-3.139\) and 1.221, with always-treat at \(-4.603\) and a difference of \(-0.243\). The two implementations differ in the third decimal, which is Monte Carlo noise from the 50,000-draw simulation, not a disagreement about the estimand.
Note the results are very close. Also note that in the covariates and treatment models, there is a time variable t0. This is because these models are pooled regression with all time points. Therefore a time fixed effect is included. In the outcome model, there is only the last period data, so the time variable is not included.
Standard errors can be calculated with bootstrapping.