Package 'forestsearch'

Description

Covariate set handed to forestsearch()/mrct_region_sims() for the paper's observed vs unobserved analysis cases.

Usage

analysis_covariates(
  dgm,
  case = c("observed", "unobserved"),
  drop = character(0)
)

Arguments

Value

Character vector of z_-prefixed covariate column names.

Description

Analyzes a subgroup and returns formatted results for summary table. Uses optimized cox_summary() and reduces redundant calculations.

Usage

analyze_subgroup(
  df_sub,
  outcome.name,
  event.name,
  treat.name,
  strata.name,
  subgroup_name,
  hr_a,
  potentialOutcome.name,
  return_medians,
  N
)

Arguments

Value

Character vector of results.

Examples

## Not run: 
library(survival)
df <- data.frame(
  tte   = gbsg$rfstime / 30.4375,
  event = gbsg$status,
  treat = gbsg$hormon
)
analyze_subgroup(
  df_sub        = df,
  outcome.name  = "tte",
  event.name    = "event",
  treat.name    = "treat",
  strata.name   = NULL,
  subgroup_name = "All",
  hr_a          = NA,
  potentialOutcome.name = NULL,
  return_medians = TRUE,
  N             = nrow(df)
)

## End(Not run)

Description

Describes the size-matched random benchmark subgroups regenerated in every simulated trial. With the default nested = TRUE, a single index draw of max(sizes) patients is taken per trial and the smaller benchmarks are its prefixes (random15 $\subset$ random20 $\subset$ random40 $\subset$ random60), exactly as in the extreme-subgroups vignettes. Membership columns are named paste0(prefix, sizes) and take values 0/1.

Usage

benchmark_spec(sizes = c(60L, 40L, 20L, 15L), nested = TRUE, prefix = "random")

Arguments

Value

An object of class "benchmark_spec".

See Also

run_subgroup_sims()

Examples

benchmark_spec()
benchmark_spec(sizes = c(100, 50), prefix = "rnd")

Description

Runs bootstrap analysis for ForestSearch, fitting Cox models and computing bias-corrected estimates and valid CIs (see vignette for references)

Usage

bootstrap_results(
  fs.est,
  df_boot_analysis,
  cox.formula.boot,
  nb_boots,
  show_three,
  H_obs,
  Hc_obs,
  seed = 8316951L,
  estimator_fn = NULL,
  effect_measure = NULL,
  boot_index_mat = NULL,
  mr_in_replicates = FALSE
)

Arguments

Value

Data.table with one row per bootstrap iteration and columns:

boot_id

Integer. Bootstrap iteration number (1 to nb_boots)

H_biasadj_1

Bias-corrected estimate for H using method 1: H_obs - (Hstar_star - Hstar_obs)

H_biasadj_2

Bias-corrected estimate for H using method 2: 2*H_obs - (H_star + Hstar_star - Hstar_obs)

Hc_biasadj_1

Bias-corrected estimate for H^c using method 1

Hc_biasadj_2

Bias-corrected estimate for H^c using method 2

max_sg_est

Numeric. Maximum subgroup hazard ratio found

L

Integer. Number of candidate factors evaluated

max_count

Integer. Maximum number of factor combinations

events_H_0

Integer. Number of events in control arm of original subgroup H on bootstrap sample

events_H_1

Integer. Number of events in treatment arm of original subgroup H on bootstrap sample

events_Hc_0

Integer. Number of events in control arm of original subgroup H^c on bootstrap sample

events_Hc_1

Integer. Number of events in treatment arm of original subgroup H^c on bootstrap sample

events_Hstar_0

Integer. Number of events in control arm of new subgroup H* on original data

events_Hstar_1

Integer. Number of events in treatment arm of new subgroup H* on original data

events_Hcstar_0

Integer. Number of events in control arm of new subgroup H^c* on original data

events_Hcstar_1

Integer. Number of events in treatment arm of new subgroup H^c* on original data

tmins_search

Numeric. Minutes spent on subgroup search in this iteration

tmins_iteration

Numeric. Total minutes for this bootstrap iteration

Pcons

Numeric. Consistency p-value for top subgroup

any_found

Integer 0/1. Identification flag; equivalent to as.integer(!is.na(Pcons)). Added in v0.2.0 to mirror the fold_summary$any_found column produced by forestsearch_tenfold, so cross-API diagnostics (sum(results$any_found == 1L)) work identically for bootstrap and CV outputs.

hr_sg

Numeric. Hazard ratio for top subgroup

N_sg

Integer. Sample size of top subgroup

E_sg

Integer. Number of events in top subgroup

K_sg

Integer. Number of factors defining top subgroup

g_sg

Numeric. Subgroup group ID

m_sg

Numeric. Subgroup index

M.1

Character. First factor label

M.2

Character. Second factor label

M.3

Character. Third factor label

M.4

Character. Fourth factor label

M.5

Character. Fifth factor label

M.6

Character. Sixth factor label

M.7

Character. Seventh factor label

grf_cuts_b

Character. GRF policy-tree cut expressions returned on this bootstrap sample, collapsed with " | " when multiple cuts are produced. NA_character_ when GRF was not used, the bootstrap forestsearch call errored, or GRF returned no cuts. Populated independently of whether a subgroup was identified: GRF may surface a candidate cut on bootstraps where the downstream consistency stage rejects all candidates. Mirrors the grf_cuts column that forestsearch_tenfold captures in its fold_summary return slot.

Rows where no valid subgroup was found will have NA for bias corrections. The returned object has a "timing" attribute with summary statistics.

Bias Correction Methods

Two bias correction approaches are implemented:

1. Method 1 (Simple Optimism):

   $H_{adj1} = H_{obs} - (H^*_{*} - H^*_{obs})$

   where $H^*_{*}$ is the new subgroup HR on bootstrap data and $H^*_{obs}$ is the new subgroup HR on original data.

2. Method 2 (Double Bootstrap):

   $H_{adj2} = 2 \times H_{obs} - (H_{*} + H^*_{*} - H^*_{obs})$

   where $H_{*}$ is the original subgroup HR on bootstrap data.

where:

• H_obs: Original subgroup HR on original data

• H_star: Original subgroup HR on bootstrap data

• Hstar_obs: New subgroup (found in bootstrap) HR on original data

• Hstar_star: New subgroup (found in bootstrap) HR on bootstrap data

Computational Details

• Uses doFuture backend for parallel execution (configured externally)

• Sets reproducible seeds: 8316951 + boot * 100 for each iteration

• Each bootstrap iteration runs full ForestSearch pipeline including variable selection, subgroup search, and consistency evaluation

• Sequential execution within each bootstrap prevents nested parallelization

• Failed bootstrap iterations generate warnings but don't stop execution

• Confounders are removed from bootstrap data to force fresh variable selection

Bootstrap Configuration

Each bootstrap iteration modifies ForestSearch arguments to:

• Suppress output: details, show_candidate_summary, plot.sg, plot.grf all set to FALSE

• Force re-selection: grf_res and grf_cuts set to NULL

• Prevent nested parallel: parallel_args$plan = "sequential", workers = 1

Performance Considerations

• Typical runtime: 1-5 seconds per bootstrap iteration

• For 1000 bootstraps with 6 workers: ~3-10 minutes total

• Memory usage scales with dataset size and number of workers

• Consider reducing nb_boots for initial testing (e.g., 100)

Error Handling

The function gracefully handles three failure modes:

1. Bootstrap sample creation fails: Returns row with all NA

2. ForestSearch fails to run: Warns and returns row with all NA

3. ForestSearch runs but finds no subgroup: Returns row with all NA

All three cases ensure the foreach loop can still combine results via rbind.

Note

This function is designed to be called within a foreach loop with %dofuture% operator. It requires:

• All functions in get_bootstrap_exports to be available in the parallel workers

• Packages listed in BOOTSTRAP_REQUIRED_PACKAGES to be installed

• Proper parallel backend setup via setup_parallel_SGcons

See Also

forestsearch_bootstrap_dofuture for the wrapper function that sets up parallelization and calls this function build_cox_formula for creating the Cox formula fit_cox_models for initial Cox model fitting get_Cox_sg for Cox model fitting on subgroups get_dfRes for processing bootstrap results into confidence intervals bootstrap_ystar for generating the Ystar matrix

Examples

## Not run: 
# Typically called via forestsearch_bootstrap_dofuture()
# Manual usage for debugging:

# 1. Fit initial ForestSearch model
fs_result <- forestsearch(
  df.analysis = mydata,
  outcome.name = "time",
  event.name = "status",
  treat.name = "treatment",
  confounders.name = c("age", "sex", "stage")
)

# 2. Build Cox formula
cox_formula <- build_cox_formula("time", "status", "treatment")

