forestsearch: Exploratory Subgroup Identification in Clinical Trials with Survival Endpoints

Description

Implements statistical methods for exploratory subgroup identification in clinical trials with survival endpoints. Provides tools for identifying patient subgroups with differential treatment effects using machine learning approaches including Generalized Random Forests (GRF), LASSO regularization, and exhaustive combinatorial search algorithms. Features bootstrap bias correction using infinitesimal jackknife methods to address selection bias in post-hoc analyses. Designed for clinical researchers conducting exploratory subgroup analyses in randomized controlled trials, particularly for multi-regional clinical trials (MRCT) requiring regional consistency evaluation. Supports both accelerated failure time (AFT) and Cox proportional hazards models with comprehensive diagnostic and visualization tools. Methods are described in León et al. (2024) doi:10.1002/sim.10163.

Author(s)

Maintainer: Larry Leon larry.leon.05@post.harvard.edu

See Also

Useful links:

• https://github.com/larry-leon/forestsearch

• https://larry-leon.github.io/forestsearch/

• Report bugs at https://github.com/larry-leon/forestsearch/issues

Cross-Validation Subgroup Match Summary

Description

Summarizes the match between cross-validation subgroups and analysis subgroups.

Usage

CV_sgs(sg1, sg2, confs, sg_analysis)

Arguments

Value

List with indicators for any match, exact match, one match, and covariate-specific matches.

Convert Factor Code to Label

Description

Converts q-indexed codes to human-readable labels using the confs_labels mapping. Supports both full format ("q1.1", "q3.0") and short format ("q1", "q3"). Handles vector input via recursion.

Usage

FS_labels(Qsg, confs_labels)

Arguments

Value

Character. Human-readable label wrapped in braces, e.g., "{age <= 50}" or "!{age <= 50}" for complement. Returns the original code if no match is found.

Subgroup summary table estimates

Description

Returns a summary table of subgroup estimates (HR, RMST, medians, etc.).

Usage

SG_tab_estimates(
  df,
  SG_flag,
  outcome.name = "tte",
  event.name = "event",
  treat.name = "treat",
  strata.name = NULL,
  hr_1a = NA,
  hr_0a = NA,
  potentialOutcome.name = NULL,
  sg1_name = NULL,
  sg0_name = NULL,
  draws = 0,
  details = FALSE,
  return_medians = TRUE,
  est.scale = "hr"
)

Arguments

Value

Data frame of subgroup summary estimates.

Violin/Boxplot Visualization of HR Estimates

Description

Creates violin plots with embedded boxplots showing the distribution of hazard ratio estimates across simulations for different analysis populations. Supports symmetric trimming to handle extreme values that can distort the display when small subgroups produce very large HR estimates.

Usage

SGplot_estimates(
  df,
  label_training = "Training",
  label_testing = "Testing",
  label_itt = "ITT (stratified)",
  label_sg = "Testing (subgroup)",
  trim_fraction = NULL,
  ylim = NULL,
  show_summary = NULL,
  title = "Distribution of HR Estimates Across Simulations",
  subtitle = NULL
)

Arguments

Value

List with components:

dfPlot_estimates

data.table formatted for plotting, with a trimmed logical column when trimming is active

plot_estimates

ggplot2 object

trim_info

List of per-group trimming diagnostics (NULL when no trimming). Each element contains: n_total, n_trimmed, n_flagged, raw_mean, raw_sd, trimmed_mean, trimmed_sd, lower_bound, upper_bound.

See Also

mrct_region_sims for generating simulation results, summaryout_mrct for tabular summaries with trimming

Disjunctive (dummy) coding for factor columns

Description

Disjunctive (dummy) coding for factor columns

Usage

acm.disjctif(df)

Arguments

Value

Data frame with dummy-coded columns.

Add ID Column to Data Frame

Description

Ensures that a data frame has a unique ID column. If id.name is not provided, a column named "id" is added. If id.name is provided but does not exist in the data frame, it is created with unique integer values.

Usage

add_id_column(df.analysis, id.name = NULL)

Arguments

Value

Data frame with the ID column added if necessary.

Add Unprocessed Variables from Original Data

Description

Add Unprocessed Variables from Original Data

Usage

add_unprocessed_vars(
  df_work,
  data,
  outcome_var,
  event_var,
  treatment_var,
  continuous_vars,
  factor_vars,
  verbose
)

Analyze subgroup for summary table (OPTIMIZED)

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.

Apply Spline Constraint to Treatment Effect Coefficients

Description

Apply Spline Constraint to Treatment Effect Coefficients

Usage

apply_spline_constraint(b0, spline_var, knot, zeta, log_hrs, k_treat, verbose)

Assemble Final Results Object

Description

Assemble Final Results Object

Usage

assemble_results(
  df_super,
  mu,
  tau,
  gamma,
  b0,
  cens_model,
  subgroup_vars,
  subgroup_cuts,
  subgroup_definitions,
  hr_results,
  continuous_vars,
  factor_vars,
  model,
  n_super,
  seed,
  spline_info = NULL
)

Assign data to subgroups based on selected node

Description

Creates treatment recommendation flags based on identified subgroup

Usage

assign_subgroup_membership(data, best_subgroup, trees, X)

Arguments

Value

Data frame with added predict.node and treat.recommend columns

Bootstrap Results for ForestSearch with Bias Correction

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
)

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

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

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, showten_subgroups, 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

Bootstrap Ystar Matrix

Description

Generates a bootstrap matrix for Ystar using parallel processing.

Usage

bootstrap_ystar(df, nb_boots, seed = 8316951L)

Arguments

Value

Matrix of bootstrap samples (nb_boots x nrow(df)).

Build Classification Rate Table from Simulation Results

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
)

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

Build Cox Model Formula

Description

Constructs a Cox model formula from variable names.

Usage

build_cox_formula(outcome.name, event.name, treat.name)

Arguments

Value

An R formula object for Cox regression.

Build Estimation Properties Table from Simulation Results

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
)

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

Calculate Covariance for Bootstrap Estimates

Description

Calculates the covariance between a vector and bootstrap estimates.

Usage

calc_cov(x, Est)

Arguments

Value

Numeric value of covariance.

Calculate counts for subgroup summary

Description

Calculates sample size, treated count, and event count for a subgroup.

Usage

calculate_counts(Y, E, Treat, N)

Arguments

Value

List with formatted counts.

Calculate Event Counts by Treatment Arm

Description

Calculate Event Counts by Treatment Arm

Usage

calculate_event_counts(dd, tt, id.x)

Calculate Hazard Ratios from Potential Outcomes

Description

Calculate Hazard Ratios from Potential Outcomes

Usage

calculate_hazard_ratios(df_super, n_super, mu, tau, model, verbose)

Arguments

Value

List of hazard ratios

Calculate Linear Predictors for Potential Outcomes

Description

Calculate Linear Predictors for Potential Outcomes

Usage

calculate_linear_predictors(
  df_super,
  covariate_cols,
  gamma,
  b0,
  spline_info = NULL
)

Calculate Maximum Combinations

Description

Calculate Maximum Combinations

Usage

calculate_max_combinations(L, maxk)

Calculate potential outcome hazard ratio

Description

Calculates the average hazard ratio from a potential outcome variable.

Usage

calculate_potential_hr(df, potentialOutcome.name)

Arguments

Value

Numeric value of average hazard ratio.

Calculate Skewness

Description

Helper function to calculate sample skewness.

Usage

calculate_skewness(x)

Arguments

Value

Numeric skewness value

Calibrate Censoring Adjustment to Match DGM Reference Distribution

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,
  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

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

Calibrate k_inter for Target Subgroup Hazard Ratio

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,
  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

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

# Verify
dgm <- setup_gbsg_dgm(model = "alt", k_inter = k, verbose = FALSE)
print(dgm)

# Calibrate to AHR instead
k_ahr <- calibrate_k_inter(target_hr_harm = 1.5, use_ahr = TRUE, verbose = TRUE)
dgm_ahr <- setup_gbsg_dgm(model = "alt", k_inter = k_ahr, verbose = FALSE)
print(dgm_ahr)

Diagnose Censoring Consistency Between DGM Source Data and Simulated Data

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

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

Confidence Interval for Estimate

