## ----config, include=FALSE----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, cache = TRUE, collapse=TRUE)

## ----setupfuture--------------------------------------------------------------
future::plan(future::multicore)

## ----eval=FALSE---------------------------------------------------------------
# progressr::handlers(global = TRUE)
# progressr::handlers(progressr::handler_progress(times = 100))

## ----message=FALSE------------------------------------------------------------
library("carts")
set.seed(42)

## -----------------------------------------------------------------------------
covariate <- function(n, p.a = 0.5) {
  return(
    data.frame(
    a = rbinom(n, 1, p.a),
    w = rnorm(n)
  )
  )
}

## -----------------------------------------------------------------------------
outcome <- setargs(
  outcome_continuous,
  mean = ~ 1 + a + w,
  par = c(10, -0.5, -1.2), # coef order as defined by the above formula
  sd = 1.5 ** 0.5
)

## -----------------------------------------------------------------------------
trial <- Trial$new(
  covariates = covariate,
  outcome = outcome,
  exclusion = identity, # no exclusion criteria (default)
  info = "Two-armed parallel design" # optional info
)

## -----------------------------------------------------------------------------
dd <- trial$simulate(n = 1e4)

## -----------------------------------------------------------------------------
table(dd$a)

## -----------------------------------------------------------------------------
head(dd, 2)

## -----------------------------------------------------------------------------
trial$simulate(n = 1e4, p.a = 0.9)$a |> table()

## -----------------------------------------------------------------------------
trial0 <- Trial$new(
  covariates = covariate,
  outcome = function(data, p.a = 0.1) rbinom(nrow(data), 1, p.a)
)
trial0$simulate(1e4)[, list(y = mean(y), a = mean(a))]
trial0$simulate(1e4, p.a = 1)[, list(y = mean(y), a = mean(a))]

## -----------------------------------------------------------------------------
trial0 <- Trial$new(
  covariates = function(n) covariate(n),
  outcome = function(data, p.a = 0.1) rbinom(nrow(data), 1, p.a)
)
trial0$simulate(1e4, p.a = 1)[, list(y = mean(y), a = mean(a))]

## -----------------------------------------------------------------------------
trial <- Trial$new(
  covariates = covariate,
  outcome = outcome_continuous
)

## -----------------------------------------------------------------------------
trial$args_model(
  mean = ~ 1 + a + w,
  par = c(10, -0.5, -1.2)
)

## -----------------------------------------------------------------------------
trial$simulate(1e4)[, list(y = mean(y)), by = "a"]

## -----------------------------------------------------------------------------
trial$simulate(1e4, par = c(0, 20, 0))[, list(y = mean(y)), by = "a"]

## -----------------------------------------------------------------------------
estimator <- est_glm(
  response = "y", # default value
  treatment = "a" # default value
)
data <- trial$simulate(n = 300)
estimator(data)

## -----------------------------------------------------------------------------
trial$estimators(marginal = est_glm()) # using def. args. of est_glm

## -----------------------------------------------------------------------------
names(trial$estimators())

## -----------------------------------------------------------------------------
print(trial)

## -----------------------------------------------------------------------------
trial$estimate_power(
  n = 300, # sample size
  R = 500 # number of Monte Carlo replicates
)

## -----------------------------------------------------------------------------
trial$estimate_power(n = 300, R = 500,
  summary.args = list(alternative = ">") # one-sided test should
  # have no power as a < 0
)

## -----------------------------------------------------------------------------
trial$estimate_power(n = 300, R = 500, p.a = 0.1)

## -----------------------------------------------------------------------------
trial$estimators(adjusted = est_glm(covariates = "w"))
trial$estimate_power(n = 300, R = 500)

## -----------------------------------------------------------------------------
# trial$estimate_samplesize(
#   power = 0.9, # default
#   estimator = trial$estimators("adjusted")
# )