# 3. Get observed estimates
cox_fits <- fit_cox_models(fs_result$df.est, cox_formula)

# 4. Set up parallel backend
library(doFuture)
plan(multisession, workers = 6)

# 5. Run bootstrap (note: this is already parallelized internally)
boot_results <- bootstrap_results(
  fs.est = fs_result,
  df_boot_analysis = fs_result$df.est,
  cox.formula.boot = cox_formula,
  nb_boots = 100,
  show_three = TRUE,
  H_obs = cox_fits$H_obs,
  Hc_obs = cox_fits$Hc_obs
)

# 6. Check results
summary(boot_results)

# Proportion of bootstraps that found a subgroup
mean(!is.na(boot_results$H_biasadj_2))

## End(Not run)

Description

Returns the $B \times N$ integer count matrix consumed by the infinitesimal-jackknife variance estimator in get_dfRes. Entry $(b, j)$ is the number of times row $j$ of df was selected on bootstrap iteration $b$.

Usage

bootstrap_ystar(df, nb_boots, seed = 8316951L)

Arguments

Value

Integer matrix with nb_boots rows and nrow(df) columns.

Implementation note (v0.2.0)

This function previously generated the count matrix via a parallel foreach pass, with counts derived by matching against a user-supplied id column. It now runs on the main process and counts by positional row index via tabulate, matching how bootstrap_results constructs bootstrap data (df[in_boot, ]). For data with unique subject ids (the standard clinical-trial case) the two implementations are equivalent. The new implementation is also robust to non-unique id values.

Reproducibility note