Description

Calculates confidence interval for an estimate, optionally on log(HR) scale.

Usage

ci_est(x, sd, alpha = 0.025, scale = "hr", est.loghr = TRUE)

Arguments

Value

List with length, lower, upper, sd, and estimate.

Compare Detection Curves Across Sample Sizes

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.

Compare Multiple Survival Regression Models

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

Compute AHR from loghr_po

Description

Computes Average Hazard Ratio from individual log hazard ratios.

Usage

compute_ahr(df, subset_indicator = NULL)

Arguments

Value

Numeric AHR value

Compute CDE from theta_0 and theta_1

Description

Computes Controlled Direct Effect as the ratio of average hazard contributions on the natural scale: CDE(S) = mean(exp(theta_1[S])) / mean(exp(theta_0[S])).

Usage

compute_cde(df, subset_indicator = NULL)

Arguments

Value

Numeric CDE value, or NA_real_ if columns are missing.

Compute Probability of Detecting True Subgroup

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 per treatment arm

• 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

# 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)")

Compute Detection Probability for Single Theta (Internal)

Description

Compute Detection Probability for Single Theta (Internal)

Usage

compute_detection_probability_single(
  theta,
  n_sg,
  prop_cens,
  k_avg,
  k_ind,
  method,
  n_mc,
  tol
)

Arguments

Value

Numeric probability

Compute and Attach CDE Values to a DGM Object

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

dgm <- setup_gbsg_dgm(model = "alt", k_inter = 2.0, verbose = FALSE)
dgm <- compute_dgm_cde(dgm)
dgm$hazard_ratios$CDE_harm   # theta-ddagger(H)
dgm$hazard_ratios$CDE         # theta-ddagger overall

Compute node metrics for a policy tree

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

Compute Hazard Ratio for a Single Subgroup

Description

Internal helper function to compute HR and CI for a subgroup. Uses robust (sandwich) standard errors for consistency with cox_summary().

Usage

compute_sg_hr(
  df,
  sg_name,
  outcome.name,
  event.name,
  treat.name,
  E.name,
  C.name,
  z_alpha = qnorm(0.975),
  conf.level = 0.95
)

Arguments

Value

Data frame with single row of HR estimates, or NULL if model fails.

Compute Hazard Ratio Estimates for Subgroups

Description

Internal function to compute Cox model hazard ratio estimates with confidence intervals for ITT, H, and Hc subgroups.

Usage

compute_sg_hr_estimates(
  df,
  df_H,
  df_Hc,
  outcome.name,
  event.name,
  treat.name,
  conf.level = 0.95,
  verbose = FALSE
)

Arguments

Value

Data frame with HR estimates

Compute Summary Statistics for Subgroups

Description

Internal function to compute summary statistics for each subgroup.

Usage

compute_sg_summary(
  df,
  df_H,
  df_Hc,
  outcome.name,
  event.name,
  treat.name,
  sg0_name,
  sg1_name
)

Arguments

Value

Data frame with summary statistics

Count ID Occurrences in Bootstrap Sample

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.

Comprehensive Wrapper for Cox Spline Analysis with AHR and CDE Plotting

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 = tempdir(),
  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

# Build a small synthetic dataset with required columns
set.seed(42)
n <- 200
df_ex <- data.frame(
  os_time   = rexp(n, rate = 0.01),
  os_event  = rbinom(n, 1, 0.6),
  treat     = rep(0:1, each = n / 2),
  biomarker = rnorm(n),
  loghr_po  = rnorm(n, mean = -0.3, sd = 0.5)
)

# With threshold - full subgroup analysis
results <- cox_ahr_cde_analysis(
  df = df_ex, z_name = "biomarker",
  hr_threshold = 1.25, plot_style = "grid",
  verbose = FALSE
)

# Without threshold - pure AHR curves
results <- cox_ahr_cde_analysis(
  df = df_ex, z_name = "biomarker",
  hr_threshold = NULL, plot_select = "ahr_only",
  verbose = FALSE
)

Fit Cox Model with Cubic Spline for Treatment Effect Heterogeneity

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

# 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
  )
)

Cox model summary for subgroup (OPTIMIZED)

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

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

Batch Cox summaries with caching

Description

For repeated calls with the same data structure but different subsets, this version pre-processes the data structure once.

Usage

cox_summary_batch(
  Y,
  E,
  Treat,
  Strata = NULL,
  subset_indices,
  return_format = c("formatted", "numeric")
)

Arguments

Value

List of results, one per subset.

Cox model summary for subgroup - vectorized version

Description

Efficiently processes multiple subgroups at once. Useful when analyzing many subgroups (e.g., in cross-validation).

Usage

cox_summary_vectorized(
  data,
  outcome_col,
  event_col,
  treat_col,
  strata_col = NULL,
  subgroup_col = "subgroup",
  return_format = c("formatted", "numeric")
)

Arguments

Value

Data frame with one row per subgroup and HR results.

Calculate Bootstrap Table Caption

Description

Generates an interpretive caption for bootstrap results table.

Usage

create_bootstrap_caption(est.scale, nb_boots, boot_success_rate)

Arguments

Value

Character string with caption

Create Bootstrap Diagnostic Plots

Description

Generates diagnostic visualization plots for bootstrap analysis.

Usage

create_bootstrap_diagnostic_plots(
  results,
  H_estimates,
  Hc_estimates,
  overall_timing = NULL
)

Arguments

Value

List of ggplot2 objects

Create Data Generating Mechanism for MRCT Simulations

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

Create Factor Summary Tables from Bootstrap Results

Description

Generates formatted GT tables summarizing factor frequencies from bootstrap subgroup analysis. Creates two complementary tables: one showing factor selection frequencies within each position (M.1, M.2, etc.), and another showing overall factor frequencies across all positions.

Usage

create_factor_summary_tables(factor_freq, n_found, min_percent = 2)

Arguments

Value

A list with up to two GT table objects:

by_position

GT table showing factor frequencies within each position. Percentages represent conditional probability of factor selection given that the position was populated. Within each position, percentages sum to approximately 100% (may not sum exactly to 100% after filtering).

overall

GT table showing total factor frequencies across all positions. Includes additional columns indicating which positions each factor appeared in and how many unique positions used the factor. Percentages represent proportion of successful iterations where the factor appeared in any position.

If no factors meet the minimum threshold, the corresponding table element will be NULL.

Note

This function requires the gt package for table creation. The overall table also requires dplyr for data aggregation. If dplyr is not available, only the position-specific table will be created and the overall element will be NULL.

Always check for NULL before using the returned tables:

if (!is.null(factor_tables$by_position)) {
  print(factor_tables$by_position)
}

If all factors have percentages below min_percent, both table elements will be NULL.

See Also

• summarize_bootstrap_subgroups for generating the factor_freq input

• format_subgroup_summary_tables for creating all subgroup summary tables

• summarize_bootstrap_results for complete bootstrap analysis workflow

• forestsearch_bootstrap_dofuture for running bootstrap analysis

Create Forest Plot Theme with Size Controls

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

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

# 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,
  ci_lwd = 2.5
)

Create Subgroup Indicator Columns from ForestSearch

Description

Internal helper to create Qrecommend and Brecommend indicator columns.

Usage

create_fs_subgroup_indicators(
  df,
  fs.est,
  col_names = c("Qrecommend", "Brecommend"),
  verbose = FALSE
)

Arguments

Value

Modified data frame with indicator columns.

Create GBSG-Based AFT Data Generating Mechanism

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

Helper Functions for GRF Subgroup Analysis

Description

This file contains helper functions used by grf.subg.harm.survival() to improve readability and modularity. Create GRF configuration object

Usage

create_grf_config(
  frac.tau,
  n.min,
  dmin.grf,
  RCT,
  sg.criterion,
  maxdepth,
  seedit
)

Arguments

Details

Creates a configuration object to organize GRF parameters

Value

List with configuration parameters

Create result object when no subgroup is found

Description

Builds result object for cases where no valid subgroup is identified

Usage

create_null_result(data, values, trees, config)

Arguments

Value

List with limited GRF results

Create Reference Subgroup Indicator Columns

Description

Creates indicator columns (0/1) in the data frame for each reference subgroup based on the provided subset expressions.

