Consider the observed data with time ordering W -> A -> Y. W is a continuous baseline covariate that affects A and Y. A is a binary treatment variable that affects Y. Y is a binary outcome variable.
W ~ N(0, 1)
A ~ binomial with P(A = 1) = expit(W)
Y ~ binomial with P(Y = 1) = expit(W + A)
where \(expit(z) = \frac{1}{1 + e^{-z}}\)
We want to know \(E[Y_1]\) the expected value of Y, intervening to set A to 1.
rexpit <- function(x) rbinom(n=length(x), size=1, prob=plogis(x))
n <- 10000
W <- rnorm(n)
A <- rexpit(W)
Y <- rexpit(W + A)
data <- data.frame(W, A, Y)
head(data)
#> W A Y
#> 1 -1.00576427 1 0
#> 2 -0.16776940 0 0
#> 3 -0.47677696 1 1
#> 4 -0.43573417 0 0
#> 5 0.35557205 0 1
#> 6 0.02219503 1 1We could try to use the simple mean of Y or the mean of Y when A = 1, but these will be biased.
Now we use ltmle().
result <- ltmle(data, Anodes = "A", Ynodes = "Y", abar = 1)
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W
#>
#> Estimate of time to completion: < 1 minute
result
#> Call:
#> ltmle(data = data, Anodes = "A", Ynodes = "Y", abar = 1)
#>
#> TMLE Estimate: 0.6954446Because we’re using simulated data, we can calculate the true value of \(E[Y_1]\)
Time ordering of data is W1 -> W2 -> W3 -> A -> Y
n <- 1000
W1 <- rnorm(n)
W2 <- rbinom(n, size=1, prob=0.3)
W3 <- rnorm(n)
A <- rexpit(-1 + 2 * W1 + W3)
Y <- rexpit(-0.5 + 2 * W1^2 + 0.5 * W2 - 0.5 * A + 0.2 * W3 * A - 1.1 * W3)
data <- data.frame(W1, W2, W3, A, Y)True value of \(E[Y_1]\) is approximately 0.5939.
SuperLearner semiparametric estimation using all parents as regressors
result <- ltmle(data, Anodes="A", Lnodes=NULL, Ynodes="Y", abar=1, SL.library="default")
#> Loading required namespace: SuperLearner
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W1 + W2 + W3 + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W1 + W2 + W3
#>
#> Loading required package: nnls
#> Loading required namespace: arm
#> Estimate of time to completion: < 1 minuteTMLE estimate:
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = "A", Lnodes = NULL, Ynodes = "Y",
#> abar = 1, SL.library = "default")
#>
#> Parameter Estimate: 0.62675
#> Estimated Std Err: 0.044853
#> p-value: <2e-16
#> 95% Conf Interval: (0.53884, 0.71466)IPTW estimate:
summary(result, estimator="iptw")
#> Estimator: iptw
#> Call:
#> ltmle(data = data, Anodes = "A", Lnodes = NULL, Ynodes = "Y",
#> abar = 1, SL.library = "default")
#>
#> Parameter Estimate: 0.61209
#> Estimated Std Err: 0.042995
#> p-value: <2e-16
#> 95% Conf Interval: (0.52782, 0.69636)SuperLearner semiparametric estimation using correctly specified regressors. This passes W1^2, W2, W3, A and W3:A as columns of matrix X to SuperLearner for the Q regression and W1 and W3 as columns of matrix X to SuperLearner for the g regression.
result <- ltmle(data, Anodes="A", Lnodes=NULL, Ynodes="Y",
Qform=c(Y="Q.kplus1 ~ I(W1^2) + W2 + W3*A"),
gform="A ~ W1 + W3", abar=1, SL.library="default")
#> Estimate of time to completion: 1 minute
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = "A", Lnodes = NULL, Ynodes = "Y",
#> Qform = c(Y = "Q.kplus1 ~ I(W1^2) + W2 + W3*A"), gform = "A ~ W1 + W3",
#> abar = 1, SL.library = "default")
#>
#> Parameter Estimate: 0.61225
#> Estimated Std Err: 0.032859
#> p-value: <2e-16
#> 95% Conf Interval: (0.54785, 0.67666)glm using correctly specified Qform and gform
result <- ltmle(data, Anodes="A", Lnodes=NULL, Ynodes="Y",
Qform=c(Y="Q.kplus1 ~ I(W1^2) + W2 + W3*A"), gform="A ~ W1 + W3",
abar=1, SL.library=NULL)
#> Estimate of time to completion: < 1 minute
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = "A", Lnodes = NULL, Ynodes = "Y",
#> Qform = c(Y = "Q.kplus1 ~ I(W1^2) + W2 + W3*A"), gform = "A ~ W1 + W3",
#> abar = 1, SL.library = NULL)
#>
#> Parameter Estimate: 0.61352
#> Estimated Std Err: 0.032575
#> p-value: <2e-16
#> 95% Conf Interval: (0.54967, 0.67736)Get summary measures (additive treatment effect, odds ratio, relative risk) for abar=1 vs abar=0
result <- ltmle(data, Anodes="A", Lnodes=NULL, Ynodes="Y",
abar=list(1, 0), SL.library="default")
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W1 + W2 + W3 + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W1 + W2 + W3
#>
#> Estimate of time to completion: < 1 minute
#> Warning in predict.lm(object, newdata, se.fit, scale = 1, type = if (type == :
#> prediction from a rank-deficient fit may be misleading
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = "A", Lnodes = NULL, Ynodes = "Y",
#> abar = list(1, 0), SL.library = "default")
#>
#> Treatment Estimate:
#> Parameter Estimate: 0.62675
#> Estimated Std Err: 0.044989
#> p-value: <2e-16
#> 95% Conf Interval: (0.53858, 0.71493)
#>
#> Control Estimate:
#> Parameter Estimate: 0.66871
#> Estimated Std Err: 0.028203
#> p-value: <2e-16
#> 95% Conf Interval: (0.61344, 0.72399)
#>
#> Additive Treatment Effect:
#> Parameter Estimate: -0.041957
#> Estimated Std Err: 0.053098
#> p-value: 0.42942
#> 95% Conf Interval: (-0.14603, 0.062113)
#>
#> Relative Risk:
#> Parameter Estimate: 0.93726
#> Est Std Err log(RR): 0.083254
#> p-value: 0.43638
#> 95% Conf Interval: (0.79614, 1.1034)
#>
#> Odds Ratio:
#> Parameter Estimate: 0.8319
#> Est Std Err log(OR): 0.23063
#> p-value: 0.42487
#> 95% Conf Interval: (0.52936, 1.3073)Time ordering of data is W -> C -> Y. C is a censoring node. The parameter of interest is Y, intervening to set C to be uncensored. Censoring nodes are very similar to treatment nodes. The main difference is that after censoring, other data is expected to be NA. It is assumed that the parameter of interest always intervenes to set censoring nodes to uncensored, so this is not specified in abar. g and Q regressions always stratify on censoring nodes, whereas g and Q regressions can either stratify or pool over treatment nodes (by using the stratify argument.)
Censoring nodes in data should be factors with two levels - “censored” and “uncensored”. The utility function BinaryToCensoring can be used to facilitate this.
n <- 100000
W <- rnorm(n)
C <- BinaryToCensoring(is.censored = rexpit(W))
summary(C)
#> censored uncensored
#> 49910 50090
Y <- rep(NA, n)
Y[C == "uncensored"] <- rexpit(W[C == "uncensored"])
data <- data.frame(W, C, Y)
head(data, 20)
#> W C Y
#> 1 -0.2987284 censored NA
#> 2 1.8769204 uncensored 1
#> 3 -0.9733347 uncensored 0
#> 4 -1.0332095 uncensored 1
#> 5 1.1169817 censored NA
#> 6 0.7518623 censored NA
#> 7 0.5282109 censored NA
#> 8 -0.6633313 uncensored 0
#> 9 -0.7219888 censored NA
#> 10 0.3909210 censored NA
#> 11 -0.4193815 censored NA
#> 12 -0.7178721 uncensored 0
#> 13 0.9187818 censored NA
#> 14 1.7185865 censored NA
#> 15 -0.4896426 censored NA
#> 16 1.2803184 censored NA
#> 17 0.3876463 censored NA
#> 18 0.1425409 censored NA
#> 19 -0.4874995 uncensored 0
#> 20 -0.2747150 censored NA
result <- ltmle(data, Anodes = NULL, Cnodes = "C", Ynodes = "Y", abar = NULL)
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W
#>
#> gform not specified, using defaults:
#> formula for C:
#> C ~ W
#>
#> Estimate of time to completion: < 1 minute
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = NULL, Cnodes = "C", Ynodes = "Y",
#> abar = NULL)
#>
#> Parameter Estimate: 0.50057
#> Estimated Std Err: 0.0023378
#> p-value: <2e-16
#> 95% Conf Interval: (0.49599, 0.50516)The naive estimate is biased (the true value is 0.5):
Time ordering of data is W -> A1 -> L -> A2 -> Y. L is time-dependent covariate because it occurs after a treatment (or censoring) variable, A1. We indicate this using the Lnodes argument.
n <- 1000
W <- rnorm(n)
A1 <- rexpit(W)
L <- 0.3 * W + 0.2 * A1 + rnorm(n)
A2 <- rexpit(W + A1 + L)
Y <- rexpit(W - 0.6 * A1 + L - 0.8 * A2)
data <- data.frame(W, A1, L, A2, Y)
head(data)
#> W A1 L A2 Y
#> 1 -1.2143715 0 -0.72291649 0 0
#> 2 0.2368600 1 0.57968819 1 0
#> 3 -0.7370930 0 0.35485803 0 1
#> 4 0.3061751 1 0.20143846 1 0
#> 5 -0.7264741 0 -0.08946743 0 0
#> 6 -1.2620376 0 0.12865984 0 0Treatment regime of interest: set A1 to 0, set A2 to 0
ltmle(data, Anodes=c("A1", "A2"), Lnodes="L", Ynodes="Y", abar=c(0, 0))
#> 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: < 1 minute
#> Call:
#> ltmle(data = data, Anodes = c("A1", "A2"), Lnodes = "L", Ynodes = "Y",
#> abar = c(0, 0))
#>
#> TMLE Estimate: 0.5235939W -> A1 -> C -> L -> A2 -> Y
n <- 1000
W <- rnorm(n)
A1 <- rexpit(W)
C <- BinaryToCensoring(is.censored = rexpit(0.6 * W - 0.5 * A1))
uncensored <- C == "uncensored"
L <- A2 <- Y <- rep(NA, n)
L[uncensored] <- (0.3 * W[uncensored] + 0.2 * A1[uncensored] + rnorm(sum(uncensored)))
A2[uncensored] <- rexpit(W[uncensored] + A1[uncensored] + L[uncensored])
Y[uncensored] <- rexpit(W[uncensored] - 0.6 * A1[uncensored] + L[uncensored] - 0.8 * A2[uncensored])
data <- data.frame(W, A1, C, L, A2, Y)
head(data)
#> W A1 C L A2 Y
#> 1 0.75424772 1 censored NA NA NA
#> 2 -0.58126094 1 uncensored -1.491608 0 0
#> 3 0.03526288 0 censored NA NA NA
#> 4 -1.56066276 0 censored NA NA NA
#> 5 0.66095190 1 censored NA NA NA
#> 6 -0.04221528 0 uncensored -1.907210 1 0Treatment regime of interest: set A1 to 1, set A2 to 0, set C to uncensored:
ltmle(data, Anodes=c("A1", "A2"), Cnodes = "C", Lnodes="L", Ynodes="Y", abar=c(1, 0))
#> 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 C:
#> C ~ W + A1
#> formula for A2:
#> A2 ~ W + A1 + L
#>
#> Estimate of time to completion: < 1 minute
#> Call:
#> ltmle(data = data, Anodes = c("A1", "A2"), Cnodes = "C", Lnodes = "L",
#> Ynodes = "Y", abar = c(1, 0))
#>
#> TMLE Estimate: 0.5521191Treatment regime of interest is: Always treat at time 1 (A1 = 1), treat at at time 2 (A2 = 1) if L > 0
abar <- matrix(nrow=n, ncol=2)
abar[, 1] <- 1
abar[, 2] <- L > 0
result.abar <- ltmle(data, Anodes=c("A1", "A2"), Cnodes = "C", Lnodes="L", Ynodes="Y", abar=abar)
#> 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 C:
#> C ~ W + A1
#> formula for A2:
#> A2 ~ W + A1 + L
#>
#> Estimate of time to completion: < 1 minute
result.abar
#> Call:
#> ltmle(data = data, Anodes = c("A1", "A2"), Cnodes = "C", Lnodes = "L",
#> Ynodes = "Y", abar = abar)
#>
#> TMLE Estimate: 0.4048778The regime can also be specified as a rule function. rule is a function applied to each row of data which returns a numeric vector of the same length as Anodes.
rule <- function(row) c(1, row["L"] > 0)
result.rule <- ltmle(data, Anodes=c("A1", "A2"), Cnodes = "C", Lnodes="L", Ynodes="Y", rule=rule)
#> 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 C:
#> C ~ W + A1
#> formula for A2:
#> A2 ~ W + A1 + L
#>
#> Estimate of time to completion: < 1 minute
result.rule
#> Call:
#> ltmle(data = data, Anodes = c("A1", "A2"), Cnodes = "C", Lnodes = "L",
#> Ynodes = "Y", rule = rule)
#>
#> TMLE Estimate: 0.4048778Specfifying the regime using abar and using rule give the same result:
summary(result.abar)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = c("A1", "A2"), Cnodes = "C", Lnodes = "L",
#> Ynodes = "Y", abar = abar)
#>
#> Parameter Estimate: 0.40488
#> Estimated Std Err: 0.042302
#> p-value: <2e-16
#> 95% Conf Interval: (0.32197, 0.48779)
summary(result.rule)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = c("A1", "A2"), Cnodes = "C", Lnodes = "L",
#> Ynodes = "Y", rule = rule)
#>
#> Parameter Estimate: 0.40488
#> Estimated Std Err: 0.042302
#> p-value: <2e-16
#> 95% Conf Interval: (0.32197, 0.48779)Consider a simple point treatment problem with observed data \(W, A, Y\). But there is a positivity problem - for small values of \(W\), \(Prob(A = 1)\) is very small.
The true parameter value, \(E[Y_1]\) is approximately 0.697.
The true TMLE standard deviation (the standard deviation of the TMLE estimate if we ran it many times on many sets of data) is approximately 0.056.
The true IPTW standard deviation (the standard deviation of the IPTW estimate if we ran it many times on many sets of data) is approximately 0.059.
The default variance.method is "tmle" - use TMLE in order to approximate the variance of the TMLE estimator. The estimated standard deviation is close to the true TMLE standard deviation.
r1 <- ltmle(df, Anodes="A", Ynodes="Y", abar = 1, estimate.time=FALSE)
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W
#>
print(summary(r1))
#> Estimator: tmle
#> Call:
#> ltmle(data = df, Anodes = "A", Ynodes = "Y", abar = 1, estimate.time = FALSE)
#>
#> Parameter Estimate: 0.67782
#> Estimated Std Err: 0.057989
#> p-value: <2e-16
#> 95% Conf Interval: (0.56417, 0.79148)If variance.method is "ic", variance is estimated using the estimated Influence Curve. This is fast to compute, but may be significantly anticonservative in data with positivity violations.
r2 <- ltmle(df, Anodes="A", Ynodes="Y", abar = 1, estimate.time=FALSE,
variance.method="ic")
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W
#>
#> Warning in CheckForVarianceWarning(inputs, g.ratio): Variance estimate is
#> based on influence curve only, which may be significantly anticonservative
#> because your data appears to contain positivity violations. It is recommended
#> to use variance.method='tmle' or variance.method='iptw' to obtain a more
#> robust variance estimate (but run time may be significantly longer). See
#> variance.method details in ?ltmle
print(summary(r2))
#> Estimator: tmle
#> Call:
#> ltmle(data = df, Anodes = "A", Ynodes = "Y", abar = 1, estimate.time = FALSE,
#> variance.method = "ic")
#>
#> Parameter Estimate: 0.67782
#> Estimated Std Err: 0.037042
#> p-value: <2e-16
#> 95% Conf Interval: (0.60522, 0.75042)If variance.method is "iptw", then use IPTW in order to approximate the variance of the TMLE estimator. This is faster to compute than variance.method = "tmle" but less accurate (and slower to compute than variance.method = "ic" but more accurate).
r3 <- ltmle(df, Anodes="A", Ynodes="Y", abar = 1, estimate.time=FALSE,
variance.method="iptw")
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W
#>
print(summary(r3))
#> Estimator: tmle
#> Call:
#> ltmle(data = df, Anodes = "A", Ynodes = "Y", abar = 1, estimate.time = FALSE,
#> variance.method = "iptw")
#>
#> Parameter Estimate: 0.67782
#> Estimated Std Err: 0.045018
#> p-value: <2e-16
#> 95% Conf Interval: (0.58959, 0.76606)If we use the IPTW estimator, variance.method does not change the estimated standard deviation (it only affects the estimated standard deviation of the TMLE estimator).
print(summary(r1, estimator="iptw"))
#> Estimator: iptw
#> Call:
#> ltmle(data = df, Anodes = "A", Ynodes = "Y", abar = 1, estimate.time = FALSE)
#>
#> Parameter Estimate: 0.73313
#> Estimated Std Err: 0.047001
#> p-value: <2e-16
#> 95% Conf Interval: (0.64101, 0.82525)
print(summary(r2, estimator="iptw")) #the same - variance.method only affects TMLE
#> Estimator: iptw
#> Call:
#> ltmle(data = df, Anodes = "A", Ynodes = "Y", abar = 1, estimate.time = FALSE,
#> variance.method = "ic")
#>
#> Parameter Estimate: 0.73313
#> Estimated Std Err: 0.047001
#> p-value: <2e-16
#> 95% Conf Interval: (0.64101, 0.82525)
print(summary(r3, estimator="iptw")) #the same - variance.method only affects TMLE
#> Estimator: iptw
#> Call:
#> ltmle(data = df, Anodes = "A", Ynodes = "Y", abar = 1, estimate.time = FALSE,
#> variance.method = "iptw")
#>
#> Parameter Estimate: 0.73313
#> Estimated Std Err: 0.047001
#> p-value: <2e-16
#> 95% Conf Interval: (0.64101, 0.82525)We can see that the values of g are very small.
summary(r1$cum.g)
#> V1
#> Min. :0.01000
#> 1st Qu.:0.06621
#> Median :0.53923
#> Mean :0.52195
#> 3rd Qu.:0.95837
#> Max. :1.00000
summary(r1$cum.g.unbounded)
#> V1
#> Min. :0.0000019
#> 1st Qu.:0.0662084
#> Median :0.5392299
#> Mean :0.5210000
#> 3rd Qu.:0.9583683
#> Max. :0.9999999
head(data.frame(df, g = r1$cum.g, g.unbounded = r1$cum.g.unbounded), 20)
#> W A Y g g.unbounded
#> 1 -1.35648823 0 0 0.0100000 3.409752e-03
#> 2 1.26884814 1 1 0.9963765 9.963765e-01
#> 3 -0.37773620 1 0 0.1873937 1.873937e-01
#> 4 0.94972013 1 1 0.9858499 9.858499e-01
#> 5 -0.24114949 0 0 0.2932962 2.932962e-01
#> 6 1.21533567 1 1 0.9954428 9.954428e-01
#> 7 0.14571385 0 0 0.6867324 6.867324e-01
#> 8 0.40461617 1 1 0.8697448 8.697448e-01
#> 9 0.84557785 1 1 0.9780277 9.780277e-01
#> 10 -0.02310326 1 1 0.5146515 5.146515e-01
#> 11 0.37353042 1 1 0.8538312 8.538312e-01
#> 12 -2.46671406 0 1 0.0100000 2.883257e-05
#> 13 -0.42082160 0 0 0.1607860 1.607860e-01
#> 14 -0.13043009 0 1 0.4005665 4.005665e-01
#> 15 -0.26469350 0 0 0.2727496 2.727496e-01
#> 16 1.69489974 1 1 0.9994187 9.994187e-01
#> 17 0.13098866 1 1 0.6729456 6.729456e-01
#> 18 0.47752803 1 1 0.9013555 9.013555e-01
#> 19 -1.29563313 0 0 0.0100000 4.425675e-03
#> 20 1.06440181 1 1 0.9913126 9.913126e-01The id argument can be used to specify hierarchical data, such as people in a household.
num.households <- 500
people.in.household <- round(runif(num.households, min = 1, max = 10))
length(people.in.household)
#> [1] 500
n <- sum(people.in.household)
n
#> [1] 2738
W.household <- rnorm(num.households)
length(W.household)
#> [1] 500
W.household.expanded <- rep(W.household, times = people.in.household)
W.indiv <- rnorm(n)
length(W.indiv)
#> [1] 2738
A <- rexpit(1.5 * W.household.expanded + 0.4 * W.indiv)
Y <- rexpit(-1 + 2.3 * W.household.expanded - 0.6 * W.indiv + 1.2 * A)id can be an integer, factor, or character (or any type that can be coerced to factor),
id <- 1:num.households
id.expanded <- rep(id, times = people.in.household)
data <- data.frame(W.household.expanded, W.indiv, A, Y)
head(cbind(id.expanded, data), 20)
#> id.expanded W.household.expanded W.indiv A Y
#> 1 1 0.3539016 -0.539311380 1 1
#> 2 1 0.3539016 0.276729241 0 1
#> 3 1 0.3539016 0.254790709 1 0
#> 4 2 -2.2184914 -1.647148651 0 0
#> 5 2 -2.2184914 -0.007023282 0 0
#> 6 3 0.1240071 0.309937561 0 0
#> 7 3 0.1240071 0.500161573 0 0
#> 8 3 0.1240071 0.057252323 1 0
#> 9 3 0.1240071 -1.227004677 0 1
#> 10 3 0.1240071 0.741631195 0 0
#> 11 4 -1.6185501 0.529298479 0 0
#> 12 4 -1.6185501 1.498849852 0 0
#> 13 4 -1.6185501 -1.353555069 0 0
#> 14 4 -1.6185501 1.808790572 0 0
#> 15 4 -1.6185501 -0.376653967 0 0
#> 16 4 -1.6185501 0.672049220 1 0
#> 17 5 0.1735576 -0.641476699 0 0
#> 18 5 0.1735576 -1.405063679 1 0
#> 19 5 0.1735576 -1.268529841 0 1
#> 20 5 0.1735576 -0.946288370 0 1result.without.id <- ltmle(data, Anodes = "A", Ynodes = "Y", abar = 0)
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W.household.expanded + W.indiv + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W.household.expanded + W.indiv
#>
#> Estimate of time to completion: < 1 minute
result.with.id <- ltmle(data, Anodes = "A", Ynodes = "Y", abar = 0, id = id.expanded)
#> Qform not specified, using defaults:
#> formula for Y:
#> Q.kplus1 ~ W.household.expanded + W.indiv + A
#>
#> gform not specified, using defaults:
#> formula for A:
#> A ~ W.household.expanded + W.indiv
#>
#> Estimate of time to completion: < 1 minuteOmitting the id argument makes the individuals seem more independent than they are, which gives artificially low variance estimates.
summary(result.without.id)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = "A", Ynodes = "Y", abar = 0)
#>
#> Parameter Estimate: 0.36207
#> Estimated Std Err: 0.012771
#> p-value: <2e-16
#> 95% Conf Interval: (0.33704, 0.38711)
summary(result.with.id)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = "A", Ynodes = "Y", abar = 0, id = id.expanded)
#>
#> Parameter Estimate: 0.36207
#> Estimated Std Err: 0.019685
#> p-value: <2e-16
#> 95% Conf Interval: (0.32349, 0.40066)The influence curve is a vector with length equal to the number of independent units.
age -> gender -> A1 -> L1a -> L1b -> Y1 -> A2 -> L2a -> L2b -> Y2
n <- 1000
age <- rbinom(n, 1, 0.5)
gender <- rbinom(n, 1, 0.5)
A1 <- rexpit(age + gender)
L1a <- 2*age - 3*gender + 2*A1 + rnorm(n)
L1b <- rexpit(age + 1.5*gender - A1)
Y1 <- plogis(age - gender + L1a + 0.7*L1b + A1 + rnorm(n))
A2 <- rexpit(age + gender + A1 - L1a - L1b)
L2a <- 2*age - 3*gender + 2*A1 + A2 + rnorm(n)
L2b <- rexpit(age + 1.5*gender - A1 - A2)
Y2 <- plogis(age - gender + L1a + L1b + A1 + 1.8*A2 + rnorm(n))
data <- data.frame(age, gender, A1, L1a, L1b, Y1, A2, L2a, L2b, Y2)Also show some different ways of specifying the nodes (either names or indexes works):
result <- ltmle(data, Anodes=c(3, 7), Lnodes=c("L1a", "L1b", "L2a", "L2b"),
Ynodes=grep("^Y", names(data)), abar=c(1, 0))
#> Qform not specified, using defaults:
#> formula for L1a:
#> Q.kplus1 ~ age + gender + A1
#> formula for L2a:
#> Q.kplus1 ~ age + gender + A1 + L1a + L1b + Y1 + A2
#>
#> gform not specified, using defaults:
#> formula for A1:
#> A1 ~ age + gender
#> formula for A2:
#> A2 ~ age + gender + A1 + L1a + L1b + Y1
#>
#> Estimate of time to completion: < 1 minute
#> Warning in CheckForVarianceWarning(inputs, g.ratio): Variance estimate is based
#> on influence curve only, which may be significantly anticonservative because
#> your data appears to contain positivity violations. Robust variance estimate is
#> not currently available with non binary outcomes but this will be addressed in a
#> future release.
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = c(3, 7), Lnodes = c("L1a", "L1b",
#> "L2a", "L2b"), Ynodes = grep("^Y", names(data)), abar = c(1,
#> 0))
#>
#> Parameter Estimate: 0.82404
#> Estimated Std Err: 0.012936
#> p-value: <2e-16
#> 95% Conf Interval: (0.79869, 0.8494)Usually you would specify a Qform for all of the Lnodes and Ynodes but in this case L1a, L1b, Y1 is a “block” of L/Y nodes not separated by Anodes or Cnodes (the same is true for L2a, L2b, Y2). Only one regression is required at the first L/Y node in a block. You can pass regression formulas for the other L/Y nodes, but they will be ignored.
result <- ltmle(data, Anodes=c(3, 7), Lnodes=c("L1a", "L1b", "L2a", "L2b"),
Ynodes=grep("^Y", names(data)), abar=c(1, 0),
Qform=c(L1a="Q.kplus1 ~ 1", L2a="Q.kplus1 ~ 1"))
#> gform not specified, using defaults:
#> formula for A1:
#> A1 ~ age + gender
#> formula for A2:
#> A2 ~ age + gender + A1 + L1a + L1b + Y1
#>
#> Estimate of time to completion: < 1 minute
#> Warning in CheckForVarianceWarning(inputs, g.ratio): Variance estimate is based
#> on influence curve only, which may be significantly anticonservative because
#> your data appears to contain positivity violations. Robust variance estimate is
#> not currently available with non binary outcomes but this will be addressed in a
#> future release.
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = c(3, 7), Lnodes = c("L1a", "L1b",
#> "L2a", "L2b"), Ynodes = grep("^Y", names(data)), Qform = c(L1a = "Q.kplus1 ~ 1",
#> L2a = "Q.kplus1 ~ 1"), abar = c(1, 0))
#>
#> Parameter Estimate: 0.82977
#> Estimated Std Err: 0.030899
#> p-value: <2e-16
#> 95% Conf Interval: (0.76921, 0.89033)Gives the same result but prints a message saying some regression formulas will be dropped:
result <- ltmle(data, Anodes=c(3, 7), Lnodes=c("L1a", "L1b", "L2a", "L2b"),
Ynodes=grep("^Y", names(data)), abar=c(1, 0),
Qform=c(L1a="Q.kplus1 ~ 1", L1b="Q.klus1~A1",
Y1="Q.kplus1~L1a", L2a="Q.kplus1 ~ 1", L2b="Q.klus1~A1", Y2="Q.kplus1~A2 + gender"))
#> L/Y nodes (after removing blocks) : L1a L2a
#> Qform names : L1a L1b Y1 L2a L2b Y2
#> The following nodes are not being considered as L/Y nodes because they are part of a block
#> of L/Y nodes. They are being dropped from Qform:
#> L1b
#> Y1
#> L2b
#> Y2
#> gform not specified, using defaults:
#> formula for A1:
#> A1 ~ age + gender
#> formula for A2:
#> A2 ~ age + gender + A1 + L1a + L1b + Y1
#>
#> Estimate of time to completion: < 1 minute
#> Warning in CheckForVarianceWarning(inputs, g.ratio): Variance estimate is based
#> on influence curve only, which may be significantly anticonservative because
#> your data appears to contain positivity violations. Robust variance estimate is
#> not currently available with non binary outcomes but this will be addressed in a
#> future release.
summary(result)
#> Estimator: tmle
#> Call:
#> ltmle(data = data, Anodes = c(3, 7), Lnodes = c("L1a", "L1b",
#> "L2a", "L2b"), Ynodes = grep("^Y", names(data)), Qform = c(L1a = "Q.kplus1 ~ 1",
#> L1b = "Q.klus1~A1", Y1 = "Q.kplus1~L1a", L2a = "Q.kplus1 ~ 1",
#> L2b = "Q.klus1~A1", Y2 = "Q.kplus1~A2 + gender"), abar = c(1,
#> 0))
#>
#> Parameter Estimate: 0.82977
#> Estimated Std Err: 0.030899
#> p-value: <2e-16
#> 95% Conf Interval: (0.76921, 0.89033)If there were a Anode or Cnode between L1b and Y1, Y1 would also need a Q regression formula.