Custom Intervention Contrasts

The current gbkmr_run() interface supports custom intervention contrasts directly. You do not need to call any function from R/old/, and you usually do not need to call the lower-level run_gbkmr_panel().

Default Contrast

By default, gbkmr_run() compares all mixture columns at the 25th percentile with all mixture columns at the 75th percentile:

fit_default <- gbkmr_run(
  data = prepared,
  time_points = 3,
  a_probs = c(0.25, 0.75),
  iter = 15000,
  K = 1000
)

This means:

low scenario:  logM1_0, logM2_0, ..., logM2_2 at their 25th percentiles
high scenario: logM1_0, logM2_0, ..., logM2_2 at their 75th percentiles

Change Quantiles

Use a_probs for a different global quantile shift:

fit_10_90 <- gbkmr_run(
  data = prepared,
  time_points = 3,
  a_probs = c(0.10, 0.90),
  iter = 15000,
  K = 1000
)

The two probabilities must be length two, with the low scenario first and the high scenario second.

Use Fixed Exposure Values

Use a_vals and astar_vals when the scientific question is about specific exposure levels. These vectors must be named with the mixture columns created by prepare_gbkmr_data().

exposure_names <- grep("^logM\\d+_\\d+$", names(prepared), value = TRUE)

# prepare_gbkmr_data(log_transform_mixtures = TRUE) stores log concentrations.
a_background <- setNames(rep(log(1), length(exposure_names)), exposure_names)
a_high <- setNames(rep(log(5), length(exposure_names)), exposure_names)

fit_fixed <- gbkmr_run(
  data = prepared,
  time_points = 3,
  a_vals = a_background,
  astar_vals = a_high,
  iter = 15000,
  K = 1000
)

If log_transform_mixtures = FALSE, pass values on the original scale. If it is TRUE, pass values on the log scale because the fitted data contain log-transformed mixture columns.

Change One Component

The names let you define component-specific and time-specific scenarios. This example increases only component 2 at every visit:

a_component <- setNames(rep(log(1), length(exposure_names)), exposure_names)
astar_component <- a_component
astar_component[c("logM2_0", "logM2_1", "logM2_2")] <- log(5)

fit_component <- gbkmr_run(
  data = prepared,
  time_points = 3,
  a_vals = a_component,
  astar_vals = astar_component,
  iter = 15000,
  K = 1000
)

What Not To Include

Only mixture columns belong in a_vals and astar_vals. Do not include baseline covariates or time-varying confounders:

# Wrong: waist_0 is a confounder, not a mixture exposure.
a_wrong <- c(logM1_0 = 0, logM2_0 = 0, waist_0 = 80)

Changing confounders inside the intervention vector would change the estimand. In g-computation, time-varying confounders are sampled forward under each mixture intervention.

Compare Multiple Contrasts

contrast_grid <- list(
  q25_q75 = c(0.25, 0.75),
  q10_q90 = c(0.10, 0.90),
  q05_q95 = c(0.05, 0.95)
)

fits <- lapply(contrast_grid, function(prob) {
  gbkmr_run(
    data = prepared,
    time_points = 3,
    a_probs = prob,
    iter = 15000,
    K = 1000
  )
})

do.call(rbind, lapply(names(fits), function(name) {
  ce <- fits[[name]]$causal_effect
  data.frame(
    contrast = name,
    estimate = ce$estimate,
    lower = ce$lower,
    upper = ce$upper
  )
}))