## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  message   = FALSE,
  warning   = FALSE,
  fig.width = 7,
  fig.height = 4
)
has_gfoRmula <- requireNamespace("gfoRmula", quietly = TRUE)

## ----pkgs---------------------------------------------------------------------
library(causalMed)
library(data.table)

## ----data-prep----------------------------------------------------------------
data("nonsurvivaldata", package = "causalMed")
dat <- as.data.table(nonsurvivaldata)
dat[, time := as.integer(time)]
setorder(dat, id, time)

# Lag recodes (applied inside the Monte Carlo loop)
init_rc <- recodes(lag1_A = 0, lag1_L1 = 0, lag1_L2 = 0)
in_rc   <- recodes(lag1_A = A, lag1_L1 = L1, lag1_L2 = L2)

## ----cm-binary----------------------------------------------------------------
models_cm <- list(
  spec_model(L1 ~ lag1_A + lag1_L1 + V + time,
             var_type = "normal", mod_type = "covariate"),
  spec_model(L2 ~ lag1_A + lag1_L2 + V + time,
             var_type = "binary", mod_type = "covariate"),
  spec_model(A  ~ lag1_A + L1 + L2 + V + time,
             var_type = "binary", mod_type = "exposure"),
  spec_model(Y_bin ~ A + L1 + L2,
             var_type = "binary", mod_type = "outcome")
)

fit_cm <- causalMed::gformula(
  data         = dat,
  id_var       = "id",
  time_var     = "time",
  base_vars    = "V",
  exposure     = "A",
  models       = models_cm,
  intervention = list(natural = NULL, always = 1, never = 0),
  ref_int      = "natural",
  init_recode  = init_rc,
  in_recode    = in_rc,
  mc_sample    = 10000,
  R            = 1,
  quiet        = TRUE,
  seed         = 20250915
)

fit_cm$effect_size

## ----gf-binary, eval = has_gfoRmula-------------------------------------------
suppressPackageStartupMessages(requireNamespace("gfoRmula"))

fit_gf <- gfoRmula::gformula(
  obs_data     = copy(dat),
  id           = "id",
  time_name    = "time",
  time_points  = length(unique(dat$time)),
  covnames     = c("L1", "L2", "A"),
  covtypes     = c("normal", "binary", "binary"),
  basecovs     = "V",
  covparams    = list(covmodels = c(
    L1 ~ lag1_A + lag1_L1 + V + time,
    L2 ~ lag1_A + lag1_L2 + V + time,
    A  ~ lag1_A + L1 + L2 + V + time
  )),
  ymodel       = Y_bin ~ A + L1 + L2,
  outcome_name = "Y_bin",
  outcome_type = "binary_eof",
  histories    = c(gfoRmula::lagged),
  histvars     = list(c("A", "L1", "L2")),
  intervention1.A = list(gfoRmula::static, rep(1, length(unique(dat$time)))),
  intervention2.A = list(gfoRmula::static, rep(0, length(unique(dat$time)))),
  int_descript = c("Always treat", "Never treat"),
  nsimul       = 10000,
  seed         = 20250915
)

fit_gf$result[, c("k", "g-form mean")]

## ----compare-binary, eval = has_gfoRmula--------------------------------------
gf_res <- fit_gf$result[k == max(fit_gf$result$k),
                         c("Interv.", "g-form mean")]
gf_res[, Intervention := c("natural", "always", "never")]
setnames(gf_res, "g-form mean", "gfoRmula")
gf_res[, `:=`(gfoRmula = round(gfoRmula, 4), Interv. = NULL)]

cm_res <- fit_cm$effect_size[, .(Intervention, causalMed = round(Est, 4))]
merge(cm_res, gf_res, by = "Intervention")

## ----cm-survival--------------------------------------------------------------
data("survivaldata", package = "causalMed")
dat_s <- as.data.table(survivaldata)
setorder(dat_s, id, time)

models_surv <- list(
  spec_model(L ~ V + lag1_L + time,
             var_type = "normal", mod_type = "covariate"),
  spec_model(A ~ V + lag1_A + lag1_L + L + time,
             var_type = "binary", mod_type = "exposure"),
  spec_model(Y ~ lag1_A + A + L + lag1_L + time,
             var_type = "binary", mod_type = "survival")
)

fit_cm_s <- causalMed::gformula(
  data         = dat_s,
  id_var       = "id",
  time_var     = "time",
  base_vars    = "V",
  exposure     = "A",
  models       = models_surv,
  intervention = list(natural = NULL, never = 0, always = 1),
  ref_int      = "natural",
  init_recode  = recodes(lag1_L = 0, lag1_A = 0),
  in_recode    = recodes(lag1_L = L, lag1_A = A),
  mc_sample    = 10000,
  R            = 1,
  quiet        = TRUE,
  seed         = 20250915
)