Usage

create_reference_subgroup_columns(df, ref_subgroups, verbose = FALSE)

Arguments

Value

List with modified df, cols, labels, and colors vectors.

Create Result Row

Description

Create Result Row

Usage

create_result_row(kk, covs.in, nx, event_counts, cox_result)

Create Sample Size Table for Multiple Scenarios

Description

Generates a table of required sample sizes for different combinations of true hazard ratios and censoring proportions.

Usage

create_sample_size_table(
  theta_values,
  prop_cens_values,
  target_power = 0.8,
  hr_threshold = 1.25,
  verbose = TRUE
)

Arguments

Value

A data.frame with columns: theta, prop_cens, n_required, achieved_power

Create Spline Variables

Description

Create Spline Variables

Usage

create_spline_variables(df_work, spline_var, knot)

Create Subgroup Indicator from Factor Definitions

Description

Parses factor definitions (e.g., "v1.1", "grade3.1") and creates a binary indicator for subgroup membership.

Usage

create_subgroup_indicator(df, sg_factors)

Arguments

Value

Integer vector (1 = in subgroup, 0 = not in subgroup)

Create Subgroup Summary Data Frame for Forest Plot

Description

Creates a data frame suitable for forestploter from multiple subgroup analyses. This is a more flexible alternative for complex subgroup configurations.

Usage

create_subgroup_summary_df(
  df_analysis,
  subgroups,
  outcome.name,
  event.name,
  treat.name,
  E.name = "E",
  C.name = "C",
  fs_bc_list = NULL,
  fs_kfold_list = NULL,
  conf.level = 0.95
)

Arguments

Value

Data frame with HR estimates for all subgroups.

Create result object for successful subgroup identification

Description

Builds comprehensive result object when a subgroup is found

Usage

create_success_result(
  data,
  best_subgroup,
  trees,
  tree_cuts,
  selected_tree,
  sg_harm_id,
  values,
  config
)

Arguments

Value

List with complete GRF results

Create Enhanced Summary Table for Baseline Characteristics

Description

Generates a formatted summary table comparing baseline characteristics between treatment arms. Supports continuous, categorical, and binary variables with p-values, standardized mean differences (SMD), and missing data summaries.

Usage

create_summary_table(
  data,
  treat_var = "treat",
  vars_continuous = NULL,
  vars_categorical = NULL,
  vars_binary = NULL,
  var_labels = NULL,
  digits = 1,
  show_pvalue = TRUE,
  show_smd = TRUE,
  show_missing = TRUE,
  table_title = "Baseline Characteristics by Treatment Arm",
  table_subtitle = NULL,
  source_note = NULL,
  font_size = 12,
  header_font_size = 14,
  footnote_font_size = 10,
  use_alternating_rows = TRUE,
  stripe_color = "#f9f9f9",
  indent_size = 20,
  highlight_pval = 0.05,
  highlight_smd = 0.2,
  highlight_color = "#fff3cd",
  compact_mode = FALSE,
  column_width_var = 200,
  column_width_stats = 120,
  show_column_borders = FALSE,
  custom_css = NULL
)

Arguments

Details

Binary variables specified via vars_binary display a single row showing the count and proportion for the "1" level. Categorical variables specified via vars_categorical that happen to be binary-coded (i.e., have exactly two levels: 0 and 1) are automatically detected and displayed in the same compact single-row format, showing only the "1" proportion.

Value

A gt table object (or data frame if gt not available)

Preset: Compact Table

Description

Preset: Compact Table

Usage

create_summary_table_compact(...)

Arguments

Preset: Minimal Table (No Highlighting, No Alternating)

Description

Preset: Minimal Table (No Highlighting, No Alternating)

Usage

create_summary_table_minimal(...)

Arguments

Preset: Presentation Table (Large Fonts)

Description

Preset: Presentation Table (Large Fonts)

Usage

create_summary_table_presentation(...)

Arguments

Preset: Publication-Ready Table

Description

Preset: Publication-Ready Table

Usage

create_summary_table_publication(...)

Arguments

Create Timing Summary Table

Description

Creates a data frame summarizing bootstrap timing information.

Usage

create_timing_summary_table(
  overall_timing,
  iteration_stats,
  fs_stats,
  overhead_stats,
  nb_boots,
  boot_success_rate
)

Arguments

Value

Data frame with timing summary

Discretize Continuous Variable into Quantile-Based Categories

Description

Discretize Continuous Variable into Quantile-Based Categories

Usage

cut_numeric(x, probs = c(0.25, 0.5, 0.75))

Arguments

Value

Integer vector with category codes (1 = lowest, max = highest)

Discretize Continuous Variable by Size Categories

Description

Discretize Continuous Variable by Size Categories

Usage

cut_size(x, breaks = c(20, 50))

Arguments

Value

Integer vector with category codes

Generate cut expressions for a variable

Description

For a continuous variable, returns expressions for mean, median, qlow, and qhigh cuts.

Usage

cut_var(x)

Arguments

Value

Character vector of cut expressions.

Compare Multiple CV Results

Description

Creates a comparison table from multiple cross-validation runs with different configurations.

Usage

cv_compare_results(
  cv_list,
  metrics = c("all", "finding", "agreement"),
  show_percentages = TRUE,
  digits = 1,
  use_gt = TRUE
)

Arguments

Value

A gt table or data.frame comparing CV results across configurations.

Create Metrics Tables for Cross-Validation Results

Description

Formats the find_summary and sens_summary outputs from forestsearch_tenfold or forestsearch_Kfold into publication-ready gt tables.

Usage

cv_metrics_tables(
  cv_result,
  sg_definition = NULL,
  title = "Cross-Validation Metrics",
  show_percentages = TRUE,
  digits = 1,
  include_raw = FALSE,
  table_style = c("combined", "separate", "minimal"),
  use_gt = TRUE
)

Arguments

Value

Depending on table_style:

• "combined": A single gt table (or data.frame)

• "separate": A list with agreement_table and finding_table

• "minimal": A single-row gt table (or data.frame)

If include_raw = TRUE, also includes sens_out and find_out matrices in the returned list.

See Also

cv_summary_tables for formatting forestsearch_KfoldOut(outall=TRUE) results

Create Summary Tables from forestsearch_KfoldOut Results

Description

Formats the detailed output from forestsearch_KfoldOut(outall=TRUE) into publication-ready gt tables. This includes ITT estimates, original subgroup estimates, and K-fold subgroup estimates.

Usage

cv_summary_tables(
  kfold_out,
  title = "Cross-Validation Summary",
  subtitle = NULL,
  show_metrics = TRUE,
  digits = 3,
  font_size = 12,
  use_gt = TRUE
)

Arguments

Value

If use_gt = TRUE, returns a list with gt table objects:

• combined_table: Combined ITT and subgroup estimates

• itt_table: ITT estimates only

• original_table: Original full-data subgroup estimates

• kfold_table: K-fold subgroup estimates

• metrics_table: Agreement and finding metrics (if show_metrics = TRUE)

If use_gt = FALSE, returns equivalent data.frames.

See Also

cv_metrics_tables for formatting forestsearch_tenfold() results

Create Compact CV Summary Text

Description

Generates a compact text string summarizing CV results, suitable for annotations in plots or reports.

Usage

cv_summary_text(
  cv_result,
  est.scale = "hr",
  include_finding = TRUE,
  include_agreement = TRUE
)

Arguments

Value

Character string with formatted CV metrics.

Default ForestSearch Parameters for GBSG Simulations

Description

Returns a list of default parameters for ForestSearch analysis in GBSG-based simulations.

Usage

default_fs_params()

Details

Default parameters are optimized for GBSG simulation scenarios with moderate sample sizes (300-1000) and typical event rates.

Variable selection defaults:

• use_lasso = TRUE: LASSO-based variable importance (default for FS)

• use_grf = FALSE: GRF-based variable importance (enable for FSlg)

The use_twostage parameter is set to FALSE by default for backward compatibility. Set to TRUE for faster exploratory analyses.

Value

List of default ForestSearch parameters

Default GRF Parameters for GBSG Simulations

Description

Returns a list of default parameters for GRF analysis in GBSG-based simulations. Parameters align with grf.subg.harm.survival() function signature.

Usage

