tmle_dynamic_ate.R

tmle_dynamic_ate <- function(logitQnA1, logitQnA0, gn1, A, Y,  d1=1, d0=0, Q_is_logit=TRUE, gbound=0.01){
  n <- length(Y)
  if (length(d1) == 1) d1 <- rep(d1, n)
  if (length(d0) == 1) d0 <- rep(d0, n)
  if (length(logitQnA1) != n || 
        length(logitQnA0) != n || 
        length(gn1) != n || 
        length(d1) != n ||
        length(d0) != n ||
        length(A) != n) stop("inputs should all have the same length") 
  if (!Q_is_logit) {
    logitQnA1 <- qlogis(logitQnA1)
    logitQnA0 <- qlogis(logitQnA0)
  }
  
  gn1 <- pmin(1.0-gbound, pmax(gbound, gn1))
  
  gnd1 <- ifelse(d1==1, gn1, 1-gn1)
  gnd0 <- ifelse(d0==1, gn1, 1-gn1)
  
  logitQnAd1 <- ifelse(d1==1, logitQnA1, logitQnA0)
  logitQnAd0 <- ifelse(d0==1, logitQnA1, logitQnA0)
  logitQnAA <- ifelse(A ==1.0, logitQnA1, logitQnA0)
  
  hAA  <-  (A==d1)/gnd1 -  (A==d0)/gnd0
  hAd1 <-      1.0/gnd1 - (d1==d0)/gnd0
  hAd0 <- (d0==d1)/gnd1 -      1.0/gnd0
  
  Qnstar <- glm(Y~h-1, data.frame(h=hAA, off=logitQnAA), family=binomial, offset=off)
  
  QnstarAd1 <- predict(Qnstar, newdata=data.frame(h=hAd1, off=logitQnAd1), type="response")
  QnstarAd0 <- predict(Qnstar, newdata=data.frame(h=hAd0, off=logitQnAd0), type="response")
  QnstarAA  <- predict(Qnstar, type="response")
  
  psi <- mean(QnstarAd1 - QnstarAd0)
  ic <- hAA * (Y - QnstarAA) + QnstarAd1 - QnstarAd0 - psi
  var <- var(ic)/length(ic)
  stderr <- sqrt(var)
  lcl <- psi - 1.96 * stderr
  ucl <- psi + 1.96 * stderr
  pvalue <- min(1, 2*(1-pnorm(psi/stderr)))
  
  data.frame(psi=psi, var=var, lcl = lcl, ucl=ucl, pvalue=pvalue)
}