testWorkingModel.Rmd

---
title: "Test for RR working models"
output: html_document
date: '2022-03-09'
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}

coverage <- c()
for(i in 1:200){
  n <- 1500
  W1 <- runif(n, -1 , 1)
  W2 <- runif(n, -1 , 1)# rbinom(n, 1 , plogis(W1))
  W <- cbind(W1, W2)
  A <- rbinom(n, 1 , plogis(0.5*(W1 + W2 )))
  Y <- rpois(n, plogis( W1 + W2 ) * exp(A*(1+ 0.5*(W1 + W2))))
  quantile(Y)
  LRR <- (1 + W1 + W2)
  true_coefs <- c(1, 0.5, 0.5)
  R <- as.numeric(1:n %in% c(which(Y>0), which(Y==0)[rbinom(sum(Y==0), size = 1, prob = 0.5)==1]))
  R <- rep(1,n)
  pR0 <- mean(R[Y==0])
  pR1 <- mean(R[Y>0])
  weights <- R / ifelse(Y>0, pR1, pR0)
  
  
  sl3_Learner_EYAW <- Lrnr_hal9001$new(family = "poisson", smoothness_orders =1, num_knots = c(10,5), max_degree = 2)
  sl3_Learner_pA1W <- Lrnr_gam$new()
  fit_msm <- npRRWorkingModel(~ 1 + W1 + W2, W, A, Y, weights, EY1W, EY0W, pA1W, sl3_Learner_EYAW, sl3_Learner_pA1W, outcome_type = "continuous") 
  
  coverage <- cbind(coverage, 
fit_msm$coefficients$ci_right >= true_coefs & fit_msm$coefficients$ci_left <= true_coefs)
}
rowMeans(coverage)
```