default_grf_params()

Value

List of default GRF parameters

Default GRF parameters (general)

Description

Default GRF parameters (general)

Usage

default_grf_params_gen()

Default ForestSearch parameters (general)

Description

Default ForestSearch parameters (general)

Usage

default_sim_params()

Define Subgroups with Flexible Cutpoints

Description

Define Subgroups with Flexible Cutpoints

Usage

define_subgroups(
  df_work,
  data,
  subgroup_vars,
  subgroup_cuts,
  continuous_vars,
  model,
  verbose
)

Bivariate Density for Split-Sample HR Threshold Detection

Description

Computes the joint density for the two-split detection criterion where both split-halves must exceed individual thresholds and their average must exceed a consistency threshold.

Usage

density_threshold_both(x, theta, prop_cens = 0.3, n_sg, k_avg, k_ind)

Arguments

Details

The detection criterion requires:

• Average of two splits: (x1 + x2)/2 >= k_avg

• Individual splits: x1 >= k_ind AND x2 >= k_ind

Under the asymptotic approximation, each split-half log(HR) estimator follows N(log(theta), 8/d) where d = n_sg * (1 - prop_cens) / 2 is the expected number of events per split.

Value

Numeric. Joint density value at x, or 0 if thresholds not met.

Vectorized Density for Integration

Description

Wrapper around density_threshold_both for use with cubature integration.

Usage

density_threshold_integrand(x, theta, prop_cens, n_sg, k_avg, k_ind)

Arguments

Value

Numeric density value.

Automatically Detect Variable Types in a Dataset

Description

Analyzes a data frame to automatically classify variables as continuous or categorical, and returns a subset of the data with specified variables excluded.

Usage

detect_variable_types(data, max_unique_for_cat = 10, exclude_vars = NULL)

Arguments

Details

The function classifies variables using the following rules:

• Numeric variables with more than max_unique_for_cat unique values are classified as continuous

• Numeric variables with max_unique_for_cat or fewer unique values are classified as categorical

• Factor, character, and logical variables are always classified as categorical

• Variables listed in exclude_vars are omitted from classification and removed from the returned dataset

Value

A list containing:

Dummy-code a data frame (numeric pass-through, factors expanded)

Description

Dummy-code a data frame (numeric pass-through, factors expanded)

Usage

dummy_encode(df)

Arguments

Value

Data frame with numeric columns unchanged and factor columns expanded via acm.disjctif.

Early Stopping Decision

Description

Evaluates whether enough evidence exists to stop early based on confidence interval for consistency proportion.

Usage

early_stop_decision(
  n_success,
  n_total,
  threshold,
  conf.level = 0.95,
  min_samples = 20
)

Arguments

Value

Character. One of "continue", "pass", or "fail".

Evaluate a Single Factor Combination with Status Tracking

Description

Tests whether a specific combination meets all criteria and returns a status code indicating how far the evaluation progressed.

Usage

evaluate_combination_with_status(
  covs.in,
  yy,
  dd,
  tt,
  zz,
  n.min,
  d0.min,
  d1.min,
  hr.threshold,
  minp,
  rmin,
  kk
)

Arguments

Value

List with:

status

Integer status code: 0 = failed variance check, 1 = passed variance, failed prevalence, 2 = passed prevalence, failed redundancy, 3 = passed redundancy, failed events, 4 = passed events, failed sample size, 5 = passed sample size, failed Cox fit, 6 = passed Cox fit, failed HR threshold, 7 = passed all criteria (success)

result

Result row if successful, NULL otherwise

Evaluate a Comparison Expression Without eval(parse())

Description

Parses a string of the form "var op value" and evaluates it directly against a data frame column using operator dispatch. Falls back to column-name lookup for bare names.

Usage

evaluate_comparison(expr, df)

Arguments

Details

Supported operators (matched longest-first to avoid partial-match ambiguity): <=, >=, !=, ==, <, >.

If no operator is found, expr is treated as a column name and the result is df[[expr]] == 1.

The value on the right-hand side is coerced to numeric when possible, otherwise kept as character for string comparisons.

Value

Logical vector of length nrow(df).

Evaluate Consistency (Two-Stage Algorithm)

Description

Evaluates a single subgroup for consistency using a two-stage approach: Stage 1 screens with fewer splits, Stage 2 uses sequential batched evaluation with early stopping for efficient evaluation.

Usage

evaluate_consistency_twostage(
  m,
  index.Z,
  names.Z,
  df,
  found.hrs,
  hr.consistency,
  pconsistency.threshold,
  pconsistency.digits = 2,
  maxk,
  confs_labels,
  details = FALSE,
  n.splits.screen = 30,
  screen.threshold = NULL,
  n.splits.max = 400,
  batch.size = 20,
  conf.level = 0.95,
  min.valid.screen = 10
)

Arguments

Value

Named numeric vector with consistency results, or NULL if not met.

Cache and validate cut expressions efficiently

Description

Evaluates all cut expressions once and caches results to avoid redundant evaluation. Much faster than evaluating repeatedly.

Usage

evaluate_cuts_once(confs, df, details = FALSE)

Arguments

Details

This replaces multiple eval(parse()) calls scattered throughout get_FSdata. By caching results, we avoid:

1. Repeated parsing of expressions

2. Repeated evaluation on dataframe

3. Redundant uniqueness checks

Value

List with:

• evaluations: List of evaluated vectors (logical TRUE/FALSE) for each cut

• is_valid: Logical vector indicating which cuts produced >1 unique value

• has_error: Logical vector indicating which cuts failed to evaluate

Evaluate Single Subgroup for Consistency (Fixed-Sample)

Description

Evaluates a single subgroup for consistency across random splits using a fixed number of splits.

Usage

evaluate_subgroup_consistency(
  m,
  index.Z,
  names.Z,
  df,
  found.hrs,
  n.splits,
  hr.consistency,
  pconsistency.threshold,
  pconsistency.digits = 2,
  maxk,
  confs_labels,
  details = FALSE
)

Arguments

Value

Named numeric vector with consistency results, or NULL if criteria not met.

Extract all cuts from fitted trees

Description

Consolidates cut information from all fitted policy trees. This is the default behavior that returns cuts from all trees regardless of which tree identified the selected subgroup.

Usage

extract_all_tree_cuts(trees, maxdepth)

Arguments

Value

List with cuts and names for each tree and combined

Extract Estimates from ForestSearch Results

Description

Extracts operating characteristics (HRs, classification metrics, etc.) from ForestSearch analysis results. Aligned with new DGM output structure.

Usage

extract_fs_estimates(
  df,
  fs_res,
  dgm,
  cox_formula = NULL,
  cox_formula_adj = NULL,
  analysis = "FS",
  fs_full = NULL,
  verbose = FALSE
)

Arguments

Value

data.table with extracted estimates including AHR metrics

Extract Subgroup Definition from ForestSearch Object

Description

Internal helper to extract human-readable subgroup definition.

Usage

extract_fs_subgroup_definition(fs.est, verbose = FALSE)

Arguments

Value

Character string describing the subgroup definition.

Extract Estimates from GRF Results

Description

Aligned with new DGM output structure including AHR metrics. Correctly handles grf.subg.harm.survival() output structure:

• sg.harm.id: subgroup definition string

• data: data frame with treat.recommend column (0 = harm, 1 = complement)

Usage

extract_grf_estimates(
  df,
  grf_est,
  dgm,
  cox_formula = NULL,
  cox_formula_adj = NULL,
  analysis = "GRF",
  frac_tau = 1,
  verbose = FALSE,
  debug = FALSE
)

Arguments

Value

data.table with extracted estimates

Extract redundancy flag for subgroup combinations

Description

Checks if adding each factor to a subgroup reduces the sample size by at least rmin.

Usage

extract_idx_flagredundancy(x, rmin)

Arguments

Value

List with id.x (membership vector) and flag.redundant (logical).

Extract cuts from selected tree only

Description

Extracts cut information only from the tree at the specified selected depth. This provides a focused set of cuts from the tree that identified the subgroup meeting the dmin.grf criterion, rather than cuts from all trees.

Usage

extract_selected_tree_cuts(trees, selected_depth, maxdepth)

Arguments

Details