fit_cm_s$effect_size

## ----gf-survival, eval = has_gfoRmula-----------------------------------------
T_max <- length(unique(dat_s$time))

fit_gf_s <- gfoRmula::gformula(
  obs_data     = copy(dat_s),
  id           = "id",
  time_name    = "time",
  time_points  = T_max,
  covnames     = c("L", "A"),              # same order as the causalMed list
  covtypes     = c("normal", "binary"),
  basecovs     = "V",
  covparams    = list(covmodels = c(
    L ~ V + lag1_L + time,
    A ~ V + lag1_A + lag1_L + L + time
  )),
  ymodel       = Y ~ lag1_A + A + L + lag1_L + time,
  outcome_name = "Y",
  outcome_type = "survival",
  histories    = c(gfoRmula::lagged),
  histvars     = list(c("A", "L")),
  intervention1.A = list(gfoRmula::static, rep(0, T_max)),
  intervention2.A = list(gfoRmula::static, rep(1, T_max)),
  int_descript = c("Never treat", "Always treat"),
  nsimul       = 10000,
  seed         = 20250915,
  nsamples     = 0
)

fit_gf_s$result[k == max(fit_gf_s$result$k), c("Interv.", "g-form risk")]

## ----compare-survival, eval = has_gfoRmula------------------------------------
gf_s_res <- fit_gf_s$result[k == max(fit_gf_s$result$k),
                              c("Interv.", "g-form risk")]
gf_s_res[, Intervention := c("natural", "never", "always")]
setnames(gf_s_res, "g-form risk", "gfoRmula")
gf_s_res[, `:=`(gfoRmula = round(gfoRmula, 4), Interv. = NULL)]

cm_s_res <- fit_cm_s$effect_size[, .(Intervention,
                                      causalMed = round(Est, 4))]
merge(cm_s_res, gf_s_res, by = "Intervention")

## ----dynamic-cm---------------------------------------------------------------
fit_dyn <- causalMed::gformula(
  data         = dat,
  id_var       = "id",
  time_var     = "time",
  base_vars    = "V",
  exposure     = "A",
  models       = models_cm,
  intervention = list(natural = NULL, treat_if_L1_pos = dyn_int(as.numeric(L1 > 0))),
  ref_int      = "natural",
  init_recode  = init_rc,
  in_recode    = in_rc,
  mc_sample    = 5000,
  R            = 1,
  quiet        = TRUE,
  seed         = 20250915
)

fit_dyn$effect_size
fit_dyn$estimate

## ----cross-val, eval = has_gfoRmula-------------------------------------------
fit_cv_cm <- causalMed::gformula(
  data         = dat, id_var = "id", time_var = "time", base_vars = "V",
  exposure     = "A", models = models_cm,
  intervention = list(natural = NULL, always = 1),
  ref_int      = "natural",
  init_recode  = init_rc, in_recode = in_rc,
  mc_sample    = 50000, R = 1, quiet = TRUE, seed = 20250915
)

fit_cv_gf <- gfoRmula::gformula(
  obs_data = copy(dat), id = "id", time_name = "time",
  time_points = length(unique(dat$time)),
  covnames = c("L1","L2","A"), covtypes = c("normal","binary","binary"),
  basecovs = "V",
  covparams = list(covmodels = c(
    L1 ~ lag1_A + lag1_L1 + V + time,
    L2 ~ lag1_A + lag1_L2 + V + time,
    A  ~ lag1_A + L1 + L2 + V + time
  )),
  ymodel = Y_bin ~ A + L1 + L2, outcome_name = "Y_bin",
  outcome_type = "binary_eof",
  histories = c(gfoRmula::lagged), histvars = list(c("A","L1","L2")),
  intervention1.A = list(gfoRmula::static, rep(1, length(unique(dat$time)))),
  int_descript = "Always treat",
  nsimul = 50000, seed = 20250915
)

cm_vals <- round(fit_cv_cm$effect_size$Est, 5)
gf_vals <- round(fit_cv_gf$result[k == max(fit_cv_gf$result$k)][["g-form mean"]], 5)

data.frame(
  Intervention = c("natural", "always"),
  causalMed    = cm_vals,
  gfoRmula     = gf_vals,
  Difference   = round(abs(cm_vals - gf_vals), 5)
)

## ----versions-----------------------------------------------------------------
c(causalMed = as.character(packageVersion("causalMed")),
  gfoRmula  = if (has_gfoRmula) as.character(packageVersion("gfoRmula")) else NA)