In Larsvanderlaan/npcausalML: Nonparametric estimation of the relative risk function

knitr::opts_chunk$set(echo = TRUE)

data_full <- sim(setD, n = 50000)
tsk1 <- sl3_Task$new(data_full, covariates = "W1", outcome = "gRtilde1")
lrnr <- Lrnr_gam$new()
lrnr <- lrnr$train(tsk1)
tsk1 <- sl3_Task$new(data, covariates = "W1", outcome = "gRtilde1")

tsk0 <- sl3_Task$new(data_full, covariates = "W1", outcome = "gRtilde0")
lrnr0 <- Lrnr_gam$new()
lrnr0 <- lrnr$train(tsk0)
tsk0 <- sl3_Task$new(data, covariates = "W1", outcome = "gRtilde1")

data.table(lrnr$predict(tsk1)/lrnr0$predict(tsk0), exp(preds))
 data.table(names(risks),risks)

risks_all <- list()

  # Generate data
  n <- 5000
#   V <- as.matrix(replicate(3,runif(n, min = -1, max = 1)))
#   X <- V
#   A <- rbinom(n, size = 1, prob = 0.12+ 0.75*plogis(sin(5*(V^2 %*% c(-1,-1,-1))) + sin(5*(V %*% c(1,-1,1))) + 0.5*cos(5*(V %*% c(-1,1,1)))^2 + sin(5*V) %*% c(1,1,-1) + cos(5*V) %*% c(1,1,1)))
#   Y <- rbinom(n, size = 1, prob = 0.05 + 0.9*plogis(0.9*(-0.5 + A + 0.3*A*(sin(3.5*V) %*% c(1,-1,1) + cos(3.5*V) %*% c(1,-1,1)) + sin(3*V) %*% c(-1,1,1) + cos(3*V) %*% c(1,1,-1))))
#  
# bound <- Vectorize(tmle3::bound)
D <- DAG.empty()
# hard_additive
theta = function(W1, W2, W3) {0.5*(W1+cos(4*W1)) + 0.5*sin(4*W2) + 0.5*sin(4*W3)}

theta <- Vectorize(theta)

key <- "hard_additive"
D <- D +
  node("W1", distr = "runif", min = -1, max = 1) +
  node("W2", distr = "runif", min = -1, max = 1) +
  node("W3", distr = "runif", min = -1, max = 1) +
  node("g", distr = "rconst", const = 0.15 + 0.75*plogis(cos((W1+W2 + W3)*5) * sin((W1+W2 + W3)*5) + cos((W1+W2 + W3)*5) + sin((W1+W2 + W3)*5) + sin(W1*5) + W1*sin(W1*5) + cos(W2*5) + 2*W1*W2 - sin(W3*5) + sin(5*W1*W3) + 2*W1*W2*W3 + W3*sin(W1*5) + cos(W2*4)*sin(W1*5) ) ) +
  node("A", distr = "rbinom", size = 1, prob = g )+
   node("phi", distr = "rconst", const = (1+W1 + W2)*sin(W1*5)  + (1+W2 + W3)*cos(5*W2) + 0.5*exp(W1*W2) + cos(6*W1*W3)+ sin(6*W3*W2) + sin(6*W1*W2) + W3*sin(5*W3) + W2*sin(5*W3)) +
  node("theta", distr = "rconst", const =  theta(W1,W2,W3))+
    node("gRtilde0", distr = "rconst",  const = (-(exp(theta) + 1)*exp(phi) + (exp(2*phi)*(exp(theta) + 1)^2 + 4*(exp(theta + phi))*(1 - exp(phi)))^(0.5))/(2*exp(theta)*(1-exp(phi)))
           ) +
   node("gRtilde1", distr = "rconst",  const =    gRtilde0*exp(theta)) +
   node("gRtilde", distr = "rconst",  const =   A*gRtilde1 + (1-A)*gRtilde0) +
   node("gR", distr = "rconst",  const =  gRtilde ) +
  node("R", distr = "rbinom", size = 1, prob = gR)+
  node("RR", distr = "rconst", const = gRtilde1/gRtilde0)