This function is used when return_selected_cuts_only = TRUE in grf.subg.harm.survival(). It returns:

• tree1, tree2, tree3: Individual tree cuts (still populated for reference)

• names1, names2, names3: Individual tree variable names

• all: Cuts from the SELECTED tree only (not union of all trees)

• all_names: Variable names from the SELECTED tree only

• selected_depth: The depth that was selected

Value

List with cuts and names, structured similarly to extract_all_tree_cuts but with only the selected tree's cuts in the all field

Extract Subgroup Information

Description

Extracts subgroup definition and membership from results.

Usage

extract_subgroup(df, top_result, index.Z, names.Z, confs_labels)

Arguments

Value

List with sg.harm, sg.harm_label, df_flag, sg.harm.id.

Extract cut information from a policy tree

Description

Extracts all split points and variables from a policy tree

Usage

extract_tree_cuts(tree)

Arguments

Value

List with cuts (expressions) and names (unique variables)

Generate Figure Note for Quarto/RMarkdown

Description

Formats the figure note from plot_sg_weighted_km() output for use in Quarto or RMarkdown documents.

Usage

figure_note(
  x,
  prefix = "*Note*: ",
  include_definition = TRUE,
  include_hr_explanation = TRUE,
  custom_text = NULL
)

Arguments

Value

Character string formatted as a figure note, or NULL if no content

Filter a vector by LASSO-selected variables

Description

Returns elements of x that are in lassokeep.

Usage

filter_by_lassokeep(x, lassokeep)

Arguments

Value

Filtered character vector or NULL.

Filter and merge arguments for function calls

Description

Simplifies the common pattern of filtering arguments from a source list to match a target function's formal parameters, then adding/overriding specific arguments.

Usage

filter_call_args(source_args, target_func, override_args = NULL)

Arguments

Details

This function:

1. Extracts formal parameter names from target_func

2. Keeps only arguments from source_args that match those names

3. Adds or overrides with any override_args provided

Reduces boilerplate and improves readability across the codebase.

Value

List of filtered arguments ready for do.call().

Find Covariate Any Match

Description

Helper function to determine if any CV fold found a subgroup involving the same covariate (not necessarily same cut).

Usage

find_covariate_any_match(sg_target, sg1, sg2, confs)

Arguments

Value

Numeric vector (0/1) indicating match for each fold.

Find k_inter Value to Achieve Target Harm Subgroup Hazard Ratio

Description

Uses numerical root-finding to determine the interaction parameter (k_inter) that achieves a specified target hazard ratio in the harm subgroup. This is the most efficient method for single target calibration.

Usage

find_k_inter_for_target_hr(
  target_hr_harm,
  data,
  continuous_vars,
  factor_vars,
  outcome_var,
  event_var,
  treatment_var,
  subgroup_vars,
  subgroup_cuts,
  k_treat = 1,
  k_inter_range = c(-10, 10),
  tol = 0.001,
  n_super = 5000,
  verbose = TRUE
)

Arguments

Details

This function uses the uniroot algorithm to solve the equation:

HR_{harm}(k_{inter}) - HR_{target} = 0

The algorithm typically converges within 5-10 iterations and achieves high precision (within the specified tolerance). If the root-finding fails, the function evaluates the boundaries and provides diagnostic information.

Value

A list of class "k_inter_result" containing:

k_inter

Numeric value of optimal k_inter parameter

achieved_hr_harm

Numeric value of achieved hazard ratio in harm subgroup

target_hr_harm

Numeric value of target hazard ratio (for reference)

error

Numeric value of absolute error between achieved and target HR

dgm

Object of class "aft_dgm_flex" containing the final DGM

convergence

Integer number of iterations to convergence

method

Character string "root-finding" indicating method used

See Also

sensitivity_analysis_k_inter for sensitivity analysis generate_aft_dgm_flex for DGM generation

Find the split that leads to a specific leaf node

Description

Identifies the split point that creates a given leaf node

Usage

find_leaf_split(tree, leaf_node)

Arguments

Value

Character string with split expression or NULL

Find Quantile for Target Subgroup Proportion

Description

Determines the quantile cutpoint that achieves a target proportion of observations in a subgroup. Useful for calibrating subgroup sizes.

Usage

find_quantile_for_proportion(
  data,
  var_name,
  target_prop,
  direction = "less",
  tol = 1e-04
)

Arguments

Details

This function uses root finding (uniroot) to determine the quantile that results in exactly the target proportion of observations being classified into the subgroup. This is particularly useful when you want to ensure a specific subgroup size regardless of the data distribution.

Value

A list containing:

quantile

The quantile value (between 0 and 1) that achieves the target proportion

cutpoint

The actual data value corresponding to this quantile

actual_proportion

The achieved proportion (should equal target_prop within tolerance)

See Also

generate_aft_dgm_flex

Find Minimum Sample Size for Target Detection Power

Description

Determines the minimum subgroup sample size needed to achieve a target detection probability for a given true hazard ratio.

Usage

find_required_sample_size(
  theta,
  target_power = 0.8,
  prop_cens = 0.3,
  hr_threshold = 1.25,
  hr_consistency = 1,
  n_range = c(20L, 500L),
  tol = 1,
  verbose = TRUE
)

Arguments

Value

A list with:

Fit AFT Model with Optional Spline Treatment Effect

Description

Fit AFT Model with Optional Spline Treatment Effect

Usage

fit_aft_model(
  df_work,
  interaction_term,
  k_treat,
  k_inter,
  verbose,
  spline_spec = NULL,
  set_var = NULL,
  beta_var = NULL
)

Fit AFT Model with Spline Treatment Effect

Description

Fit AFT Model with Spline Treatment Effect

Usage

fit_aft_model_spline(
  df_work,
  covariate_cols,
  interaction_term,
  k_treat,
  k_inter,
  spline_spec,
  verbose,
  set_var = NULL,
  beta_var = NULL
)

Fit Standard AFT Model (Non-Spline)

Description

Fit Standard AFT Model (Non-Spline)

Usage

fit_aft_model_standard(
  df_work,
  covariate_cols,
  interaction_term,
  k_treat,
  k_inter,
  verbose,
  set_var = NULL,
  beta_var = NULL
)

Fit causal survival forest

Description

Wrapper function to fit GRF causal survival forest with appropriate settings

Usage

fit_causal_forest(X, Y, W, D, tau.rmst, RCT, seedit)

Arguments

Value

Causal survival forest object

Fit Cox Model for Subgroup

Description

Fit Cox Model for Subgroup

Usage

fit_cox_for_subgroup(yy, dd, tt, id.x)

Fit Cox Models for Subgroups

Description

Fits Cox models for two subgroups defined by treatment recommendation.

Usage

fit_cox_models(df, formula)

Arguments

Value

List with HR and SE for each subgroup.

Fit policy trees up to specified depth

Description

Fits policy trees of depths 1 through maxdepth and computes metrics

Usage

fit_policy_trees(X, data, dr.scores, maxdepth, n.min)

Arguments

Value

List with trees and combined values

ForestSearch: Exploratory Subgroup Identification

Description

Identifies subgroups with differential treatment effects in clinical trials using a combination of Generalized Random Forests (GRF), LASSO variable selection, and exhaustive combinatorial search with split-sample validation.

Usage

forestsearch(
  df.analysis,
  outcome.name = "tte",
  event.name = "event",
  treat.name = "treat",
  id.name = "id",
  potentialOutcome.name = NULL,
  flag_harm.name = NULL,
  confounders.name = NULL,
  parallel_args = list(plan = "callr", workers = 6, show_message = TRUE),
  df.predict = NULL,
  df.test = NULL,
  is.RCT = TRUE,
  seedit = 8316951,
  est.scale = "hr",
  use_lasso = TRUE,
  use_grf = TRUE,
  grf_res = NULL,
  grf_cuts = NULL,
  max_n_confounders = 1000,
  grf_depth = 2,
  dmin.grf = 12,
  frac.tau = 0.6,
  return_selected_cuts_only = TRUE,
  conf_force = NULL,
  defaultcut_names = NULL,
  cut_type = "default",
  exclude_cuts = NULL,
  replace_med_grf = FALSE,
  cont.cutoff = 4,
  conf.cont_medians = NULL,
  conf.cont_medians_force = NULL,
  n.min = 60,
  hr.threshold = 1.25,
  hr.consistency = 1,
  sg_focus = "hr",
  fs.splits = 1000,
  m1.threshold = Inf,
  pconsistency.threshold = 0.9,
  stop_threshold = 0.95,
  showten_subgroups = FALSE,
  d0.min = 12,
  d1.min = 12,
  max.minutes = 3,
  minp = 0.025,
  details = FALSE,
  maxk = 2,
  by.risk = 12,
  plot.sg = FALSE,
  plot.grf = FALSE,
  max_subgroups_search = 10,
  vi.grf.min = -0.2,
  use_twostage = TRUE,
  twostage_args = list()
)

