inst/example_SDE_glm.R

In jlstiles/SDE_transport: SDEtransportsim

devtools::install_github("jeremyrcoyle/sl3")
devtools::install_github("jlstiles/SDE_transport")
library("SDEtransport")
library(sl3)

#Set up data generating process:
# data generating process for 1-d W
f_W = function(n) {
  W1 = rnorm(n)
  W2 = rbinom(n, 1, 0.72)
  data.frame(W1=W1, W2=W2)
}
f_S = function(W) {
  with(W, plogis(W1 - W2 + 0.3))
}
# make a pscore model
f_A = function(S,W) {
  df = as.data.frame(cbind(S=S, W))
  with(df, plogis(-.6*S-.7*W1 + 0.1*W2 +.17))
}
# make a intermediate confounder model
f_Z = function(A,S,W) {
  df = as.data.frame(cbind(S=S, W, A = A))
  with(df, plogis(.1*S-.4*W1 - W2 + 1*A-.3))
}
# make an M model according to the restrictions
f_M = function(Z,W,S) {
  df = as.data.frame(cbind(Z=Z, W, S=S))
  with(df, plogis(-.14*S + 1*W1 - W2 + 1.2*Z +.1))
}
# make a Y model according to the restrictions, main terms linear logistic reg.
# plug-in is biased and not robust like tmle
f_Y = function(M,Z,W) {
  df = as.data.frame(cbind(M = M, Z = Z, W))
  with(df, plogis(W2*M + 3*cos(W1)*Z-.4))
}

# generate n random samples
n = 1e2
set.seed(1)
data = gendata.SDEtransport(n, f_W = f_W, f_S = f_S, f_A = f_A, f_Z = f_Z, f_M = f_M, f_Y = f_Y)
mean(data$S==0)

# define the forms for all the regressions
forms=list(Sform = formula("S ~ W1"),
           Aform = formula("A ~ S + W2"),
           Zstarform = formula("Z ~ S + A + W1"),
           QZform = formula("Qstar_Mg ~ Z + S + W1"),
           Mstarform = formula("M ~ Z + W2"),
           Yform = formula("Y ~ M + Z + W1"))

# generate n random samples
n = 5e3
set.seed(1)
data = gendata.SDEtransport(n, f_W = f_W, f_S = f_S, f_A = f_A, f_Z = f_Z, f_M = f_M, f_Y = f_Y)
mean(data$S==0)

W = data[,grep("W", colnames(data))]
data = list(W=W, S=data$S, A=data$A, Z=data$Z, M=data$M, Y=data$Y)

# run the tmle, specifying 2 bootstraps for bootstrap inference but obviously 2 is way too few.
testres = SDE_tmle_glm(data=data, truth = NULL, truncate = list(lower =.0001, upper = .9999), 
                        B = 2, forms, RCT = NULL) 

testres

# run the tmle, with a known data generating process, in case of simulation
truth = list(f_W = f_W, f_S = f_S, f_A = f_A, f_Z = f_Z, f_M = f_M, f_Y = f_Y)
testres1 = SDE_tmle_glm(data=data, truth = truth, truncate = list(lower =.0001, upper = .9999), 
                       B = 2, forms, RCT = NULL) 

testres1