Marginal Structural Models - ltmleMSM()

Marginal Structural Models (MSMs) - Dynamic Treatment

In this example there are two treatment nodes and one outcome: W A1 L A2 Y W is normally distributed and L is continuous in (0, 1). We are interested in treatments where A1 is set to either 0 or 1 and A2 is set dynamically. The treatment for A2 is indexed by theta between 0 and 1. If $$L > theta$$, set A2 to 1, otherwise set A2 to 0.
Here is a function that can be used to generate observed data (if abar = NULL) or generate counterfactual truth (if abar is a list with a1 and theta):

GenerateData <- function(n, abar = NULL) {
W <- rnorm(n)
if (is.null(abar)) {
A1 <- rexpit(W)
} else {
A1 <- abar$a1 } L <- plogis(rnorm(n) + 0.3 * W + 0.5 * A1) if (is.null(abar)) { A2 <- rexpit(-0.5 * W + A1 - 0.6 * L) } else { A2 <- as.integer(L > abar$theta)
}
Y <- rexpit(-2 + W + A1 + L + 2 * A2)
if (is.null(abar)) {
return(data.frame(W, A1, L, A2, Y))
} else {
return(mean(Y))
}
}

Set up regimes and summary.measures:

set.seed(11)
n <- 10000
data <- GenerateData(n)
regimes <- array(dim = c(n, 2, 10)) #n x num.Anodes x num.regimes
theta.set <- seq(0, 1, length.out = 5)
summary.measures <- array(theta.set, dim = c(10, 2, 1))
colnames(summary.measures) <- c("a1", "theta")
cnt <- 0
for (a1 in 0:1) {
for (theta.index in 1:5) {
cnt <- cnt + 1
regimes[, 1, cnt] <- a1
regimes[, 2, cnt] <- data\$L > theta.set[theta.index]
summary.measures[cnt, , 1] <- c(a1, theta.set[theta.index])
}
}
summary.measures
#> , , 1
#>
#>       a1 theta
#>  [1,]  0  0.00
#>  [2,]  0  0.25
#>  [3,]  0  0.50
#>  [4,]  0  0.75
#>  [5,]  0  1.00
#>  [6,]  1  0.00
#>  [7,]  1  0.25
#>  [8,]  1  0.50
#>  [9,]  1  0.75
#> [10,]  1  1.00
#>             W A1          L A2 Y
#> 1 -0.59103110  1 0.77137613  1 1
#> 2  0.02659437  0 0.47561600  1 0
#> 3 -1.51655310  0 0.08392365  1 0
regimes[1:3, , ]
#> , , 1
#>
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    0    1
#> [3,]    0    1
#>
#> , , 2
#>
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    0    1
#> [3,]    0    0
#>
#> , , 3
#>
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    0    0
#> [3,]    0    0
#>
#> , , 4
#>
#>      [,1] [,2]
#> [1,]    0    1
#> [2,]    0    0
#> [3,]    0    0
#>
#> , , 5
#>
#>      [,1] [,2]
#> [1,]    0    0
#> [2,]    0    0
#> [3,]    0    0
#>
#> , , 6
#>
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    1    1
#> [3,]    1    1
#>
#> , , 7
#>
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    1    1
#> [3,]    1    0
#>
#> , , 8
#>
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    1    0
#> [3,]    1    0
#>
#> , , 9
#>
#>      [,1] [,2]
#> [1,]    1    1
#> [2,]    1    0
#> [3,]    1    0
#>
#> , , 10
#>
#>      [,1] [,2]
#> [1,]    1    0
#> [2,]    1    0
#> [3,]    1    0
working.msm <- "Y ~ a1*theta"
summary(ltmleMSM(data, Anodes = c("A1", "A2"), Lnodes = "L", Ynodes = "Y",
regimes = regimes, summary.measures = summary.measures,
working.msm = working.msm))
#> Qform not specified, using defaults:
#> formula for L:
#> Q.kplus1 ~ W + A1
#> formula for Y:
#> Q.kplus1 ~ W + A1 + L + A2
#>
#> gform not specified, using defaults:
#> formula for A1:
#> A1 ~ W
#> formula for A2:
#> A2 ~ W + A1 + L
#>
#> Estimate of time to completion: 2 to 3 minutes
#> Estimator:  tmle
#>              Estimate Std. Error   CI 2.5% CI 97.5% p-value
#> (Intercept)  0.531670   0.050294  0.433095    0.630  <2e-16 ***
#> a1           1.039499   0.074409  0.893660    1.185  <2e-16 ***
#> theta       -1.817514   0.071766 -1.958172   -1.677  <2e-16 ***
#> a1:theta    -0.006772   0.110068 -0.222502    0.209   0.951
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Let’s compare to the true coefficients of the MSM. First we find the true value of $$E[Y_{a1, theta}]$$ for 5 values of theta.

truth <- rep(NA_real_, 10)
cnt <- 0
for (a1 in 0:1) {
for (theta.index in 1:5) {
cnt <- cnt + 1
truth[cnt] <- GenerateData(n = 1e6,
abar = list(a1 = a1, theta = theta.set[theta.index]))
}
}

Fit a working MSM to the true values of $$E[Y_{a1, theta}]$$.

m.true <- glm(working.msm,
data = data.frame(Y = truth, summary.measures[, , 1]),
family = "quasibinomial")
m.true
#>
#> Call:  glm(formula = working.msm, family = "quasibinomial", data = data.frame(Y = truth,
#>     summary.measures[, , 1]))
#>
#> Coefficients:
#> (Intercept)           a1        theta     a1:theta
#>     0.51497      0.95904     -1.76125     -0.02284
#>
#> Degrees of Freedom: 9 Total (i.e. Null);  6 Residual
#> Null Deviance:       1.341
#> Residual Deviance: 0.02342   AIC: NA

The estimated MSM coefficients are close to the true MSM coefficients.