Arguments

Details

Algorithm Overview:

1. Variable Selection: GRF identifies variables with treatment effect heterogeneity; LASSO selects most predictive

2. Subgroup Discovery: Exhaustive search over factor combinations up to maxk

3. Consistency Validation: Split-sample validation ensures reproducibility

4. Selection: Choose subgroup based on sg_focus criterion

Two-Stage Consistency Algorithm: When use_twostage = TRUE, the consistency evaluation uses an optimized algorithm that can provide 3-10x speedup:

• Stage 1: Quick screening with n.splits.screen splits eliminates clearly non-viable candidates

• Stage 2: Sequential batched evaluation with early stopping for candidates passing Stage 1

The two-stage algorithm is recommended for:

• Exploratory analyses with many candidate subgroups

• Large fs.splits values (>200)

• Iterative model development

For final regulatory submissions, use_twostage = FALSE may be preferred for exact reproducibility.

Value

A list of class "forestsearch" containing:

grp.consistency

Consistency evaluation results including:

• out_sg: Selected subgroup based on sg_focus

• sg_focus: Focus criterion used

• df_flag: Treatment recommendations

• algorithm: "twostage" or "fixed"

• n_candidates_evaluated: Number evaluated

• n_passed: Number passing threshold

find.grps

Subgroup search results

confounders.candidate

Candidate confounders considered

confounders.evaluated

Confounders after variable selection

df.est

Analysis data with treatment recommendations

df.predict

Prediction data with recommendations (if provided)

df.test

Test data with recommendations (if provided)

minutes_all

Total computation time

grf_res

GRF results object

sg_focus

Subgroup focus criterion used

sg.harm

Selected subgroup definition

grf_cuts

GRF cut points used

prop_maxk

Proportion of max combinations searched

max_sg_est

Maximum subgroup HR estimate

grf_plot

GRF plot object (if plot.grf = TRUE)

args_call_all

All arguments for reproducibility

References

• FDA Guidance for Industry: Enrichment Strategies for Clinical Trials

• Athey & Imbens (2016). Recursive partitioning for heterogeneous causal effects. PNAS.

• Wager & Athey (2018). Estimation and inference of heterogeneous treatment effects using random forests. JASA.

See Also

subgroup.consistency for consistency evaluation details forestsearch_bootstrap_dofuture for bootstrap inference forestsearch_Kfold for cross-validation

Package website: https://larry-leon.github.io/forestsearch/

Source code: https://github.com/larry-leon/forestsearch

ForestSearch K-Fold Cross-Validation

Description

This function assesses the stability and reproducibility of ForestSearch subgroup identification through cross-validation. For each fold:

1. Train ForestSearch on (K-1) folds

2. Apply the identified subgroup to the held-out fold

3. Compare predictions to the original full-data analysis

Usage

forestsearch_Kfold(
  fs.est,
  Kfolds = nrow(fs.est$df.est),
  seedit = 8316951L,
  parallel_args = list(plan = "multisession", workers = 6, show_message = TRUE),
  sg0.name = "Not recommend",
  sg1.name = "Recommend",
  details = FALSE
)

Arguments

Details

Performs K-fold cross-validation for ForestSearch, evaluating subgroup identification and agreement between training and test sets.

Value

List with components:

resCV

Data frame with CV predictions for each observation

cv_args

Arguments used for CV ForestSearch calls

timing_minutes

Execution time in minutes

prop_SG_found

Percentage of folds where a subgroup was found

sg_analysis

Original subgroup definition from full-data analysis

sg0.name, sg1.name

Subgroup labels

Kfolds

Number of folds used

sens_summary

Named vector of sensitivity metrics (sens_H, sens_Hc, ppv_H, ppv_Hc)

find_summary

Named vector of subgroup-finding metrics (Any, Exact, etc.)

Cross-Validation Types

• Leave-One-Out (LOO): When Kfolds = nrow(df), each observation is held out once. Most thorough but computationally intensive.

• K-Fold: When Kfolds < nrow(df), data is split into K roughly equal folds. Good balance of bias-variance tradeoff.

Output Metrics

The returned resCV data frame contains:

• treat.recommend: Prediction from CV model

• treat.recommend.original: Prediction from full-data model

• cvindex: Fold assignment

• sg1, sg2: Subgroup definitions found in each fold

See Also

forestsearch for initial subgroup identification forestsearch_KfoldOut for summarizing CV results forestsearch_tenfold for repeated K-fold simulations

ForestSearch K-Fold Cross-Validation Output Summary

Description

Summarizes cross-validation results for ForestSearch, including subgroup agreement and performance metrics.

Usage

forestsearch_KfoldOut(res, details = FALSE, outall = FALSE)

Arguments

Value

If outall=FALSE, a list with sens_metrics_original and find_metrics. If outall=TRUE, a list with summary tables and metrics.

ForestSearch Bootstrap with doFuture Parallelization

Description

Orchestrates bootstrap analysis for ForestSearch using doFuture parallelization. Implements bias correction methods to adjust for optimism in subgroup selection.

Usage

forestsearch_bootstrap_dofuture(
  fs.est,
  nb_boots,
  seed = 8316951L,
  details = FALSE,
  show_three = FALSE,
  parallel_args = list()
)

Arguments

Value

List with the following components:

results

Data.table with bias-corrected estimates for each bootstrap iteration

SG_CIs

List of confidence intervals for H and Hc (raw and bias-corrected)

FSsg_tab

Formatted table of subgroup estimates

Ystar_mat

Matrix (nb_boots x n) of bootstrap sample indicators

H_estimates

Detailed estimates for subgroup H

Hc_estimates

Detailed estimates for subgroup Hc

summary

(If create_summary=TRUE) Enhanced summary with tables and diagnostics

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.

Variable Naming Convention

• H: Original subgroup (harm/questionable, treat.recommend == 0)

• Hc: Complement subgroup (recommend, treat.recommend == 1)

• _obs: Estimate from original data

• _star: Estimate from bootstrap data

• _biasadj_1: Bias correction method 1

• _biasadj_2: Bias correction method 2

Performance

Typical runtime: 1-5 seconds per bootstrap iteration. For 1000 bootstraps with 6 workers, expect 3-10 minutes total. Memory usage scales with dataset size and number of workers.

Requirements

• Original fs.est must have identified a valid subgroup

• Requires packages: data.table, foreach, doFuture, survival

• For plots: requires ggplot2

See Also

forestsearch for initial subgroup identification bootstrap_results for the core bootstrap worker function build_cox_formula for Cox formula construction fit_cox_models for Cox model fitting

ForestSearch Repeated K-Fold Cross-Validation

Description

This function performs multiple independent K-fold cross-validations to assess the variability in subgroup identification. Each simulation:

1. Randomly shuffles the data

2. Performs K-fold CV

3. Records sensitivity and agreement metrics

Results are summarized across all simulations.

Usage

forestsearch_tenfold(
  fs.est,
  sims,
  Kfolds = 10,
  details = TRUE,
  seed = 8316951L,
  parallel_args = list(plan = "multisession", workers = 6, show_message = TRUE)
)

Arguments

Details

Runs repeated K-fold cross-validation simulations for ForestSearch and summarizes subgroup identification stability across repetitions.

Value

List with components:

sens_summary

Named vector of median sensitivity metrics across simulations

find_summary

Named vector of median subgroup-finding metrics

sens_out

Matrix of sensitivity metrics (sims x metrics)

find_out

Matrix of finding metrics (sims x metrics)

timing_minutes

Total execution time

sims

Number of simulations run

Kfolds