Bootstrap indices generated by this version are not bit-identical to those produced by v0.1.x for the same seed (the generator now uses the main-process Mersenne-Twister stream rather than per-task L'Ecuyer streams). Statistical properties of bootstrap summaries are preserved.

Examples

## Not run: 
df <- data.frame(id = 1:50, tte = rexp(50), event = rbinom(50, 1, 0.6),
                 treat = rep(0:1, 25))
ystar <- bootstrap_ystar(df, nb_boots = 10)
dim(ystar)         # 10 x 50
rowSums(ystar)     # all equal to 50

## End(Not run)

Description

Constructs a publication-quality gt table summarizing subgroup identification and classification rates across one or more data generation scenarios and analysis methods. The layout mirrors Table 4 of Leon et al. (2024) with metrics grouped by model scenario (null / alt) and columns for each analysis method.

Usage

build_classification_table(
  scenario_results,
  analyses = NULL,
  digits = 2,
  title = "Subgroup Identification and Classification Rates",
  n_sims = NULL,
  bold_threshold = 0.05,
  font_size = 12,
  subgroup_notation = c("harm", "benefit")
)

Arguments

Details

For each scenario the function computes:

• any(H): Proportion of simulations identifying any subgroup.

• sens(H): Mean sensitivity (only under alternative).

• sens(Hc): Mean specificity.

• ppv(H): Mean positive predictive value (only under alternative).

• ppv(Hc): Mean negative predictive value.

• avg|H|: Mean size of identified subgroup (when found).

Under the null hypothesis the rows are reduced to any(H), sens(Hc), ppv(Hc), and avg|H|.

Value

A gt table object.

See Also

format_oc_results, summarize_simulation_results

Examples

## Not run: 
# Assemble results from H0 and H1 simulations
scenarios <- list(
  null = list(
    results = results_null, label = "M1",
    n_sample = 700, dgm = dgm_null, hypothesis = "null"
  ),
  alt = list(
    results = results_alt, label = "M1",
    n_sample = 700, dgm = dgm_calibrated, hypothesis = "alt"
  )
)

build_classification_table(scenarios, n_sims = 100)

## End(Not run)

Description

Constructs a Cox model formula from variable names, optionally adjusted for additional covariates.

Usage

build_cox_formula(
  outcome.name,
  event.name,
  treat.name,
  adjust_covariates = NULL
)

Arguments

Value

An R formula object for Cox regression. When adjust_covariates is supplied, the treatment term is always placed first on the right-hand side.

Examples

build_cox_formula("time_months", "event", "treat")
build_cox_formula("time_months", "event", "treat",
                  adjust_covariates = c("strata(site)", "age"))

Description

Constructs a publication-quality gt table summarizing estimation properties for hazard ratios in the identified subgroup and its complement. The layout mirrors Table 5 of Leon et al. (2024), showing average estimate, empirical SD, min, max, and relative bias for each estimator.

Usage

build_estimation_table(
  results,
  dgm,
  analysis_method = "FSlg",
  n_boots = NULL,
  digits = 2,
  title = "Estimation Properties",
  subtitle = NULL,
  font_size = 12,
  cde_H = NULL,
  cde_Hc = NULL,
  subgroup_notation = c("harm", "benefit"),
  trim_threshold = 1000,
  trim_fraction = 0.01
)

Arguments

Details

Uses the paper's notation conventions:

• theta-dagger: Marginal (causal) HR truth

• theta-ddagger: Controlled direct effect (CDE) truth

• theta-hat(H-hat): Plugin Cox estimate in identified subgroup

• theta-hat*(H-hat): Bootstrap bias-corrected estimate

Includes both Cox-based HR and AHR (Average Hazard Ratio from loghr_po) estimators when AHR columns are present in the results.

For each subgroup (H and Hc) the function reports:

• Avg: Mean of the estimates across estimable simulations.

• SD: Empirical standard deviation.

• Min / Max: Range.

• b-dagger: Relative bias (percent) vs marginal truth, 100 * (Avg - theta_dagger) / theta_dagger.

• b-ddagger (conditional): Relative bias (percent) vs CDE truth, shown when CDE values are available.

When bootstrap-corrected columns (hr.H.bc, hr.Hc.bc) are present in results, an additional bias-corrected row (theta-hat*(H-hat)) is added per subgroup.

When AHR columns (ahr.H.hat, ahr.Hc.hat) are present, AHR estimation rows are appended using the DGM's true AHR values for relative bias calculation.

When CDE columns (cde.H.hat, cde.Hc.hat) are present and CDE truth values are available, CDE estimation rows (theta-ddagger(H-hat)) are appended. The b-dagger column for CDE rows reports bias relative to the CDE truth rather than the marginal HR.

Value

A gt table object, or NULL if no estimable realizations exist.

See Also

build_classification_table, format_oc_results, get_dgm_hr

Examples

## Not run: 
# Basic usage (marginal truth only)
build_estimation_table(
  results = results_alt,
  dgm = dgm_calibrated,
  analysis_method = "FSlg"
)

# With CDE truth (full Table 5 alignment)
build_estimation_table(
  results = results_alt,
  dgm = dgm_calibrated,
  cde_H = 2.25,
  cde_Hc = 0.60
)

## End(Not run)

Description

Uses root-finding to select a value of cens_adjust for simulate_from_dgm such that a chosen censoring summary statistic in the simulated data matches the corresponding statistic from the DGM reference data (dgm$df_super).

Usage

calibrate_cens_adjust(
  dgm,
  target = c("rate", "km_median"),
  n = 1000,
  rand_ratio = 1,
  analysis_time = 48,
  max_entry = 24,
  seed = 42,
  interval = c(-3, 3),
  tol = 1e-04,
  tol_rel = 2.5,
  n_eval = 2000,
  verbose = TRUE,
  ...
)

Arguments

Details

Two calibration targets are supported:

"rate"

Overall censoring rate (proportion censored). Finds cens_adjust such that mean(event_sim == 0) in simulated data equals mean(event == 0) in dgm$df_super.

"km_median"

KM-based median censoring time, estimated by reversing the event indicator so censored observations become the "event" of interest. Finds cens_adjust such that the simulated KM median matches the reference KM median.

How the objective function works

At each candidate cens_adjust value, the objective function:

1. Calls simulate_from_dgm() with n = n_eval and the candidate cens_adjust.

2. Calls check_censoring_dgm() with verbose = FALSE to extract the target metric.

3. Returns sim_metric - ref_metric.

uniroot finds the zero crossing, i.e. the cens_adjust at which simulated and reference metrics are equal.

Monotonicity

The objective is monotone in cens_adjust for both targets:

• Larger cens_adjust → longer censoring times → lower censoring rate and higher KM median.

• Smaller cens_adjust → shorter censoring times → higher censoring rate and lower KM median.

If uniroot fails (the target lies outside the search interval), the boundary values are printed and a wider interval should be tried.

Stochastic noise

Because the objective function involves simulation, there is Monte Carlo noise. Setting a fixed seed and a sufficiently large n_eval (>= 2000) reduces noise enough for reliable root-finding. The tol argument controls the root-finding tolerance on the cens_adjust scale (not the metric scale).

Value

A named list with elements:

cens_adjust

Calibrated cens_adjust value.

target

Calibration target used.

ref_value

Reference metric value from dgm$df_super.

sim_value

Achieved metric value in simulated data at the calibrated cens_adjust.

residual

Absolute difference between sim_value and ref_value.

iterations

Number of uniroot iterations.

diagnostic

Output of check_censoring_dgm at the calibrated value (invisibly).

See Also

simulate_from_dgm, check_censoring_dgm, generate_aft_dgm_flex

Examples

## Not run: 
library(survival)

# Build DGM on months scale
gbsg$time_months <- gbsg$rfstime / 30.4375

dgm <- generate_aft_dgm_flex(
  data            = gbsg,
  continuous_vars = c("age", "size", "nodes", "pgr", "er"),
  factor_vars     = c("meno", "grade"),
  outcome_var     = "time_months",
  event_var       = "status",
  treatment_var   = "hormon",
  subgroup_vars   = c("er", "meno"),
  subgroup_cuts   = list(er = 20, meno = 0)
)

# Calibrate so simulated censoring rate matches reference
cal_rate <- calibrate_cens_adjust(
  dgm           = dgm,
  target        = "rate",
  n             = 1000,
  analysis_time = 84,
  max_entry     = 24
)
cat("Calibrated cens_adjust (rate):", cal_rate$cens_adjust, "\n")

# Calibrate to KM median censoring time instead
cal_km <- calibrate_cens_adjust(
  dgm           = dgm,
  target        = "km_median",
  n             = 1000,
  analysis_time = 84,
  max_entry     = 24
)
cat("Calibrated cens_adjust (km_median):", cal_km$cens_adjust, "\n")

# Use calibrated value in simulation
sim <- simulate_from_dgm(
  dgm           = dgm,
  n             = 1000,
  analysis_time = 84,
  max_entry     = 24,
  cens_adjust   = cal_rate$cens_adjust,
  seed          = 123
)
mean(sim$event_sim)   # event rate
mean(sim$event_sim == 0)  # censoring rate — should match ref

## End(Not run)

Description

Searches over a grid of interaction multipliers (k_inter) to find the value that produces a target effect size in the subgroup $Q$. Returns a modified "glm_dgm" object with the calibrated interaction and updated true effects.

Usage

calibrate_glm_interaction(
  data,
  factor_vars,
  continuous_vars = NULL,
  outcome_var,
  treatment_var,
  target_effect,
  outcome_type = c("binary", "continuous", "count"),
  effect_measure = NULL,
  offset_var = NULL,
  subgroup_vars = NULL,
  subgroup_cuts = NULL,
  k_treat = 1,
  adverse_outcome = FALSE,
  k_inter_range = c(-10, 10),
  grid_step = 0.05,
  tol_rel = 2.5,
  n_super = 5000L,
  seed = 8316951L,
  verbose = FALSE
)

Arguments

Details

The function constructs a grid of k_inter values, calls generate_glm_dgm for each, and selects the value whose subgroup effect is closest to target_effect.

For binary outcomes with effect_measure = "OR", the target is on the switched-treatment OR scale (OR > 1 = treatment increases the outcome).

Value

An object of class "glm_dgm" with the calibrated k_inter and updated hazard_ratios.

See Also

generate_glm_dgm, simulate_from_glm_dgm

Examples

## Not run: 
dgm <- calibrate_glm_interaction(
  data          = actg_df,
  factor_vars   = paste0("z", 1:12),
  outcome_var   = "y_binary",
  treatment_var = "treat",
  target_effect = 2.0,
  outcome_type  = "binary",
  subgroup_vars = c("z1", "z2"),
  subgroup_cuts = list(z1 = 1L, z2 = 1L),
  verbose       = TRUE
)
print(dgm)

## End(Not run)

Description

Finds the interaction effect multiplier (k_inter) that achieves a target hazard ratio in the harm subgroup.

Usage

calibrate_k_inter(
  target_hr_harm,
  model = "alt",
  k_treat = 1,
  cens_type = "weibull",
  k_inter_range = c(-100, 100),
  tol = 1e-06,
  tol_rel = 2.5,
  use_ahr = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

This function uses uniroot to find the k_inter value such that the empirical HR (or AHR) in the harm subgroup equals target_hr_harm.

Value

Numeric value of k_inter that achieves the target HR

Examples

## Not run: 
# Find k_inter for HR = 1.5 in harm subgroup
k <- calibrate_k_inter(target_hr_harm = 1.5, verbose = TRUE)

# Verify
dgm <- create_gbsg_dgm(model = "alt", k_inter = k, verbose = TRUE)
print(dgm$hr_H_true)  # Should be close to 1.5

# Calibrate to AHR instead
k_ahr <- calibrate_k_inter(target_hr_harm = 1.5, use_ahr = TRUE, verbose = TRUE)
dgm_ahr <- create_gbsg_dgm(model = "alt", k_inter = k_ahr, verbose = TRUE)
print(dgm_ahr$AHR_H_true)  # Should be close to 1.5

## End(Not run)

Description

Finds the treatment effect multiplier (k_treat) that achieves a target overall hazard ratio in a DGM built by generate_aft_dgm_flex. The root is found by uniroot over k_treat with every other argument to generate_aft_dgm_flex() held fixed.

Usage

calibrate_k_treat(
  target_hr_overall,
  base_args,
  k_treat_range = c(-5, 5),
  tol = 1e-06,
  tol_rel = 2.5,
  use_ahr = FALSE,
  verbose = FALSE
)

Arguments

Details

This is the natural sibling of calibrate_k_inter (which targets the harm-subgroup HR under model = "alt") and calibrate_cens_adjust (which targets the censoring rate). It is typically used with model = "null" to fix the overall HR for a uniform-effect DGM, but also works with model = "alt" when the marginal HR is the calibration target rather than the subgroup-specific HR.

Reproducibility requires the uniroot objective to be deterministic in k_treat. This is achieved by including an integer seed field in base_args so that every call to generate_aft_dgm_flex() during the search uses the same random-number stream. If base_args$seed is NULL the objective becomes stochastic and uniroot may not converge; the function emits a warning in this case.

Value

Numeric scalar. The calibrated k_treat value. Stops with an error if the target cannot be reached – the search interval fails to bracket it even after automatic extension, or the achieved HR lies outside tol_rel percent of the target.

See Also

calibrate_k_inter, calibrate_cens_adjust, generate_aft_dgm_flex

Examples

## Not run: 
library(survival)
data(gbsg)
gbsg$time_months <- gbsg$rfstime / 30.4375

base_args <- list(
  data            = gbsg,
  continuous_vars = c("age", "size", "nodes", "pgr", "er"),
  factor_vars     = c("meno", "grade"),
  outcome_var     = "time_months",
  event_var       = "status",
  treatment_var   = "hormon",
  model           = "null",
  n_super         = 5000,
  seed            = 99,
  verbose         = FALSE
)

# Calibrate k_treat so that the overall HR is exactly 0.70
k <- calibrate_k_treat(target_hr_overall = 0.70,
                       base_args         = base_args,
                       verbose           = TRUE)

# Verify
dgm <- do.call(generate_aft_dgm_flex,
               c(base_args, list(k_treat = k)))
print(dgm$hazard_ratios$overall)   # ~ 0.70

# Calibrate to AHR instead
k_ahr <- calibrate_k_treat(target_hr_overall = 0.70,
                           base_args         = base_args,
                           use_ahr           = TRUE)
dgm_ahr <- do.call(generate_aft_dgm_flex,
                   c(base_args, list(k_treat = k_ahr)))
print(dgm_ahr$hazard_ratios$AHR)   # ~ 0.70

## End(Not run)

Description

Given per-subgroup approximation probabilities and simulated FPR values at multiple sample sizes, fits the power-law model $L_{\text{eff}}(N) = C \cdot (N / n_{\min})^\alpha$.

Usage

calibrate_L_eff(N, P1, sim_fpr, n_min = 60)

Arguments

Details

At each sample size, the implied $L_{\text{eff}}$ is:

$L_{\text{eff}} = \frac{\log(1 - \text{FPR}_{\text{sim}})} {\log(1 - P_1)}$

A log-linear model is then fit: $\log(L_{\text{eff}}) = \log(C) + \alpha \log(N / n_{\min})$.

Value

A list of class "leff_calibration" with components:

C

Numeric. Fitted constant.

alpha

Numeric. Fitted power-law exponent.

n_min

Integer. Minimum subgroup size.

r_squared

Numeric. R-squared of the log-linear fit.

data

Data frame with N, P1, sim_fpr, L_eff, and fitted values.

fit

The lm object.

Examples

# From binary threshold calibration results
cal <- calibrate_L_eff(
  N       = c(200, 500, 700, 1000),
  P1      = c(0.152, 0.109, 0.091, 0.071),
  sim_fpr = c(0.185, 0.220, 0.270, 0.280),
  n_min   = 60
)
print(cal)

Description

Compares the censoring distribution observed in the data used to build the DGM against the censoring generated by simulate_from_dgm. Reports censoring rates, time quantiles, KM-based median censoring times, and flags substantial discrepancies.

Usage

check_censoring_dgm(
  sim_data,
  dgm,
  treat_var = "treat_sim",
  rate_tol = 0.1,
  median_tol = 0.25,
  verbose = TRUE
)

Arguments

Details

The reference censoring distribution is derived from dgm$df_super, sampled with replacement from the data passed to generate_aft_dgm_flex(). Columns y (observed time) and event (event indicator) in df_super reflect the original observed censoring process on the DGM time scale.

The KM median censoring time is estimated by reversing the event indicator (1 - event), treating events as censored and censored observations as the event of interest. This gives a non-parametric estimate of the censoring time distribution unconfounded by event occurrence.

Common causes of discrepancy: (1) time-scale mismatch (DGM built on days, analysis_time in months); check exp(dgm$model_params$mu) against your analysis_time. (2) Large cens_adjust shifting censoring substantially from the fitted model. (3) Short analysis_time or time_eos making administrative censoring dominate the censoring process.

Value

Invisibly returns a named list. Elements are: rates (data frame of censoring rates overall and by arm); quantiles (data frame of censoring-time quantiles among censored subjects); km_medians (data frame of KM-based median censoring times); and flags (character vector of triggered warnings, empty if none).

See Also

simulate_from_dgm, generate_aft_dgm_flex

Examples

## Not run: 
dgm <- setup_gbsg_dgm(model = "null", verbose = FALSE)
sim_data <- simulate_from_dgm(dgm, n = 200)
check_censoring_dgm(sim_data, dgm = dgm)

## End(Not run)

Description

Separates successful results from errors in the output of a foreach loop run with .errorhandling = "pass". Successful results are combined via rbindlist (with an rbind fallback if data.table is not available) and failed tasks are counted, reported, and optionally escalated to a hard error.

Usage

collect_results(
  raw_list,
  label = "",
  stop_on_all_fail = TRUE,
  max_error_messages = 3L,
  keep_diagnostics = FALSE
)

Arguments

Details

This replaces the fragile pattern of calling nrow() or rbind() directly on a raw foreach result list, which fails silently when error objects are mixed in with data-frame rows.

Design rationale for the reporting strategy:

• A warning is emitted so that calling code (and tryCatch / withCallingHandlers blocks) can capture the failure condition programmatically.

• The error summary is also printed with cat to stderr, because warning and message output is sometimes buffered or suppressed inside future/foreach parallel workers or when options(warn = -1) is set by user code upstream. The cat() output is the visible-no-matter-what fallback.

• When stop_on_all_fail = TRUE (the default), a completely empty result list is escalated to an error rather than silently returning an empty data frame. This avoids downstream code proceeding with zero rows and producing misleading summaries. Calibration code that may legitimately tolerate some failures (e.g. Monte Carlo FPR estimation) can set stop_on_all_fail = FALSE.

The returned object is always a data frame – never NULL – because attributes (n_failed, n_total) attached to NULL are silently discarded by R, whereas attributes on an empty data frame are preserved and inspectable by the caller.

Value

A data frame of row-bound successful results. Always a data frame; never NULL – an empty data frame is returned when every task failed and stop_on_all_fail = FALSE. Two attributes are always attached:

n_failed

Integer. Number of tasks that returned errors.

n_total

Integer. Total number of tasks (successes + failures).

When keep_diagnostics = TRUE, a third attribute "diagnostics" may be attached – a named list of per-replicate diagnostic objects, keyed by "sim<ID>" when the source element carries an "sim_id" attribute (run_simulation_analysis sets this automatically).

Preserving per-replicate diagnostics

data.table::rbindlist() drops per-element attributes during row-binding. When the caller has populated each replicate with heavy diagnostic objects (e.g. via run_simulation_analysis with a non-empty keep argument), set keep_diagnostics = TRUE to retain them. The diagnostics attribute on each input element must carry attr(..., "sim_id") for stable keying; that is set automatically by run_simulation_analysis().

See Also

reset_workers for resetting the future plan and releasing worker memory between simulation phases.

Examples

library(foreach)
library(doFuture)
future::plan("sequential")

raw <- foreach(i = 1:5, .errorhandling = "pass") %dofuture% {
  if (i == 3L) stop("simulated failure")
  data.frame(task = i, value = rnorm(1))
}
res <- collect_results(raw, label = "demo")
# Warning: demo: 1 of 5 tasks failed.
nrow(res)
# [1] 4
attr(res, "n_failed")
# [1] 1

Description

Generates and compares detection probability curves for multiple subgroup sample sizes.

Usage

compare_detection_curves(
  n_sg_values,
  prop_cens = 0.3,
  hr_threshold = 1.25,
  hr_consistency = 1,
  theta_range = c(0.5, 3),
  n_points = 40L,
  verbose = TRUE
)

Arguments

Value

A data.frame with all curves combined, including n_sg as a factor.

Examples

## Not run: 
comparison <- compare_detection_curves(
  n_sg_values = c(40, 60, 80, 100),
  prop_cens = 0.2
)

# Plot with ggplot2
library(ggplot2)
ggplot(comparison, aes(x = theta, y = probability, color = factor(n_sg))) +
  geom_line(linewidth = 1) +
  labs(x = "True HR", y = "P(Detect)", color = "n_sg") +
  theme_minimal()

## End(Not run)

Description

Cross-tabulates subgroup assignments from two analyses and produces a formatted gt summary table with concordance statistics. Accepts forestsearch result objects, GRF result lists, or raw indicator vectors. Non-detection (no subgroup found) is handled automatically: all subjects are classified to the complement.

Usage

compare_fs_subgroups(
  fs1,
  fs2,
  label1 = "Analysis 1",
  label2 = "Analysis 2",
  id.name = NULL,
  sg0_label = NULL,
  sg1_label = NULL,
  subgroup_notation = c("harm", "benefit"),
  title = "Subgroup Membership Comparison",
  subtitle = NULL,
  font_size = 12
)

Arguments

Value

A list with:

table

A gt cross-tabulation table.

crosstab

The raw cross-tabulation matrix.

concordance

Agreement, kappa, counts.

membership

Per-subject assignments from both analyses.

Examples

## Not run: 
# Harm search: Cox vs Poisson (both forestsearch objects)
result <- compare_fs_subgroups(fs_cox, fs_pois,
  label1 = "Cox PH", label2 = "Poisson",
  subgroup_notation = "harm")

# Benefit search: FS vs GRF (handles non-detection automatically)
result <- compare_fs_subgroups(fs_result, grf_result,
  label1 = "ForestSearch", label2 = "GRF",
  subgroup_notation = "benefit")

# Raw vectors
result <- compare_fs_subgroups(
  c(rep(0L, 30), rep(1L, 70)),
  c(rep(0L, 25), rep(1L, 75)),
  label1 = "Method A", label2 = "Method B")

## End(Not run)

Description

Performs comprehensive comparison of multiple survreg models including convergence checking, information criteria comparison, and model selection.

Usage

compare_multiple_survreg(
  ...,
  model_names = NULL,
  verbose = TRUE,
  criteria = c("AIC", "BIC")
)

Arguments

Value

A list of class "multi_survreg_comparison" containing:

models

Named list of input models

convergence

Convergence status for each model

comparison

Model comparison statistics

rankings

Model rankings by different criteria

best_model

Name of the best model

recommendation

Text recommendation

Examples

## Not run: 
fit1 <- survreg(Surv(time, status) ~ x, dist = "weibull")
fit2 <- survreg(Surv(time, status) ~ x, dist = "lognormal")
comparison <- compare_multiple_survreg(fit1, fit2)

## End(Not run)

Description

Runs forestsearch once per combination of sg_focus[i] + selection_rule[i] (tuple semantics; the two vectors are paired element-by-element and must have the same length). All other arguments are held fixed across runs and passed through via ....

Usage

compare_selection_rules(
  df.analysis,
  sg_focus,
  selection_rule,
  compute_cis = TRUE,
  n_splits = 1000L,
  ci_seed = 1L,
  plot_xlim = NULL,
  show_band = TRUE,
  combo_labels = NULL,
  verbose = TRUE,
  ...
)

Arguments

Details

For each combination the wrapper:

1. Captures the in-flight console output from forestsearch (including the pre-consistency Candidate Evaluation Preview and post-consistency Candidate Evaluation Summary when show_candidate_summary = TRUE is passed).

2. Optionally computes the frontier CIs via compute_frontier_cis.

3. Builds the Pareto-frontier plot via plot_pareto_frontier.

The captured stdout is returned per combination so the caller can cat() it later for side-by-side inspection in a Quarto chunk. Plot objects are returned individually and (when patchwork is installed) composed into a single side-by-side figure. The full fs object for each combination is returned so downstream diagnostics (pareto_frontier_table, explain_pareto_selection) can be run without re-fitting.

Value

A list with class "forestsearch_comparison" and components:

combos

data.frame with columns sg_focus, selection_rule, label.

fs

Named list of forestsearch result objects, one per combination. NULL entries mark fits that errored.

ci_tab

Named list of frontier-CI tables (or all NULL if compute_cis = FALSE).

plots

Named list of ggplot objects, one per combination.

plot_grid

A single patchwork object placing the plots side by side, or NULL if patchwork isn't installed.

plot_combined

A single ggplot composing all configurations onto one Pareto plot, with each winner labeled S1: <combo_label>, S2: <combo_label>, etc. Returned only when ALL successful configurations share an identical passing set (see plot_pareto_combined for the equality criterion); NULL otherwise. See plot_combined_subsets below for the per-subset fallback when some – but not all – configurations match.

plot_combined_subsets

Named list of ggplot objects, one per equivalence class of configurations sharing an identical passing set. Only groups of size $\ge 2$ are included (a singleton group has no peer to combine with). Names are derived from the shared sg_focus when constant within a group (e.g., "effMaxSG"); from the concatenated combo labels otherwise. When all valid combos share one passing set this list contains exactly one element, equal to plot_combined. The list is empty when no two combos match.

combined_skip_reason

Character scalar or NULL. When plot_combined is NULL because the equality precondition failed, this captures the specific reason – size mismatch, definition-set mismatch, or value drift on a named column – as a single string. NULL when the combined plot was built or when fewer than two fits succeeded.

console

Named list of character vectors – the captured stdout from each forestsearch run (full output).

diagnostics

Named list of per-combo slices of the captured output: each element is a list with character-scalar fields preview (the CANDIDATE EVALUATION PREVIEW block), summary (the CANDIDATE EVALUATION SUMMARY block), and full (the entire captured stdout, same content as console[[i]] collapsed to a single string). preview and summary are NA_character_ when the corresponding banner is absent (e.g.\ show_candidate_summary was FALSE, or the run errored before reaching it). This slice is what a comparison document should display as the primary diagnostic; full (or console[[i]]) remains available for an expandable / on-demand view of the complete run output.

errors

Named list of error messages (or NULL) for combos that failed.

Tuple semantics

With sg_focus = c("hrMaxSG", "hrMaxSG") and selection_rule = c("pareto", "both"), two runs are performed: one with (hrMaxSG, pareto) and one with (hrMaxSG, both). The length-mismatch case raises an error.

Console capture

All stdout from each forestsearch() call is captured (not just the diagnostic tables) via capture.output(). The PREVIEW / SUMMARY blocks are clearly delineated by banner separators within the captured text. Use cat(out$console[[i]]) to replay the output in a Quarto chunk; results = "asis" works but is not required since the captured text already includes newlines.

See Also

forestsearch, plot_pareto_frontier, plot_pareto_combined, pareto_frontier_table, explain_pareto_selection, compute_frontier_cis.

Examples

## Not run: 
out <- compare_selection_rules(
  df.analysis      = actg_df,
  sg_focus         = c("hrMaxSG", "hrMaxSG"),
  selection_rule   = c("pareto",  "both"),
  # All other forestsearch args:
  confounders.name = confounders.name,
  outcome.name     = adverse_outcome,
  treat.name       = treat.name,
  id.name          = id.name,
  outcome_type     = "continuous",
  effect_measure   = "MD",
  adverse_outcome  = TRUE,
  hr.threshold     = 10,
  hr.consistency   = 5,
  pconsistency.threshold = 0.90,
  show_candidate_summary = TRUE,
  details          = TRUE
)

# Inspect captured console output for combo 1
cat(out$console[[1]])

# Side-by-side plot
print(out$plot_grid)

# Use fs object for downstream diagnostics
pareto_frontier_table(out$fs[[1]], ci_table = out$ci_tab[[1]])
explain_pareto_selection(out$fs[[2]], ci_table = out$ci_tab[[2]])

## End(Not run)

Description

Builds the per-subgroup comparison frame between two run_subgroup_sims() results (or RDS payloads loaded from the vignettes – any list carrying sim_hrs / sim_ubs / sim_ns matrices with identical column names). For each input the per-design statistics replicate the vignettes' Section 6.6.1 definitions exactly: tail probabilities are percentages denominated over converged fits (na.rm = TRUE), medians are unconditional, N is the across-trial mean subgroup size, and structurally empty subgroups (all-NA columns) are mapped from NaN to NA.

Usage

compare_subgroup_sims(
  x,
  y,
  expect_designs = NULL,
  suffixes = c("_r", "_f"),
  est_thresholds = NULL,
  ub_thresholds = NULL
)

Arguments

Details

Column names are structural and retained across outcome types (hr05, hr1, ub2, ub3, mHR, mUB), exactly as sim_hrs serves generic duty on GLM results: for subgroup_glm() fits they hold estimate-scale statistics at the resolved thresholds. Threshold resolution: explicit arguments win; otherwise the inputs' effect metadata supplies them; otherwise the HR legacy c(0.5, 1.0) / c(2, 3). An NA threshold yields an all-NA column. The two inputs must carry compatible metadata (both none, or the same effect measure) – comparing an MD study against an HR study is a hard error, raised before the panel-alignment guard so the categorical incompatibility is reported first.

Two validations guard the row-wise alignment (formerly the memo's Guard 1 and Guard 2): the subgroup panels must be identical, and, when expect_designs is supplied, each input's design label must match it exactly.

The medians here use median(), matching the memo; they are equal to summary.subgroup_sims()'s type-7 quantile(..., 0.50) values up to floating-point formulation ((a + b) / 2 versus a + 0.5 * (b - a)), which is why the memo's numbers and the vignettes' tables have always agreed to the printed digit.

Value

A plain data.frame (one row per subgroup, matching the memo's static-mode frame type) with columns subgroup, then for each of N, ub2, ub3, mUB, hr05, hr1, mHR the x and y columns interleaved. The inputs' n_sims and design values are attached as attributes "n_sims" and "designs"; when the inputs carry effect metadata it is attached as attribute "effect" (informational; attributes are dropped by most transformations).

See Also

run_subgroup_sims(), summary.subgroup_sims(), subgroup_glm()

Examples

## Not run: 
pl_r <- readRDS("results/extreme_sims_resample_10000_payload.rds")
pl_f <- readRDS("results/extreme_sims_fixed_10000_payload.rds")
cmp <- compare_subgroup_sims(pl_r, pl_f,
                             expect_designs = c("resample", "fixed"))

## End(Not run)

Description

Calculates the probability that a true subgroup with given hazard ratio will be detected using the ForestSearch consistency-based criteria.

Usage

compute_detection_probability(
  theta,
  n_sg,
  prop_cens = 0.3,
  hr_threshold = 1.25,
  hr_consistency = 1,
  method = c("cubature", "monte_carlo"),
  n_mc = 100000L,
  tol = 1e-04,
  verbose = FALSE
)

Arguments

Details

This function computes P(detect | theta) using the asymptotic normal approximation for the log hazard ratio estimator. The detection criterion is based on ForestSearch's split-sample consistency evaluation:

1. The subgroup HR estimate must exceed hr_threshold on average

2. Each split-half must individually exceed hr_consistency

The approximation assumes:

• Large sample sizes (CLT applies)

• Var(log(HR)) ~ 4/d for the full subgroup (d = total subgroup events); each 50/50 split-half (d/2 events) has variance 8/d

• Independence between split-halves (conditional on true effect)

Value

If theta is scalar, returns a single probability. If theta is a vector, returns a data.frame with columns: theta, probability.

Examples

## Not run: 
# Single HR value
prob <- compute_detection_probability(
  theta = 1.5,
  n_sg = 60,
  prop_cens = 0.2,
  hr_threshold = 1.25
)

# Vector of HR values for power curve
hr_values <- seq(1.0, 2.5, by = 0.1)
results <- compute_detection_probability(
  theta = hr_values,
  n_sg = 60,
  prop_cens = 0.2,
  hr_threshold = 1.25,
  verbose = TRUE
)

# Plot detection probability curve
plot(results$theta, results$probability, type = "l",
     xlab = "True HR", ylab = "P(detect)")

## End(Not run)

Description

Extends the Section 2.1 approximation (Leon et al., 2024, eq. 3) to arbitrary GLM outcomes via the effective information $d_{\text{eff}}$. This is the GLM generalization of compute_detection_probability.

Usage

compute_detection_probability_glm(
  theta,
  d_eff,
  c1 = 1.25,
  c2 = 1,
  effect_scale = c("ratio", "difference"),
  method = c("cubature", "monte_carlo"),
  n_mc = 500000L,
  seed = 42L,
  verbose = FALSE
)

Arguments

Details

The detection criterion (Leon et al., 2024, Section 2.1):

1. Screening: $W_1 + W_2 \ge 2k_1$

2. Consistency: $\min(W_1, W_2) \ge k_2$

with $W_1, W_2 \sim N(\beta, 8/d_{\text{eff}})$ independently.

For ratio-scale effects, $\beta = \log(\theta)$, $k_j = \log(c_j)$. For difference-scale effects, $\beta = \theta$, $k_j = c_j$.

Value

If theta is scalar, a single probability. If vector, a data.frame with columns theta, probability.

See Also

compute_detection_probability for the survival-specific version.

Examples

# Survival: reproduce Figure 2
d <- d_eff_survival(n_sg = 60, prop_cens = 0.45)
compute_detection_probability_glm(theta = 2.0, d_eff = d)

# Binary: OR = 2.0, event rate 30%
d <- d_eff_binary(n_sg = 100, p_event = 0.30)
compute_detection_probability_glm(theta = 2.0, d_eff = d)

# Count: IRR = 2.0, 80 total events
d <- d_eff_count(total_events = 80)
compute_detection_probability_glm(theta = 2.0, d_eff = d)

# Continuous: mean difference = 0.5, sigma = 1.5
d <- d_eff_continuous(n_sg = 100, sigma_y = 1.5)
compute_detection_probability_glm(
  theta = 0.5, d_eff = d,
  c1 = 0.2, c2 = 0.0,
  effect_scale = "difference"
)

# Power curve
d <- d_eff_binary(100, 0.30)
result <- compute_detection_probability_glm(
  theta = seq(0.5, 3.0, by = 0.1), d_eff = d
)
plot(result$theta, result$probability, type = "l")

Description

Calculates Controlled Direct Effect (CDE) hazard ratios from the super-population potential outcomes (theta_0, theta_1) and attaches them to the DGM's hazard_ratios list. This enables automatic CDE detection by build_estimation_table.

Usage

compute_dgm_cde(dgm, harm_col = NULL)

Arguments

Details

The CDE for subgroup S is defined as:

CDE(S) = mean(exp(theta_1[S])) / mean(exp(theta_0[S]))

which is the ratio of average hazard contributions on the natural scale. This differs from the AHR (exp(mean(loghr_po))) due to Jensen's inequality. In the notation of Leon et al. (2024), CDE corresponds to theta-ddagger.

The function detects the subgroup indicator column automatically, checking for flag.harm, flag_harm, and H in the super-population data frame.

Value

The DGM object with CDE values added to dgm$hazard_ratios (CDE, CDE_harm, CDE_no_harm) and to top-level fields (dgm$CDE, dgm$cde_H, dgm$cde_Hc).

See Also

build_estimation_table, get_dgm_hr

Examples

## Not run: 
dgm <- create_gbsg_dgm(model = "alt", k_inter = 2.0)
dgm <- compute_dgm_cde(dgm)
dgm$hazard_ratios$CDE_harm   # theta-ddagger(H)
dgm$hazard_ratios$CDE         # theta-ddagger overall

## End(Not run)

Description

For each candidate subgroup on fs$grp.consistency$out_sg$pareto_frontier, computes three 95\

Usage

compute_frontier_cis(
  fs,
  n_splits = 1000L,
  ci_level = 0.95,
  seed = NULL,
  verbose = FALSE
)

Arguments

Details

• Naive CI: full-sample Wald CI from a Cox or GLM model refit on the subgroup's data. This is the textbook CI that ignores the multiple-testing / subgroup-search context. Anti-conservative by construction.

• Split CI: a 95\ estimate, with SE estimated by the empirical standard deviation of the $2S$ individual half-sample effect estimates produced by n_splits random 50/50 splits. This is a half-jackknife SE (Shao1996). Naming: the $\sim$ in column labels marks this as a resampling-derived approximation, not a model-based CI.

• FSBC-mimic CI: a bias-corrected interval inspired by the bootstrap bias-correction algorithm of Leon2024fs (eq 7), with the selection-on-bootstrap-data term $\eta_b^*(\hat H_b^*)$ zeroed out (the selected subgroup is treated as fixed across half-jackknife replicates). The bias-corrected estimate is $\hat\beta^{\mathrm{FSBC}} = 2\hat\beta - \overline{\hat\beta^{(h)}}$ where $\overline{\hat\beta^{(h)}}$ is the mean of the $2S$ half-sample estimates; the SE is the same half-jackknife SE as for the Split CI. See “Details — FSBC-mimic interpretation” below.

All three CIs are computed post-hoc on the returned forestsearch object; this function does not modify any internal state.

Value

A data.table keyed by frontier row index (m) with columns:

m

Candidate id (matches pareto_frontier$m).

estimate

Refit effect estimate on the natural scale.

naive_lcl, naive_ucl

Naive Wald CI bounds.

split_sd

Half-jackknife SE estimate: empirical SD of the $2S$ per-half effect estimates on the SE scale (log for ratio measures, natural for differences).

split_lcl, split_ucl

Split-derived CI bounds.

fsbc_estimate

Bias-corrected effect estimate on the natural scale ($2\hat\beta - \overline{\hat\beta^{(h)}}$ on the SE scale, then back-transformed).

fsbc_lcl, fsbc_ucl

FSBC-mimic CI bounds.

n_valid_splits

Number of splits that produced valid per-half estimates.

Rows where the refit fails (e.g., zero events in one arm) have NA in all CI columns.

Details — Half-jackknife SE

For each of n_splits random 50/50 splits we fit the effect model separately on each half, obtaining $2S$ per-half estimates $\{\hat\beta^{(s,h)} : s=1,\ldots,S; h=1,2\}$. The half-jackknife SE estimator is

$\widehat{\mathrm{SE}} = \mathrm{sd}\bigl\{\hat\beta^{(s,h)}\bigr\},$

i.e., the plain empirical SD of the $2S$ estimates. No $\sqrt{2}$ correction is applied: simulation studies confirm that this estimator's CI achieves nominal coverage in a no-selection setting (see the package's pressure-test script). Earlier versions of this function used the average $(\hat\beta^{(s,1)} + \hat\beta^{(s,2)})/2$ per split; that approach is invalid because the two halves are complementary and their average is nearly constant.

Details — FSBC-mimic interpretation

The FSBC-mimic CI is a quick approximation to the bootstrap bias-corrected estimator of Leon2024fs. Substituting the $2S$ per-half estimates for bootstrap replicates, and treating the selected subgroup $\hat H$ as fixed (i.e., omitting the $\eta_b^*(\hat H_b^*) = \hat\beta_b^*(\hat H_b^*) - \hat\beta(\hat H_b^*)$ term that captures variability in FS selection on the bootstrap data), the bias-corrected estimator from paper eq (7) collapses to

$\hat\beta^{\mathrm{FSBC}}(\hat H) = \hat\beta(\hat H) - \frac{1}{2S}\sum_{s,h}\bigl[\hat\beta^{(s,h)} - \hat\beta(\hat H)\bigr] = 2\hat\beta(\hat H) - \overline{\hat\beta^{(h)}}.$

The SE is the same half-jackknife SE as for the Split CI.

The FSBC-mimic CI is not a full FSBC CI: it omits the selection-on-bootstrap-data variance. Use this column as a quick bias-corrected diagnostic, not for inference. For full FSBC CIs use forestsearch_bootstrap_dofuture.

See Also

pareto_frontier_table, frontier_member_flags, plot_pareto_frontier.

Description

Aggregates scores by leaf node and calculates treatment effect differences

Usage

compute_node_metrics(data, dr.scores, tree, X, n.min)

Arguments

Value

Data frame with node metrics

Examples

## Not run: 
# compute_node_metrics() is called internally by grf.subg.harm.survival().
# See grf.subg.harm.survival() for the standard entry point.

## End(Not run)

Description

Computes normal-approximation confidence intervals $\hat\beta_j \pm z_{1 - \alpha/2} \sqrt{\widehat{V}_{jj}}$ using the sandwich (or Cox model-based) variance returned by dina_fit().

Usage

## S3 method for class 'dina'
confint(object, parm, level = 0.95, ...)

Arguments

Value

a matrix with two columns giving lower and upper limits, named in standard confint() style (e.g. "2.5 %", "97.5 %").

Description

Approximates the forestsearch consistency rate for a candidate subgroup without repeated sample splits and refits. The subgroup is fit once; the random-split pair of half-sample treatment effects is represented as ⁠beta_hat +/- D⁠, a multiplier sum of the per-subject treatment dfbeta contributions (see the file header). Handles the Cox (survival) path and the GLM paths (binary OR/RR/RD, continuous MD, count IRR).

Usage

consistency_resample(
  df,
  hr.consistency = 1,
  consistency_threshold = NULL,
  comparison_threshold = NULL,
  outcome_type = c("survival", "binary", "continuous", "count"),
  effect_measure = NULL,
  method = c("both", "closed", "mc"),
  multiplier = c("rademacher", "normal", "poisson"),
  draws = 1000L,
  tte.name = "Y",
  event.name = "Event",
  treat.name = "Treat",
  outcome.name = "Y",
  offset.name = NULL,
  cox_init = 0,
  adjust_covariates = NULL,
  adverse_outcome = TRUE,
  seed = NULL
)

Arguments

Details

multiplier = "rademacher" (default) reproduces the 50/50 Bernoulli split of run_single_consistency_split(); "normal" is the Lin (1997) Gaussian multiplier and "poisson" a centred unit-Poisson (DBP) multiplier.

Value

List with rate_closed, rate_mc, beta_hat, sigma_D, n, d, delta, measure; all-NA rates (with available summaries) if the fit fails or the configuration is unsupported.

See Also

run_single_consistency_split(), consistency_resample_compare()

Examples

## Not run: 
# survival (Cox)
consistency_resample(sub, hr.consistency = 1.0)
# binary odds ratio
consistency_resample(sub, outcome_type = "binary", effect_measure = "OR",
                     consistency_threshold = 1.0,
                     outcome.name = "Y", treat.name = "Treat")
# count incidence-rate ratio
consistency_resample(sub, outcome_type = "count", effect_measure = "IRR",
                     consistency_threshold = 1.0,
                     outcome.name = "events", offset.name = "fu_time")

## End(Not run)

Description

Calls the package's run_single_consistency_split() R_true times and compares the Monte Carlo rate to the closed-form and multiplier approximations from consistency_resample(). Survival path only; for the GLM paths use the standalone GLM validator.

Usage

consistency_resample_compare(
  df,
  hr.consistency = 1,
  R_true = 400L,
  draws = 1000L,
  multiplier = "rademacher",
  tte.name = "Y",
  event.name = "Event",
  treat.name = "Treat",
  cox_init = 0,
  adjust_covariates = NULL,
  seed = NULL
)

Arguments

Value

One-row data.frame with n, d, beta_hat, sigma_D, rate_true, valid_true, rate_closed, rate_mc, err_closed, err_mc.

See Also

consistency_resample(), run_single_consistency_split()

Description

Counts the number of times an ID appears in a bootstrap sample.

Usage

count_boot_id(x, dfb)

Arguments

Value

Integer count of occurrences.

Examples

df_boot <- data.frame(id = c(1, 2, 1, 3, 1), id_boot = 1:5)
count_boot_id(1, df_boot)  # returns 3
count_boot_id(4, df_boot)  # returns 0

Description

This wrapper function combines Cox spline fitting with comprehensive visualization of Average Hazard Ratios (AHRs) and Controlled Direct Effects (CDEs) as described in the MRCT subgroups analysis documentation.

Usage

cox_ahr_cde_analysis(
  df,
  tte_name = "os_time",
  event_name = "os_event",
  treat_name = "treat",
  z_name = "biomarker",
  loghr_po_name = "loghr_po",
  theta1_name = "theta_1",
  theta0_name = "theta_0",
  spline_df = 3,
  alpha = 0.2,
  hr_threshold = 0.7,
  plot_style = c("combined", "separate", "grid"),
  plot_select = c("all", "profile_ahr", "ahr_only"),
  save_plots = FALSE,
  output_dir = "plots",
  verbose = TRUE
)

Arguments

Value

List of class "cox_ahr_cde" containing:

cox_fit

Results from cox_cs_fit function.

ahr_results

AHR calculations for different subgroup definitions.

cde_results

CDE calculations if theta variables available.

optimal_cutpoint

Optimal biomarker cutpoint, or NULL when hr_threshold is NULL.

subgroup_stats

Statistics for recommended and questionable subgroups, or overall-only when hr_threshold is NULL.

data

List with z_values, loghr_po, and subgroup assignments.

Examples

## Not run: 
# With threshold - full subgroup analysis
results <- cox_ahr_cde_analysis(
  df = df_large, z_name = "z_bm",
  hr_threshold = 1.25, plot_style = "grid"
)

# Without threshold - pure AHR/CDE curves
results <- cox_ahr_cde_analysis(
  df = df_large, z_name = "z_bm",
  hr_threshold = NULL, plot_style = "grid"
)

# Compact two-panel without threshold
results <- cox_ahr_cde_analysis(
  df = df_large, z_name = "z_bm",
  hr_threshold = NULL,
  plot_select = "profile_ahr"
)

# Single AHR panel only
results <- cox_ahr_cde_analysis(
  df = df_large, z_name = "z_bm",
  hr_threshold = 1.25,
  plot_select = "ahr_only"
)

## End(Not run)

Description

Estimates treatment effects as a function of a continuous covariate using a Cox proportional hazards model with natural cubic splines. The function models treatment-by-covariate interactions to detect effect modification.

Usage

cox_cs_fit(
  df,
  tte_name = "os_time",
  event_name = "os_event",
  treat_name = "treat",
  strata_name = NULL,
  z_name = "bm",
  alpha = 0.2,
  spline_df = 3,
  z_max = Inf,
  z_by = 1,
  z_window = 0,
  z_quantile = 0.9,
  show_plot = TRUE,
  plot_params = NULL,
  truebeta_name = NULL,
  verbose = TRUE
)

Arguments

Details

Model Structure

The function fits:

$h(t|Z,A) = h_0(t) \exp(\beta_0 A + f(Z) + g(Z) \cdot A)$

Where:

• A is treatment (0/1)

• Z is the continuous effect modifier

• f(Z) is modeled with natural splines (main effect)

• g(Z) is modeled with natural splines (interaction)

• The log hazard ratio is: $\beta(Z) = \beta_0 + g(Z)$

Plot Parameters

The plot_params argument accepts a list with:

• xlab: x-axis label

• main_title: plot title

• ylimit: y-axis limits c(min, max)

• y_pad_zero: padding below zero line

• y_delta: extra space for count labels

• cex_legend: legend text size

• cex_count: count text size

• show_cox_primary: show standard Cox estimate line

• show_null: show null effect line (log(HR)=0)

• show_target: show target effect line (e.g., log(0.80))

Value

List containing:

z_profile

Vector of z values where treatment effect is estimated

loghr_est

Point estimates of log(HR) at each z value

loghr_lower

Lower confidence bound

loghr_upper

Upper confidence bound

se_loghr

Standard errors of log(HR) estimates

counts_profile

Number of observations near each z value

cox_primary

Log(HR) from standard Cox model (no interaction)

model_fit

The fitted coxph model object

spline_basis

The natural spline basis object

Examples

## Not run: 
# Simulate data
set.seed(123)
df <- data.frame(
  os_time = rexp(500, 0.01),
  os_event = rbinom(500, 1, 0.7),
  treat = rbinom(500, 1, 0.5),
  bm = rnorm(500, 50, 10)
)

# Fit model
result <- cox_cs_fit(df, z_name = "bm", alpha = 0.20)

# Custom plotting
result <- cox_cs_fit(
  df,
  z_name = "bm",
  plot_params = list(
    xlab = "Biomarker Level",
    main_title = "Treatment Effect by Biomarker",
    cex_legend = 1.2
  )
)

## End(Not run)

Description

Called in analyze_subgroup() <– SG_tab_estimates

Usage

cox_summary(
  Y,
  E,
  Treat,
  Strata = NULL,
  use_strata = !is.null(Strata),
  return_format = c("formatted", "numeric")
)

Arguments

Details

Calculates hazard ratio and confidence interval for a subgroup using Cox regression. Optimized version with reduced overhead and better error handling.

Value

Character string with formatted HR and CI (or numeric vector if return_format="numeric").

Examples

## Not run: 
library(survival)
cox_summary(
  Y     = gbsg$rfstime / 30.4375,
  E     = gbsg$status,
  Treat = gbsg$hormon
)

## End(Not run)

Description

Wrapper function to create a data generating mechanism (DGM) for MRCT simulation scenarios using generate_aft_dgm_flex.

Usage

create_dgm_for_mrct(
  df_case,
  model_type = c("alt", "null"),
  log_hrs = NULL,
  confounder_var = NULL,
  confounder_effect = NULL,
  include_regA = TRUE,
  verbose = FALSE
)

Arguments

Details

Model Types

alt

Alternative hypothesis: Treatment effect varies by biomarker level (heterogeneous treatment effect). Default log_hrs create HR ranging from 2.0 (harm) to 0.5 (benefit) across biomarker range

null

Null hypothesis: Uniform treatment effect regardless of biomarker level. Default log_hrs = log(0.7) uniformly

Confounder Effects

By default, NO prognostic confounder effect is forced. The confounder_var and confounder_effect parameters allow optionally specifying ANY baseline covariate to have a fixed prognostic effect in the outcome model.

The regA variable (region indicator) is included as a factor by default but without a forced effect - its coefficient is estimated from data.

Value

An object of class "aft_dgm_flex" for use with simulate_from_dgm and mrct_region_sims

See Also

generate_aft_dgm_flex for underlying DGM creation mrct_region_sims for running simulations with the DGM

Examples

## Not run: 
# Prepare data
df_case <- read.csv("data/dfsynthetic.csv")
df_case$regA <- df_case$region_asia

# Alternative hypothesis (heterogeneous treatment effect)
dgm_alt <- create_dgm_for_mrct(
  df_case = df_case,
  model_type = "alt",
  log_hrs = log(c(3, 1.25, 0.50)),
  verbose = TRUE
)

# Null hypothesis (uniform effect)
dgm_null <- create_dgm_for_mrct(
  df_case = df_case,
  model_type = "null",
  verbose = TRUE
)

# With forced confounder effect
dgm_conf <- create_dgm_for_mrct(
  df_case = df_case,
  model_type = "alt",
  confounder_var = "prior_treat",
  confounder_effect = log(1.5),
  verbose = TRUE
)

## End(Not run)

Description

Creates a forestploter theme with parameters that control overall plot sizing and appearance. This is the primary way to control how large the forest plot renders.

Usage

create_forest_theme(
  base_size = 10,
  scale = 1,
  row_padding = NULL,
  ci_pch = 15,
  ci_lwd = NULL,
  ci_Theight = NULL,
  ci_col = "black",
  header_fontsize = NULL,
  body_fontsize = NULL,
  footnote_fontsize = NULL,
  footnote_col = "darkcyan",
  title_fontsize = NULL,
  cv_fontsize = NULL,
  cv_col = "gray30",
  refline_lwd = NULL,
  refline_lty = "dashed",
  refline_col = "gray30",
  vertline_lwd = NULL,
  vertline_lty = "dashed",
  vertline_col = "gray20",
  arrow_type = "closed",
  arrow_col = "black",
  summary_fill = "black",
  summary_col = "black"
)

Arguments

Details

The base_size parameter is the primary way to control plot size. When you change base_size, the following are automatically scaled:

• All font sizes (body, header, footnote, CV, title)

• Row padding (vertical and horizontal)

• CI line width and T-bar height

• Reference and vertical line widths

The scaling formula uses base_size = 10 as the reference point:

• base_size = 10: Default sizing

• base_size = 12: 20% larger

• base_size = 14: 40% larger

• base_size = 16: 60% larger

You can override any individual parameter by specifying it explicitly.

The theme does NOT set row background colors - those are determined automatically by plot_subgroup_results_forestplot() based on row types (ITT, reference, posthoc, etc.).

Value

A list of class "fs_forest_theme" containing all theme parameters.

See Also

plot_subgroup_results_forestplot, render_forestplot

Examples

## Not run: 
# Simple: just increase base_size for larger plot
large_theme <- create_forest_theme(base_size = 14)

# Or use scale for quick adjustment
large_theme <- create_forest_theme(base_size = 10, scale = 1.4)

# Fine-tune specific elements
custom_theme <- create_forest_theme(
  base_size = 14,
  cv_fontsize = 12,  # Override auto-calculated CV font size
  ci_lwd = 2.5       # Override auto-calculated CI line width
)

# Use with plot_subgroup_results_forestplot
result <- plot_subgroup_results_forestplot(
  fs_results = list(fs.est = fs, fs_bc = fs_bc),
  df_analysis = df.analysis,
  outcome.name = "time",
  event.name = "status",
  treat.name = "treatment",
  theme = large_theme
)

render_forestplot(result)

## End(Not run)

Description

Creates a data generating mechanism (DGM) for survival simulations based on the German Breast Cancer Study Group (GBSG) dataset. Supports heterogeneous treatment effects via treatment-subgroup interactions.

Usage

create_gbsg_dgm(
  model = c("alt", "null"),
  k_treat = 1,
  k_inter = 1,
  k_z3 = 1,
  z1_quantile = 0.25,
  n_super = DEFAULT_N_SUPER,
  cens_type = c("weibull", "uniform"),
  use_rand_params = FALSE,
  seed = SEED_BASE,
  verbose = FALSE
)

Arguments

Details

This version is aligned with generate_aft_dgm_flex() and calculate_hazard_ratios() methodology, computing individual-level potential outcomes and average hazard ratios (AHR).

Subgroup Definition

The harm subgroup H is defined as: z1 = 1 AND z3 = 1, where:

• z1: Low estrogen receptor (ER <= 25th percentile by default)

• z3: Premenopausal status (meno == 0)

Model Specification

The AFT model uses covariates: treat, z1, z2, z3, z4, z5, and (for "alt") the interaction zh = treat * z1 * z3.

Interaction Effect (k_inter)

The k_inter parameter modifies the zh coefficient in the AFT model:

gamma[zh] <- k_inter * gamma[zh]

This affects the hazard ratio for the harm subgroup:

• HR(H) = exp(-gamma[treat]/sigma - gamma[zh]/sigma)

• HR(Hc) = exp(-gamma[treat]/sigma)

When k_inter = 0, HR(H) = HR(Hc) (no heterogeneity).

Alignment with generate_aft_dgm_flex

This function now computes:

• theta_0: Log-hazard contribution under control

• theta_1: Log-hazard contribution under treatment

• loghr_po: Individual causal log hazard ratio (theta_1 - theta_0)

• AHR metrics: exp(mean(loghr_po)) for overall and subgroups

Value

A list of class "gbsg_dgm" containing:

df_super_rand

Data frame with randomized super-population including potential outcomes (theta_0, theta_1, loghr_po)

hr_H_true

Empirical hazard ratio in harm subgroup (Cox-based)

hr_Hc_true

Empirical hazard ratio in complement subgroup (Cox-based)

hr_causal

Overall causal (ITT) hazard ratio (Cox-based)

AHR

Overall average hazard ratio (from loghr_po)

AHR_H_true

Average hazard ratio in harm subgroup

AHR_Hc_true

Average hazard ratio in complement subgroup

hazard_ratios

List matching generate_aft_dgm_flex output format

model_params

List with AFT model parameters (mu, sigma, gamma, etc.)

cens_params

List with censoring model parameters

subgroup_info

List with subgroup definitions and true factor names

analysis_vars

Character vector of analysis variable names

model_type

Character indicating "alt" or "null"

See Also

simulate_from_gbsg_dgm for generating data from the DGM calibrate_k_inter for finding k_inter to achieve target HR

Examples

## Not run: 
# Alternative hypothesis with default parameters
dgm_alt <- create_gbsg_dgm(model = "alt", verbose = TRUE)

# Null hypothesis
dgm_null <- create_gbsg_dgm(model = "null", verbose = TRUE)

# Custom subgroup HR via k_inter
dgm_custom <- create_gbsg_dgm(
  model = "alt",
  k_treat = 1.2,
  k_inter = 2.0,
  verbose = TRUE
)

# Access AHR metrics (aligned with generate_aft_dgm_flex)
dgm_alt$hazard_ratios$AHR_harm
dgm_alt$hazard_ratios$AHR_no_harm

## End(Not run)

Description

Fits a GLM within the specified subgroup and returns a one-row data frame containing the effect estimate, confidence interval, sample size, and event count. Supports binary, continuous, and count (with offset) outcomes.

Usage

create_glm_row(
  df,
  outcome.name,
  treat.name = "treat",
  effect_measure = "log_OR",
  offset.name = NULL,
  overdispersion = "none",
  conf.level = 0.95,
  min_arm_n = 5L,
  min_arm_events = 3L,
  verbose = FALSE
)

Arguments

Value

A one-row data frame with columns: N, n_treat, n_control, events_treat, events_control (for binary/count), est, lower, upper, se, effect_measure.

See Also

grf.subg.harm.glm