## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

## ----setup, message = FALSE---------------------------------------------------
library(postcard)
withr::local_seed(1395878)
withr::local_options(list(postcard.verbose = 0))

## -----------------------------------------------------------------------------
n_train <- 2000
n_test <- 200
b1 <- 0.9
b2 <- 0.2
b3 <- -0.4
b4 <- -0.6

train_pois <- glm_data(
    Y ~ b1*log(X1)+b2*X2+b3*X3+b4*X2*X3,
    X1 = runif(n_train, min = 1, max = 10),
    X2 = rnorm(n_train),
    X3 = rgamma(n_train, shape = 1),
    family = poisson()
  )
test_pois <- glm_data(
    Y ~ b1*log(X1)+b2*X2+b3*X3+b4*X2*X3,
    X1 = runif(n_test, min = 1, max = 10),
    X2 = rnorm(n_test),
    X3 = rgamma(n_test, shape = 1),
    family = poisson()
  )

## -----------------------------------------------------------------------------
lrnr <- fit_best_learner(list(mod = Y ~ X1 + X2 + X3), data = train_pois)
preds <- dplyr::pull(predict(lrnr, new_data = test_pois))

power_marginaleffect(
  response = test_pois$Y,
  predictions = preds,
  var1 = function(var0) 1.1 * var0,
  kappa1_squared = 2,
  estimand_fun = "rate_ratio",
  target_effect = 1.4,
  exposure_prob = 1/2,
  margin = 0.8
)

## -----------------------------------------------------------------------------
ancova_prog_list <- list(
  ANCOVA = glm(Y ~ X1 + X2 + X3, data = train_pois, family = poisson),
  # "Null model" = glm(Y ~ 1, data = train_pois, family = poisson),
  "ANCOVA with prognostic score" = fit_best_learner(list(mod = Y ~ X1 + X2 + X3), data = train_pois)
)
test_pois_fun <- function(n) {
  glm_data(
    Y ~ b1*log(X1)+b2*X2+b3*X3+b4*X2*X3,
    X1 = runif(n, min = 1, max = 10),
    X2 = rnorm(n),
    X3 = rgamma(n, shape = 1),
    family = poisson()
  )
}

## -----------------------------------------------------------------------------
rpm <- repeat_power_marginaleffect(
  model_list = ancova_prog_list,
  test_data_fun = test_pois_fun,
  ns = seq(20, 500, by = 10), n_iter = 20,
  var1 = function(var0) 1.1 * var0,
  kappa1_squared = function(kap0) 1.1 * kap0,
  estimand_fun = "rate_ratio",
  target_effect = 1.4,
  exposure_prob = 1/2,
  margin = 0.8
)

plot(rpm)