## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  fig.width = 7,
  fig.height = 4
)

## ----load, message=FALSE, warning=FALSE---------------------------------------
library(causalMed)
library(data.table)

## ----data---------------------------------------------------------------------
data("survivaldata")
dat <- as.data.table(survivaldata)
head(dat, 8)

## ----models-------------------------------------------------------------------
init_s <- recodes(lag1_A = 0, lag1_L = 0, lag1_M = 0)
in_s   <- recodes(lag1_A = A, lag1_L = L, lag1_M = M)

models_surv <- list(
  spec_model(A ~ V + lag1_A + lag1_L + time,
             var_type = "binary", mod_type = "exposure"),
  spec_model(L ~ V + A + lag1_L + time,
             var_type = "normal", mod_type = "covariate"),
  spec_model(M ~ V + A + L + lag1_M + time,
             var_type = "binary", mod_type = "mediator"),
  spec_model(C ~ V + L + time,
             var_type = "binary", mod_type = "censor"),    # loss to follow-up
  spec_model(Y ~ V + A + M + L + A:M + time,
             var_type = "binary", mod_type = "survival")   # discrete-time hazard
)

## ----fit-single---------------------------------------------------------------
fit_surv <- mediation(
  data           = dat,
  id_var         = "id",
  time_var       = "time",
  base_vars      = "V",
  exposure       = "A",
  outcome        = "Y",
  models         = models_surv,
  init_recode    = init_s,
  in_recode      = in_s,
  mediation_type = "I",
  mc_sample      = 20000,
  R              = 1,      # set R > 1 for bootstrap CIs (see below)
  quiet          = TRUE,
  seed           = 2026
)

## ----interventions-single-----------------------------------------------------
fit_surv$effect_size

## ----decomp-single------------------------------------------------------------
fit_surv$estimate

## ----nvw, eval=FALSE----------------------------------------------------------
# mediation(..., n_vw = 1)

## ----regime-------------------------------------------------------------------
fit_window <- mediation(
  data             = dat,
  id_var           = "id",
  time_var         = "time",
  base_vars        = "V",
  exposure         = "A",
  outcome          = "Y",
  models           = models_surv,
  init_recode      = init_s,
  in_recode        = in_s,
  mediation_type   = "I",
  exposure_regime  = c(0, 1, 1, 0, 0),   # exposed at t = 1 and 2 only
  reference_regime = 0,                  # never exposed
  mc_sample        = 20000,
  R                = 1,
  quiet            = TRUE,
  seed             = 2026
)
fit_window$estimate

## ----regime-support-----------------------------------------------------------
fit_window$data_summary$regime_support

## ----yam-data-----------------------------------------------------------------
data("yamamurodata")
yam <- as.data.table(yamamurodata)
head(yam, 6)

## ----yam-models---------------------------------------------------------------
models_yam <- list(
  spec_model(A ~ V + lag1_A + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 +
               I(as.integer(time == 2)),
             var_type = "binary", mod_type = "exposure", subset = time > 0),
  spec_model(L ~ V + A + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 +
               I(as.integer(time == 2)),
             var_type = "normal", mod_type = "covariate", subset = time > 0),
  spec_model(M1 ~ V + A + L + I(L^2) + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 +
               I(as.integer(time == 2)),
             var_type = "normal", mod_type = "mediator", subset = time > 0),
  spec_model(M2 ~ V + A + L + I(L^2) + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 +
               I(as.integer(time == 2)),
             var_type = "normal", mod_type = "mediator", subset = time > 0),
  spec_model(Y ~ V + A + L + I(L^2) + M1 + M2 + M1:M2 +
               I(as.integer(time == 1)) + I(as.integer(time == 2)),
             var_type = "binary", mod_type = "survival")
)

fit_yam <- mediation(
  data           = yam,
  id_var         = "id",
  time_var       = "time",
  base_vars      = c("V", "L0base", "M10base", "M20base"),
  exposure       = "A",
  outcome        = "Y",
  models         = models_yam,
  init_recode    = recodes(L = L0base, M1 = M10base, M2 = M20base,
                           lag1_A = 0, lag1_L = 0, lag1_M1 = 0, lag1_M2 = 0),
  in_recode      = recodes(lag1_A = A, lag1_L = L, lag1_M1 = M1, lag1_M2 = M2),
  mediation_type = "I",
  mc_sample      = 20000,
  R              = 1,
  quiet          = TRUE,
  seed           = 2026
)

