If you edit this Rmd, re-knit and refresh every mirror: this folder’s
ltmle-doodle.htmland zip, and the demo snapshots (copy the freshly writtenltmle-results.jsonoverresults.jsonin roadmap’ssrc/pages/analytics/components/LtmleDoodleDemo/and statistics’components/demo/LtmleDoodleDemo/). The statistics repo also mirrors this whole kit folder.
library(ltmle)
library(jsonlite)Two-stage longitudinal data with one binary treatment per stage and
one terminal binary outcome, one observed record per subject:
O = (L1, A1, L2, A2, Y)
get_data <- function(n) {
L1 <- rnorm(n)
A1 <- rbinom(n, 1, plogis(0.01 * L1))
L2 <- 0.5 * L1 + 0.7 * A1 + rnorm(n)
A2 <- rbinom(n, 1, plogis(0.01 * (L2 + A1 + L1)))
Y <- rbinom(n, 1, plogis(3 * L1 * A1 + 3 * L2 * A2))
data.frame(L1, A1, L2, A2, Y)
}
get_data_d <- function(n, abar = 1) {
L1 <- rnorm(n)
A1 <- if (identical(abar, "DTR")) as.numeric(L1 < 0) else abar[1]
L2 <- 0.5 * L1 + 0.7 * A1 + rnorm(n)
A2 <- if (identical(abar, "DTR")) as.numeric(L2 < 0) else abar[length(abar)]
Y <- rbinom(n, 1, plogis(3 * L1 * A1 + 3 * L2 * A2))
data.frame(L1, A1, L2, A2, Y)
}csv_path <- "ltmle-doodle.csv"
if (!file.exists(csv_path)) {
set.seed(123)
df <- get_data(1000)
write.csv(df, csv_path, row.names = FALSE)
} else {
df <- read.csv(csv_path)
}
head(df)## L1 A1 L2 A2 Y
## 1 -0.56047565 0 -1.1012245 0 0
## 2 -0.23017749 0 -0.4223460 1 0
## 3 1.55870831 1 0.5772561 0 1
## 4 0.07050839 0 0.6623229 0 1
## 5 0.12928774 1 1.8849989 1 1
## 6 1.71506499 0 2.9847460 1 1The naive estimator filters observed data to rows where the regime
was followed and takes the empirical mean of Y. It ignores
time-varying confounding by L2.
naive_static <- function(df, a1, a2) mean(df$Y[df$A1 == a1 & df$A2 == a2])
naive_dtr <- function(df) {
match_d <- df$A1 == as.numeric(df$L1 < 0) & df$A2 == as.numeric(df$L2 < 0)
mean(df$Y[match_d])
}
naive <- list(
always_on = naive_static(df, 1, 1),
always_off = naive_static(df, 0, 0),
dtr = naive_dtr(df)
)
naive## $always_on
## [1] 0.6798419
##
## $always_off
## [1] 0.4703833
##
## $dtr
## [1] 0.2741313The kit uses SL.library = c("SL.glm") for a sub-2-minute
runtime. To match the Session 14b3 doodle exactly, swap in the commented
SL.library line below; expect a 2-5 minute runtime.
sl_lib <- c("SL.glm")
# sl_lib <- c("SL.glm", "SL.earth", "SL.xgboost")
set.seed(123)
res_on <- ltmle(df,
Anodes = c("A1", "A2"), Lnodes = c("L1", "L2"), Ynodes = "Y",
abar = c(1, 1), SL.library = sl_lib, SL.cvControl = list(V = 10))
sum_on <- summary(res_on)
sum_on
set.seed(123)
res_off <- ltmle(df,
Anodes = c("A1", "A2"), Lnodes = c("L1", "L2"), Ynodes = "Y",
abar = c(0, 0), SL.library = sl_lib, SL.cvControl = list(V = 10))
sum_off <- summary(res_off)
sum_off
set.seed(123)
d <- cbind(df$L1 < 0, df$L2 < 0)
res_dtr <- ltmle(df,
Anodes = c("A1", "A2"), Lnodes = c("L1", "L2"), Ynodes = "Y",
abar = d, SL.library = sl_lib, SL.cvControl = list(V = 10))
sum_dtr <- summary(res_dtr)
sum_dtr## Estimator: tmle
## Call:
## ltmle(data = df, Anodes = c("A1", "A2"), Lnodes = c("L1", "L2"),
## Ynodes = "Y", abar = c(1, 1), SL.library = sl_lib, SL.cvControl = list(V = 10))
##
## Parameter Estimate: 0.66119
## Estimated Std Err: 0.027388
## p-value: <2e-16
## 95% Conf Interval: (0.60751, 0.71487)
##
## Estimator: tmle
## Call:
## ltmle(data = df, Anodes = c("A1", "A2"), Lnodes = c("L1", "L2"),
## Ynodes = "Y", abar = c(0, 0), SL.library = sl_lib, SL.cvControl = list(V = 10))
##
## Parameter Estimate: 0.46011
## Estimated Std Err: 0.031603
## p-value: <2e-16
## 95% Conf Interval: (0.39817, 0.52205)
##
## Estimator: tmle
## Call:
## ltmle(data = df, Anodes = c("A1", "A2"), Lnodes = c("L1", "L2"),
## Ynodes = "Y", abar = d, SL.library = sl_lib, SL.cvControl = list(V = 10))
##
## Parameter Estimate: 0.25559
## Estimated Std Err: 0.027516
## p-value: <2e-16
## 95% Conf Interval: (0.20166, 0.30952)Because we control the DGP we can read off the true counterfactual
means by simulating from get_data_d at large
n.
set.seed(999)
oracle <- list(
always_on = mean(get_data_d(1e6, abar = c(1, 1))$Y),
always_off = mean(get_data_d(1e6, abar = c(0, 0))$Y),
dtr = mean(get_data_d(1e6, abar = "DTR")$Y)
)
oracle## $always_on
## [1] 0.643711
##
## $always_off
## [1] 0.500445
##
## $dtr
## [1] 0.261327This chunk writes the JSON the analytics page reads to display the
LTMLE bar. The page recomputes oracle and naive live in the browser; the
LTMLE point estimate is a snapshot from this kit at
seed=123, n=1000.
extract_estimate <- function(s) {
est <- s$treatment$estimate
ci <- s$treatment$CI
list(est = unname(est), ci_low = unname(ci[1]), ci_high = unname(ci[2]))
}
results <- list(
seed = 123,
n = nrow(df),
sl_library = sl_lib,
naive = naive,
ltmle = list(
always_on = extract_estimate(sum_on),
always_off = extract_estimate(sum_off),
dtr = extract_estimate(sum_dtr)
),
oracle = oracle
)
write_json(results, "ltmle-results.json", auto_unbox = TRUE, digits = 6, pretty = TRUE)
results## $seed
## [1] 123
##
## $n
## [1] 1000
##
## $sl_library
## [1] "SL.glm"
##
## $naive
## $naive$always_on
## [1] 0.6798419
##
## $naive$always_off
## [1] 0.4703833
##
## $naive$dtr
## [1] 0.2741313
##
##
## $ltmle
## $ltmle$always_on
## $ltmle$always_on$est
## [1] 0.6611929
##
## $ltmle$always_on$ci_low
## [1] 0.6075133
##
## $ltmle$always_on$ci_high
## [1] 0.7148725
##
##
## $ltmle$always_off
## $ltmle$always_off$est
## [1] 0.4601124
##
## $ltmle$always_off$ci_low
## [1] 0.3981709
##
## $ltmle$always_off$ci_high
## [1] 0.5220539
##
##
## $ltmle$dtr
## $ltmle$dtr$est
## [1] 0.2555881
##
## $ltmle$dtr$ci_low
## [1] 0.2016575
##
## $ltmle$dtr$ci_high
## [1] 0.3095186
##
##
##
## $oracle
## $oracle$always_on
## [1] 0.643711
##
## $oracle$always_off
## [1] 0.500445
##
## $oracle$dtr
## [1] 0.261327