knitr::opts_chunk$set(echo = TRUE)
library(simcausal)
library(npOddsRatio)
passes <- c()
for(i in 1:200){
  #n <- 100000
    # lrnr <- Lrnr_hal9001_custom$new(max_degree = 3, num_knots = c(30,20,1), fit_control = list(parallel = TRUE))
  # task <- sl3_Task$new(data, covariates =c("W1", "W2") ,outcome = "Y") 
  # lrnr_Y0 <- lrnr$train(task[data$A==0])
  # Q0 <- lrnr_Y0$predict(task)
  # true <- coef(glm.fit(cbind(rep(1,n), data$W1), data$Y, offset = qlogis(Q0), weights = data$A, family  = binomial()))
  # 

  n <- 5000
  D <- DAG.empty()

  D <- D +
    node("W1", distr = "runif", min = -1, max = 1) +
    node("W2", distr = "runif", min = -1, max = 1) +
    node("A", distr = "rbinom", size = 1, prob = plogis(W1 + W2) )+
    node("Z1", distr = "rgamma", shape=3, rate = 2) +
    node("Z", distr = "rconst", const=min(Z1,4)) +
    node("Q", distr = "rconst", const = plogis(-1+A + A*(-1.5*W1  -2*(Z-1.5)*(Z<=1.5) + 2*W1*Z   ) + 0.5*W1+ (Z-1.3) + (Z-1.5)*(Z>=1.5) *(1.5+W1) + (Z-2.5)*(Z>=2.5)+ 0.5*W2))+

    node("Y", distr = "rbinom", size = 1, prob = Q) +
    node("G", distr = "rconst",const = 1- 0.6*plogis(0.4*(W1 + W2  +A-0.5) + (Z-1.3)*(2.5+W1) + 2*(Z-2)*(Z>=2) ))+
    node("Delta", distr = "rbinom", size = 1, prob = G)

  setD <- set.DAG(D )
  data <- sim(setD, n = n)
  true <- c(1.3511903, 0.9620972)
  library(doMC)
  registerDoMC()
 lrnr <- Lrnr_hal9001_custom$new(max_degree = 3, num_knots = c(20,10,1), fit_control = list(parallel = TRUE))

  fit <- npORMissing(~ poly(W1,degree=1, raw = T)   , W = data[,c("W1", "W2")], A = data$A, Y = data$Y , Z = data$Z, sl3_learner_A = Lrnr_gam$new(),
                     Delta = data$Delta, sl3_learner_default = lrnr )
  ci <- fit$coefs

  passes <- cbind(passes, ci[,3] <= true & ci[,4] >= true)
 print(rowMeans(passes))

  print(apply(passes,1,table))
}


Larsvanderlaan/npOddsRatio documentation built on May 3, 2022, 12:05 p.m.