fit_yam$effect_size

## ----yam-truth----------------------------------------------------------------
truth <- data.table(
  Effect = c("Total effect", "Direct effect", "Indirect effect (M1)",
             "Indirect effect (M2)", "TE - (Direct + Indirect)"),
  True   = c(-6.36, -3.20, -2.29, -0.97, 0.10)   # from ?yamamurodata
)
est <- as.data.table(fit_yam$estimate)[, .(Effect, Estimate = RD * 100)]
cmp <- merge(truth, est, by = "Effect", sort = FALSE)
cmp[, Difference := Estimate - True]
cmp

## ----censor-data--------------------------------------------------------------
set.seed(11)
dat_c <- copy(dat)
dat_c[, C := rbinom(.N, 1, plogis(-3.5 + 0.8 * L))]
dat_c[, after_cens := cumsum(shift(C, fill = 0)), by = id]
dat_c <- dat_c[after_cens == 0][, after_cens := NULL]
dat_c[C == 1, Y := 0L]   # censored before the event in that period

## ----fit-censor---------------------------------------------------------------
models_cens <- list(
  spec_model(A ~ V + lag1_A + lag1_L + time,
             var_type = "binary", mod_type = "exposure"),
  spec_model(L ~ V + A + lag1_L + time,
             var_type = "normal", mod_type = "covariate"),
  spec_model(M ~ V + A + L + lag1_M + time,
             var_type = "binary", mod_type = "mediator"),
  spec_model(C ~ V + A + L + time,
             var_type = "binary", mod_type = "censor"),     # censoring model
  spec_model(Y ~ V + A + M + L + A:M + time,
             var_type = "binary", mod_type = "survival")
)

fit_cens <- mediation(
  data           = dat_c,
  id_var         = "id",
  time_var       = "time",
  base_vars      = "V",
  exposure       = "A",
  outcome        = "Y",
  models         = models_cens,
  init_recode    = init_s,
  in_recode      = in_s,
  mediation_type = "I",
  mc_sample      = 20000,
  R              = 1,
  quiet          = TRUE,
  seed           = 2026
)

fit_cens$estimate

## ----subset, eval=FALSE-------------------------------------------------------
# # the exposure can occur only while gvhdm1 == 0 ...
# spec_model(gvhd ~ all + cmv + male + age + ...,
#            var_type = "binary", mod_type = "exposure",
#            subset = gvhdm1 == 0)
# 
# # ... and stays at 1 afterwards
# out_recode = recodes(gvhd = ifelse(gvhdm1 == 1, 1, gvhd))

## ----boot, eval=FALSE---------------------------------------------------------
# future::plan(future::multisession)
# 
# fit_ci <- mediation(
#   data           = dat,
#   id_var         = "id",
#   time_var       = "time",
#   base_vars      = "V",
#   exposure       = "A",
#   outcome        = "Y",
#   models         = models_surv,
#   init_recode    = init_s,
#   in_recode      = in_s,
#   mediation_type = "I",
#   mc_sample      = 20000,
#   R              = 500,
#   seed           = 2026
# )
# 
# future::plan(future::sequential)
# 
# fit_ci$estimate                  # with Sd, percentile and normal limits
# fit_ci$boot_estimates$effects    # the per-replicate estimates

## ----natural, warning=TRUE----------------------------------------------------
fit_nat <- mediation(
  data           = dat,
  id_var         = "id",
  time_var       = "time",
  base_vars      = "V",
  exposure       = "A",
  outcome        = "Y",
  models         = models_surv,
  init_recode    = init_s,
  in_recode      = in_s,
  mediation_type = "N",
  mc_sample      = 20000,
  R              = 1,
  quiet          = TRUE,
  seed           = 2026
)

## ----tmle-survival------------------------------------------------------------
fit_tmle_s <- mediation(
  data           = dat_c,
  id_var         = "id",
  time_var       = "time",
  base_vars      = "V",
  exposure       = "A",
  outcome        = "Y",
  models         = models_cens,
  init_recode    = init_s,
  in_recode      = in_s,
  mediation_type = "N",
  estimator      = "tmle",
  quiet          = TRUE,
  seed           = 2026
)

fit_tmle_s$estimate