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

## ----setup--------------------------------------------------------------------
library(TrialEmulation)

## -----------------------------------------------------------------------------
trial_pp <- trial_sequence(estimand = "PP") # Per-protocol
trial_itt <- trial_sequence(estimand = "ITT") # Intention-to-treat

## -----------------------------------------------------------------------------
trial_pp_dir <- file.path(tempdir(), "trial_pp")
dir.create(trial_pp_dir)
trial_itt_dir <- file.path(tempdir(), "trial_itt")
dir.create(trial_itt_dir)

## -----------------------------------------------------------------------------
data("data_censored")
trial_pp <- trial_pp |>
  set_data(
    data = data_censored,
    id = "id",
    period = "period",
    treatment = "treatment",
    outcome = "outcome",
    eligible = "eligible"
  )

# Function style without pipes
trial_itt <- set_data(
  trial_itt,
  data = data_censored,
  id = "id",
  period = "period",
  treatment = "treatment",
  outcome = "outcome",
  eligible = "eligible"
)

## -----------------------------------------------------------------------------
trial_itt

## -----------------------------------------------------------------------------
ipw_data(trial_itt)

## -----------------------------------------------------------------------------
trial_pp <- trial_pp |>
  set_switch_weight_model(
    numerator = ~age,
    denominator = ~ age + x1 + x3,
    model_fitter = stats_glm_logit(save_path = file.path(trial_pp_dir, "switch_models"))
  )
trial_pp@switch_weights

## -----------------------------------------------------------------------------
trial_pp <- trial_pp |>
  set_censor_weight_model(
    censor_event = "censored",
    numerator = ~x2,
    denominator = ~ x2 + x1,
    pool_models = "none",
    model_fitter = stats_glm_logit(save_path = file.path(trial_pp_dir, "switch_models"))
  )
trial_pp@censor_weights

## -----------------------------------------------------------------------------
trial_itt <- set_censor_weight_model(
  trial_itt,
  censor_event = "censored",
  numerator = ~x2,
  denominator = ~ x2 + x1,
  pool_models = "numerator",
  model_fitter = stats_glm_logit(save_path = file.path(trial_itt_dir, "switch_models"))
)
trial_itt@censor_weights

## -----------------------------------------------------------------------------
trial_pp <- trial_pp |> calculate_weights()
trial_itt <- calculate_weights(trial_itt)

## -----------------------------------------------------------------------------
show_weight_models(trial_itt)

## -----------------------------------------------------------------------------
trial_pp <- set_outcome_model(trial_pp)
trial_itt <- set_outcome_model(trial_itt, adjustment_terms = ~x2)

## -----------------------------------------------------------------------------
trial_pp <- set_expansion_options(
  trial_pp,
  output = save_to_datatable(),
  chunk_size = 500
)
trial_itt <- set_expansion_options(
  trial_itt,
  output = save_to_datatable(),
  chunk_size = 500
)

## ----eval = FALSE-------------------------------------------------------------
# trial_pp <- trial_pp |>
#   set_expansion_options(
#     output = save_to_csv(file.path(trial_pp_dir, "trial_csvs")),
#     chunk_size = 500
#   )
# 
# trial_itt <- set_expansion_options(
#   trial_itt,
#   output = save_to_csv(file.path(trial_itt_dir, "trial_csvs")),
#   chunk_size = 500
# )
# 
# trial_itt <- set_expansion_options(
#   trial_itt,
#   output = save_to_duckdb(file.path(trial_itt_dir, "trial_duckdb")),
#   chunk_size = 500
# )
# 
# # For the purposes of this vignette the previous `save_to_datatable()` output
# # option is used in the following code.

## -----------------------------------------------------------------------------
trial_pp <- expand_trials(trial_pp)
trial_itt <- expand_trials(trial_itt)

## -----------------------------------------------------------------------------
trial_pp@expansion

## -----------------------------------------------------------------------------
trial_itt <- load_expanded_data(trial_itt, seed = 1234, p_control = 0.5)

## -----------------------------------------------------------------------------
x2_sq <- outcome_data(trial_itt)$x2^2
outcome_data(trial_itt)$x2_sq <- x2_sq
head(outcome_data(trial_itt))

## -----------------------------------------------------------------------------
trial_itt <- fit_msm(
  trial_itt,
  weight_cols = c("weight", "sample_weight"),
  modify_weights = function(w) {
    q99 <- quantile(w, probs = 0.99)
    pmin(w, q99)
  }
)

## -----------------------------------------------------------------------------
trial_itt@outcome_model

## -----------------------------------------------------------------------------
trial_itt@outcome_model@fitted@model$model
trial_itt@outcome_model@fitted@model$vcov

## -----------------------------------------------------------------------------
trial_itt

## -----------------------------------------------------------------------------
preds <- predict(
  trial_itt,
  newdata = outcome_data(trial_itt)[trial_period == 1, ],
  predict_times = 0:10,
  type = "survival",
)
plot(preds$difference$followup_time, preds$difference$survival_diff,
  type = "l", xlab = "Follow up", ylab = "Survival difference"
)
lines(preds$difference$followup_time, preds$difference$`2.5%`, type = "l", col = "red", lty = 2)
lines(preds$difference$followup_time, preds$difference$`97.5%`, type = "l", col = "red", lty = 2)