Number of folds per simulation

Parallelization Strategy

Unlike the single K-fold function which parallelizes across folds, this function parallelizes across simulations for better efficiency when running many repetitions. Each simulation runs its K-fold CV sequentially.

See Also

forestsearch_Kfold for single K-fold CV forestsearch_KfoldOut for summarizing CV results

Format Confidence Interval for Estimates

Description

Formats confidence interval for estimates.

Usage

format_CI(estimates, col_names)

Arguments

Value

Character string formatted as \"estimate (lower, upper)\".

Format Bootstrap Diagnostics Table with gt

Description

Creates a publication-ready diagnostics table from bootstrap results.

Usage

format_bootstrap_diagnostics_table(
  diagnostics,
  nb_boots,
  results,
  H_estimates = NULL,
  Hc_estimates = NULL
)

Arguments

Value

A gt table object

Format Bootstrap Results Table with gt

Description

Creates a publication-ready table from ForestSearch bootstrap results, with bias-corrected confidence intervals, informative formatting, and optional subgroup definition footnote.

Usage

format_bootstrap_table(
  FSsg_tab,
  nb_boots,
  est.scale = "hr",
  boot_success_rate = NULL,
  sg_definition = NULL,
  title = NULL,
  subtitle = NULL
)

Arguments

Value

A gt table object

Format Bootstrap Timing Table with gt

Description

Creates a publication-ready timing summary table from bootstrap results.

Usage

format_bootstrap_timing_table(timing_list, nb_boots, boot_success_rate)

Arguments

Value

A gt table object

Format Continuous Variable Definition for Display

Description

Format Continuous Variable Definition for Display

Usage

format_continuous_definition(var_data, cut_spec, var_name)

Format ForestSearch Details Output for Beamer Two-Column Display

Description

Captures forestsearch(details = TRUE) console output and splits it into two columns for readable beamer slides. Left column shows variable selection (GRF, LASSO, candidate factors); right column shows subgroup search, consistency evaluation, and results.

Usage

format_fs_details(
  fs_output,
  split_after = "Candidate factors",
  fontsize = "scriptsize",
  col_widths = c(0.48, 0.52),
  max_width = 48
)

Arguments

Value

Invisibly returns a list with left and right character vectors. Side effect: emits LaTeX via cat() for use in a chunk with results='asis'.

Quarto Setup

No special LaTeX packages required. Works in any beamer frame without the fragile option.

Usage

In a Quarto beamer chunk with results='asis' and echo=FALSE, first capture the forestsearch output with capture.output(), then call format_fs_details(fs_output).

Format Operating Characteristics Results as GT Table

Description

Creates a formatted gt table from simulation operating characteristics results.

Usage

format_oc_results(
  results,
  analyses = NULL,
  metrics = "all",
  digits = 3,
  digits_hr = 3,
  title = "Operating Characteristics Summary",
  subtitle = NULL,
  use_gt = TRUE
)

Arguments

Details

The function summarizes simulation results across multiple metrics:

• Found: Proportion of simulations finding a subgroup (any.H)

• Classification: Sen, spec, PPV, NPV

• HR Estimates: Mean Cox hazard ratios in true (H) and identified (H-hat) subgroups and their complements

• AHR Estimates: Mean average hazard ratios (from loghr_po) in true and identified subgroups

• CDE Estimates: Controlled direct effects (from theta_0/theta_1) in true and identified subgroups

• Subgroup Size: Average, min, max sizes

Column notation aligns with build_estimation_table and Leon et al. (2024): H = true (oracle) subgroup, H-hat = identified subgroup. The asterisk (*) is reserved for bootstrap bias-corrected estimates and is not used in this summary table.

Value

A gt table object (if use_gt = TRUE and gt package available) or data.frame

Format results for subgroup summary

Description

Formats results for subgroup summary table.

Usage

format_results(
  subgroup_name,
  n,
  n_treat,
  d,
  m1,
  m0,
  drmst,
  hr,
  hr_a = NA,
  hr_po = NA,
  return_medians = TRUE
)

Arguments

Value

Character vector of results.

Format Search Results

Description

Format Search Results

Usage

format_search_results(
  results_list,
  Z,
  details,
  t.sofar,
  L,
  max_count,
  filter_counts = NULL
)

Arguments

Format Subgroup Summary Tables with gt

Description

Creates publication-ready gt tables for bootstrap subgroup analysis

Usage

format_subgroup_summary_tables(subgroup_summary, nb_boots)

Arguments

Value

List of gt table objects

Generate Synthetic Survival Data using AFT Model with Flexible Subgroups

Description

Creates a data generating mechanism (DGM) for survival data using an Accelerated Failure Time (AFT) model with Weibull distribution. Supports flexible subgroup definitions and treatment-subgroup interactions.

Usage

generate_aft_dgm_flex(
  data,
  continuous_vars,
  factor_vars,
  continuous_vars_cens = NULL,
  factor_vars_cens = NULL,
  set_beta_spec = list(set_var = NULL, beta_var = NULL),
  outcome_var,
  event_var,
  treatment_var = NULL,
  subgroup_vars = NULL,
  subgroup_cuts = NULL,
  draw_treatment = FALSE,
  model = "alt",
  k_treat = 1,
  k_inter = 1,
  n_super = 5000,
  select_censoring = TRUE,
  cens_type = "weibull",
  cens_params = list(),
  seed = 8316951,
  verbose = TRUE,
  standardize = FALSE,
  spline_spec = NULL
)

Arguments

Details

Subgroup Cutpoint Specifications

The subgroup_cuts parameter accepts multiple flexible specifications:

Fixed Value

subgroup_cuts = list(er = 20)  # er <= 20

Quantile-based

subgroup_cuts = list(
  er = list(type = "quantile", value = 0.25)  # er <= 25th percentile
)

Function-based

subgroup_cuts = list(
  er = list(type = "function", fun = median)  # er <= median
)

Range

subgroup_cuts = list(
  age = list(type = "range", min = 40, max = 60)  # 40 <= age <= 60
)

Greater than

subgroup_cuts = list(
  nodes = list(type = "greater", quantile = 0.75)  # nodes > 75th percentile
)

Multiple values (for categorical)

subgroup_cuts = list(
  grade = list(type = "multiple", values = c(2, 3))  # grade in (2, 3)
)

Custom function

subgroup_cuts = list(
  er = list(
    type = "custom",
    fun = function(x) x <= quantile(x, 0.3) | x >= quantile(x, 0.9)
  )
)

Model Structure

The AFT model with Weibull distribution is specified as:

\log(T) = \mu + \gamma' X + \sigma \epsilon

Where:

• T is the survival time

• \mu is the intercept

• \gamma contains the covariate effects

• X includes treatment, covariates, and treatment x subgroup interaction

• \sigma is the scale parameter

• \epsilon follows an extreme value distribution

Interaction Term

The model creates a SINGLE interaction term representing the treatment effect modification when ALL subgroup conditions are simultaneously satisfied. This is not multiple separate interactions but one combined indicator.

Value

A named list of class aft_dgm containing:

Author(s)

Your Name

References

Leon, L.F., et al. (2024). Statistics in Medicine.

Kalbfleisch, J.D. and Prentice, R.L. (2002). The Statistical Analysis of Failure Time Data (2nd ed.). Wiley.

Examples

df <- survival::gbsg
dgm <- generate_aft_dgm_flex(
  data            = df,
  outcome_var     = "rfstime",
  event_var       = "status",
  treatment_var   = "hormon",
  continuous_vars = c("age", "size", "nodes", "pgr", "er"),
  factor_vars     = "meno",
  model           = "null",
  verbose         = FALSE
)
str(dgm)

Generate Synthetic Data using Bootstrap with Perturbation

Description

Generate Synthetic Data using Bootstrap with Perturbation

Usage

generate_bootstrap_synthetic(
  data,
  continuous_vars,
  cat_vars,
  n = NULL,
  seed = 123,
  noise_level = 0.1,
  id_var = NULL,
  cat_flip_prob = NULL,
  preserve_bounds = TRUE,
  ordinal_vars = NULL
)

Arguments

Value

A data frame with synthetic data

Generate Bootstrap Sample with Added Noise

Description