setD <- set.DAG(D, vecfun = c("bound", "round", "theta"))
data <- sim(setD, n = n)

  ##################
  ###################   WORKS
  #A <- rbinom(n, size = 1, prob = 0.2+ 0.6*plogis(sin(5*(V %*% c(1,1,1))) + 0.5*cos(5*(V %*% c(1,1,1)))^2 + sin(5*V) %*% c(1,1,1) + cos(5*V) %*% c(1,1,1)))
  #Y <- rbinom(n, size = 1, prob = 0.05 + 0.9*plogis(0.9*(-0.5 + A + 0.3*A*(sin(3.5*V) %*% c(1,-1,1) + cos(3.5*V) %*% c(1,-1,1)) + sin(3*V) %*% c(-1,1,1) + cos(3*V) %*% c(1,1,-1))))
    ##################
    ##################
  Y <- data$R
  A <- data$A
  X <- as.matrix(data[, c("W1", "W2", "W3")])
  V <- X[,1]

  Q1 <- data$gRtilde1
  Q0 <- data$gRtilde0
  Q <- data$gRtilde
  g1 <- data$g

  lrnr_A <- Lrnr_hal9001$new(max_degree = 1, smoothness_orders = 1, num_knots = c(10, 5))
  lrnr_Y <- Lrnr_hal9001$new(max_degree = 1, smoothness_orders = 1, num_knots = c(10, 5))
  lrnr_Y <- Stack$new(Lrnr_xgboost$new(max_depth = 7), Lrnr_xgboost$new(max_depth = 3), Lrnr_xgboost$new(max_depth = 6), Lrnr_xgboost$new(max_depth = 5), Lrnr_xgboost$new(max_depth = 4))
  lrnr_A <- lrnr_Y

  task <- make_task(V, X, A, Y, folds = 5)
  likelihood <- make_likelihood(task, lrnr_A ,lrnr_Y, lrnr_V = Lrnr_gam$new(family = binomial()), cv = T)
  genr <- make_generator(likelihood)
  task_RR <- genr(task, "validation")



  basis1 <- fourier_basis$new(orders = c(1,0,0), max_degrees = c(1,2,3))
  basis2 <- fourier_basis$new(orders = c(2,0,0), max_degrees = c(1,2,3))
  basis3 <- fourier_basis$new(orders = c(3,0,0), max_degrees = c(1,2,3))

  basis_list <- list("k=0" = NULL,"k=1"= basis1,
                     "k=2" = basis2, 
                     "k=3" =basis3
                    )

  lrnr <- Stack$new(Lrnr_weight_helper$new(lrnr = Lrnr_LRR_xgboost$new(max_depth = 4, nrounds = 20), name = "xgboost_4"),
                    Lrnr_weight_helper$new(lrnr =Lrnr_LRR_xgboost$new(max_depth = 3, nrounds = 20), name = "xgboost_3"),
            make_learner(Pipeline,lrnr_SL.gam1, Lrnr_chainer_link$new()),
   Lrnr_weight_helper$new(lrnr = Lrnr_LRR_hal9001$new(max_degree = 1, smoothness_orders =1, num_knots = c(15)), name = "hal9001_1"),
                       make_learner(Pipeline,lrnr_SL.gam2, Lrnr_chainer_link$new()),
                                                   make_learner(Pipeline,lrnr_SL.glm, Lrnr_chainer_link$new()),
   make_learner(Pipeline,lrnr_SL.inter, Lrnr_chainer_link$new())

  )
  lrnrs <- list()
  for(name in names(basis_list)) {
    basis <- basis_list[[name]]

    lrr_plugin <- LRR_plugin_task_generator$new(sieve_basis = basis, name = name)

    lrnrs <- c(lrnrs, list(Pipeline$new(lrr_plugin, lrnr$clone())))


  }
  lrnr_lrr <- make_learner(Stack, lrnrs)

  lrnr_lrr <- lrnr_lrr$train(task_RR)
  preds <- lrnr_lrr$predict(task_RR)


  risks <- apply(preds, 2 , function(preds) {
  f <- preds
  gRtilde1 <- Q1
  gRtilde0 <- Q0
  gRtilde0 <- bound(gRtilde0, 0.01)
  gRtilde1 <- bound(gRtilde1, 0.01)
  f0 = log(lrnr$predict(tsk1)/lrnr0$predict(tsk0))
  loss <- (f0-f)^2
  return(mean(loss))
})

Larsvanderlaan/npcausalML documentation built on July 30, 2023, 4:32 p.m.