Creates a bootstrap sample from a dataset with controlled noise added to both continuous and categorical variables. This function is useful for generating synthetic datasets that maintain the general structure of the original data while introducing controlled variation.

Usage

generate_bootstrap_with_noise(
  data,
  n = NULL,
  continuous_vars = NULL,
  cat_vars = NULL,
  id_var = "pid",
  seed = 123,
  noise_level = 0.1
)

Arguments

Details

The function performs the following operations:

Bootstrap Sampling

Samples n observations with replacement from the original dataset.

Continuous Variable Noise

• Adds Gaussian noise with standard deviation = original SD × noise_level

• Constrains values to remain within original variable bounds

• Preserves integer type for variables that appear to be integers

Categorical Variable Perturbation

• Changes values with probability = noise_level / 2

• Binary variables: flips to opposite value

• Multi-level unordered: randomly selects from other levels

• Ordered factors: weights selection toward adjacent levels

• Preserves factor levels and ordering from original data

Value

A data frame with the same structure as the input data, containing bootstrap sampled observations with added noise.

Note

• The function assumes that categorical variables with numeric encoding should maintain their numeric type unless they are factors in the input

• Missing values (NA) are handled appropriately in calculations but are not imputed

• For ordered factors or variables named "grade", the perturbation favors transitions to adjacent levels over distant levels

See Also

sample for bootstrap sampling, rnorm for noise generation

Generate Combination Indices

Description

Creates indices for all factor combinations up to maxk

Usage

generate_combination_indices(L, maxk)

Generate Complement Expression

Description

Creates the logical complement of a subgroup expression. Handles common patterns like "var <= x" -> "var > x".

Usage

generate_complement_expression(expr)

Arguments

Value

Character string with negated expression.

Generate Detection Probability Curve

Description

Computes detection probability across a range of hazard ratios to create a power-like curve for subgroup detection.

Usage

generate_detection_curve(
  theta_range = c(0.5, 3),
  n_points = 50L,
  n_sg,
  prop_cens = 0.3,
  hr_threshold = 1.25,
  hr_consistency = 1,
  include_reference = TRUE,
  method = "cubature",
  verbose = TRUE
)

Arguments

Value

A data.frame with columns:

Generate Synthetic GBSG Data using Generalized Bootstrap

Description

Generate Synthetic GBSG Data using Generalized Bootstrap

Usage

generate_gbsg_bootstrap_general(n = 686, seed = 123, noise_level = 0.1)

Arguments

Value

Synthetic GBSG dataset

Generate Readable Subgroup Labels from ForestSearch Object

Description

Extracts human-readable subgroup labels that are also valid R expressions for use with plotKM.band_subgroups(). Attempts to extract the actual subgroup definition (e.g., "er <= 0") rather than column references.

Usage

generate_readable_sg_labels(fs.est, verbose = FALSE)

Arguments

Value

Character vector of length 2: c(harm_label, benefit_label)

Generate Super Population and Calculate Linear Predictors

Description

Generate Super Population and Calculate Linear Predictors

Usage

generate_super_population(
  df_work,
  n_super,
  draw_treatment,
  gamma,
  b0,
  mu,
  tau,
  verbose,
  spline_info = NULL
)

Fit Cox Model for Subgroup

Description

Fits a Cox model for a subgroup and returns estimate and standard error.

Usage

get_Cox_sg(df_sg, cox.formula, est.loghr = TRUE, cox_initial = log(1))

Arguments

Details

Function is utilized throughout codebase

Value

List with estimate and standard error.

ForestSearch Data Preparation and Feature Selection

Description

Prepares a dataset for ForestSearch, including options for LASSO-based dimension reduction, GRF cuts, forced cuts, and flexible cut strategies. Returns a list with the processed data, subgroup factor names, cut expressions, and LASSO selection results.

Usage

get_FSdata(
  df.analysis,
  use_lasso = FALSE,
  use_grf = FALSE,
  grf_cuts = NULL,
  confounders.name,
  cont.cutoff = 4,
  conf_force = NULL,
  conf.cont_medians = NULL,
  conf.cont_medians_force = NULL,
  replace_med_grf = TRUE,
  defaultcut_names = NULL,
  cut_type = "default",
  exclude_cuts = NULL,
  outcome.name = "tte",
  event.name = "event",
  details = TRUE
)

Arguments

Value

A named list containing:

Get Best Model from Comparison

Description

Extracts the best fitting model object from a comparison result. If no single best model can be determined, returns the Weibull model if selected by either AIC or BIC. Defaults to Weibull0 model if no model can be determined.

Usage

get_best_survreg(comparison_result)

Arguments

Value

A survreg model object (defaults to Weibull0 model)

Get all exported functions from ForestSearch namespace

Description

Get all exported functions from ForestSearch namespace

Usage

get_bootstrap_exports()

Get all combinations of subgroup factors up to maxk

Description

Generates all possible combinations of subgroup factors up to a specified maximum size.

Usage

get_combinations_info(L, maxk)

Arguments

Value

List with max_count (total combinations) and indices_list (indices for each k).

Get forced cut expressions for variables

Description

For each variable in conf.force.names, returns cut expressions if continuous.

Usage

get_conf_force(df, conf.force.names, cont.cutoff = 4)

Arguments

Value

Character vector of cut expressions.

Get indicator vector for selected subgroup factors

Description

Returns a vector indicating which factors are included in a subgroup combination.

Usage

get_covs_in(
  kk,
  maxk,
  L,
  counts_1factor,
  index_1factor,
  counts_2factor = NULL,
  index_2factor = NULL,
  counts_3factor = NULL,
  index_3factor = NULL
)

Arguments

Value

Numeric vector indicating selected factors (1 = included, 0 = not included).

Get variable name from cut expression

Description

Extracts the variable name from a cut expression.

Usage

get_cut_name(thiscut, confounders.name)

Arguments

Value

Character vector of variable names.

Bootstrap Confidence Interval and Bias Correction Results

Description

Calculates confidence intervals and bias-corrected estimates for bootstrap results.

Usage

get_dfRes(
  Hobs,
  seHobs,
  H1_adj,
  H2_adj = NULL,
  ystar,
  cov_method = "standard",
  cov_trim = 0,
  est.scale = "hr",
  est.loghr = TRUE
)

Arguments

Value

Data.table with confidence intervals and estimates.

Generate Prediction Dataset with Subgroup Treatment Recommendation

Description

Creates a prediction dataset with a treatment recommendation flag based on the subgroup definition. Supports both label expressions (e.g., "\{er <= 0\}") and bare column names (e.g., "q3.1").

Usage

get_dfpred(df.predict, sg.harm, version = 1)

Arguments

Details

Each element of sg.harm is processed as follows:

1. Outer braces and leading ! are stripped.

2. If the result matches "var op value" (where op is one of <=, <, >=, >, ==, !=), the comparison is executed directly on df.predict[[var]].

3. Otherwise the expression is treated as a column name and membership is df.predict[[name]] == 1.

Value

Data frame with treatment recommendation flag (treat.recommend): 0 for harm subgroup, 1 for complement.

See Also

evaluate_comparison for the operator-dispatch logic, forestsearch for the main analysis function.

Extract HR from DGM (Backward Compatible)

Description

Extracts hazard ratios from DGM object, supporting both old and new formats. Also supports CDE (controlled direct effect) extraction for Table 5 of Leon et al. (2024) alignment (theta-ddagger).

Usage

get_dgm_hr(dgm, which = "hr_H")

Arguments

Value

Numeric hazard ratio value

Create DGM with Output File Path

Description

Wrapper function that creates a GBSG DGM and generates a standardized output file path for saving results.

Usage

get_dgm_with_output(
  model_harm,
  n,
  k_treat = 1,
  target_hr_harm = NULL,
  cens_type = "weibull",
  out_dir = NULL,
  file_prefix = "sim",
  file_suffix = "",
  include_hr_in_name = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Value

List with components:

dgm

The gbsg_dgm object

out_file

Character path to output file (NULL if out_dir is NULL)

k_inter

The k_inter value used (either calibrated or default)

Get Parameter with Default Fallback

Description

Safely retrieves a named element from a list, returning a default value if the element is missing or NULL.

Usage

get_param(args_list, param_name, default_value)

Arguments