Package 'forestsearch'

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

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

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

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)

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)

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

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

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
)

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

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
)

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

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
)

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)

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

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

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

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

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

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

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.

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

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

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)

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:

Description

Creates a publication-quality forest plot using ggplot2 for the CI panel and patchwork to assemble label and annotation columns alongside it. Unlike forestploter, fig.height maps directly to row density — row_height = fig.height / n_rows with no hidden scaling.

Usage

gg_forest(
  subgroups,
  est,
  lo,
  hi,
  cat_vec = NULL,
  cat_colours = NULL,
  annot = NULL,
  ref_line = 1,
  vert_lines = NULL,
  ref_col = "firebrick",
  ref_lty = "dashed",
  vert_col = "grey50",
  vert_lty = "dotted",
  xlim = NULL,
  ticks_at = NULL,
  tick_labels = NULL,
  xlog = TRUE,
  xlab = "Hazard Ratio",
  title = NULL,
  subtitle = NULL,
  footnote = NULL,
  point_size = 2.5,
  line_size = 0.8,
  point_shape = 21,
  base_size = 11,
  widths = NULL,
  row_expand = 0.6
)

Arguments

Value

A patchwork object. Render with print() or plot(). Control dimensions entirely via knitr chunk options fig.width / fig.height: row height = fig.height / n_rows.

Description

Identifies subgroups with differential treatment effect using generalized random forests (GRF) and policy trees. This function uses causal survival forests to identify heterogeneous treatment effects and policy trees to create interpretable subgroup definitions.

Usage

grf.subg.harm.survival(
  data,
  confounders.name,
  outcome.name,
  event.name,
  id.name,
  treat.name,
  frac.tau = 1,
  n.min = 60,
  dmin.grf = 0,
  RCT = TRUE,
  details = FALSE,
  sg.criterion = "mDiff",
  maxdepth = 2,
  seedit = 8316951,
  return_selected_cuts_only = FALSE
)

Arguments

Details

The return_selected_cuts_only parameter controls which cuts are returned:

FALSE (default)

Returns all cuts from all fitted trees (depths 1 to maxdepth). This provides the full set of candidate splits for downstream exploration and is the original behavior for backward compatibility.

TRUE

Returns only cuts from the tree at the depth that identified the "winning" subgroup meeting the dmin.grf criterion. This is useful when you want a focused set of cuts associated with the selected subgroup, reducing noise from non-selected trees.

When return_selected_cuts_only = TRUE and no subgroup meets the criteria, tree.cuts will be empty (character(0)).

Value

A list with GRF results, including:

Additional tree-specific cuts and objects (tree1, tree2, tree3) based on maxdepth

Description

Produces a templated text summary of the estimation properties table, automatically populating numerical results from the simulation output. Useful for reproducible vignettes where interpretation paragraphs should update when simulations are re-run.

Usage

interpret_estimation_table(
  results,
  dgm,
  analysis_method = "FSlg",
  n_sims = NULL,
  n_boots = 300,
  digits = 2,
  scenario = NULL,
  cat = TRUE
)

Arguments

Value

Invisibly returns the interpretation as a character string.

See Also

build_estimation_table, format_oc_results, get_dgm_hr

Description

Simulates multi-regional clinical trials and evaluates ForestSearch subgroup identification. Splits data by region into training and testing populations, identifies subgroups using ForestSearch on training data, and evaluates performance on the testing region.

Usage

mrct_region_sims(
  dgm,
  n_sims,
  n_sample = NULL,
  region_var = "z_regA",
  sg_focus = "minSG",
  maxk = 1,
  hr.threshold = 0.9,
  hr.consistency = 0.8,
  pconsistency.threshold = 0.9,
  confounders.name = NULL,
  conf_force = NULL,
  fs_args = list(),
  sim_args = list(rand_ratio = 1, draw_treatment = TRUE),
  analysis_time = 60,
  cens_adjust = 0,
  parallel_args = list(plan = "multisession", workers = NULL, show_message = TRUE),
  details = FALSE,
  verbose_n_sims = 2L,
  seed = NULL
)

Arguments

Details

Simulation Process

For each simulation:

1. Sample from super-population using simulate_from_dgm

2. Split by region_var into training and testing populations

3. Estimate HRs in ITT, training, and testing populations

4. Run forestsearch on training population

5. Apply identified subgroup to testing population

6. Calculate subgroup-specific estimates

Region Variable

The region_var parameter is used ONLY for splitting data into training/testing populations. It does not imply any prognostic effect. To include prognostic confounder effects, specify them when creating the DGM using create_dgm_for_mrct or generate_aft_dgm_flex.

Value

A data.table with simulation results containing:

sim

Simulation index

n_itt

ITT sample size

hr_itt

ITT hazard ratio (stratified if strat variable present)

hr_ittX

ITT hazard ratio stratified by region

n_train

Training (non-region A) sample size

hr_train

Training population hazard ratio

n_test

Testing (region A) sample size

hr_test

Testing population hazard ratio

any_found

Indicator: 1 if subgroup identified, 0 otherwise

sg_found

Character description of identified subgroup

n_sg

Subgroup sample size

hr_sg

Subgroup hazard ratio in testing population

POhr_sg

Potential outcome hazard ratio in subgroup (testing)

prev_sg

Subgroup prevalence (proportion of testing population)

n_sg_train

Subgroup sample size in training population

hr_sg_train

Subgroup hazard ratio in training population

POhr_sg_train

Potential outcome hazard ratio in subgroup (training)

hr_sg_null

Subgroup HR when found, NA otherwise

See Also

forestsearch for subgroup identification algorithm generate_aft_dgm_flex for DGM creation simulate_from_dgm for data simulation create_dgm_for_mrct for MRCT-specific DGM wrapper summaryout_mrct for summarizing simulation results

Description

Creates a visualization of the detection probability curve.

Usage

plot_detection_curve(
  curve_data,
  add_reference_lines = TRUE,
  add_threshold_line = TRUE,
  title = NULL,
  ...
)

Arguments

Value

Invisibly returns the input data.

Description

Creates Kaplan-Meier survival difference band plots comparing the identified ForestSearch subgroup (sg.harm) and its complement against the ITT population. This function wraps plotKM.band_subgroups() from the weightedsurv package, automatically extracting subgroup definitions from ForestSearch results.

Usage

plot_km_band_forestsearch(
  df,
  fs.est = NULL,
  sg_cols = NULL,
  sg_labels = NULL,
  sg_colors = NULL,
  itt_color = "azure3",
  outcome.name = "tte",
  event.name = "event",
  treat.name = "treat",
  xlabel = "Time",
  ylabel = "Survival differences",
  yseq_length = 5,
  draws_band = 1000,
  tau_add = NULL,
  by_risk = 6,
  risk_cex = 0.75,
  risk_delta = 0.035,
  risk_pad = 0.015,
  ymax_pad = 0.11,
  show_legend = TRUE,
  legend_pos = "topleft",
  legend_cex = 0.75,
  ref_subgroups = NULL,
  verbose = FALSE
)

Arguments

Details

This function simplifies the workflow of creating KM survival difference band plots for ForestSearch-identified subgroups. It can work in two modes:

Mode 1: With ForestSearch result (fs.est provided)

1. Extracts the subgroup definition from the ForestSearch result

2. Creates binary indicator columns (Qrecommend, Brecommend) in df

3. Generates appropriate labels from the subgroup definition

4. Calls plotKM.band_subgroups() with configured parameters

Mode 2: With pre-defined columns (sg_cols provided)

1. Uses existing indicator columns in df

2. Requires sg_labels and sg_colors to match sg_cols

The sg.harm subgroup (Qrecommend) represents patients with questionable treatment benefit (where treat.recommend == 0 in ForestSearch output). The complement (Brecommend) represents patients recommended for treatment.

Value

Invisibly returns a list containing:

df

The modified data frame with subgroup indicators

sg_cols

Character vector of subgroup column names used

sg_labels

Character vector of subgroup labels used

sg_colors

Character vector of colors used

sg_harm_definition

The subgroup definition extracted from fs.est

ref_subgroups

The reference subgroups list (if provided)

Subgroup Extraction

When fs.est is provided, the subgroup definition is extracted from:

• fs.est$grp.consistency$out_sg$sg.harm_label - Human-readable labels

• fs.est$sg.harm - Technical factor names (fallback)

• fs.est$df.est$treat.recommend - Subgroup membership indicator

Note

This function requires the weightedsurv package, which can be installed from GitHub: devtools::install_github("larry-leon/weightedsurv")

See Also

forestsearch for running the subgroup analysis plot_sg_weighted_km for weighted KM plots plot_sg_results for comprehensive subgroup visualization

Description

Creates weighted Kaplan-Meier survival curves for the identified subgroups (H and Hc) using the weightedsurv package, matching the pattern used in sg_consistency_out().

Usage

plot_sg_weighted_km(
  fs.est,
  fs_bc = NULL,
  outcome.name = "Y",
  event.name = "Event",
  treat.name = "Treat",
  by.risk = NULL,
  sg0_name = NULL,
  sg1_name = NULL,
  conf.int = TRUE,
  show.logrank = TRUE,
  show.cox = TRUE,
  show.cox.bc = TRUE,
  put.legend.lr = "topleft",
  ymax = 1.05,
  xmed.fraction = 0.65,
  hr_bc_position = "bottomright",
  hr_bc_cex = 0.725,
  title = NULL,
  verbose = FALSE
)

Arguments

Details

This function uses the exact same calling pattern as plot_subgroup() in the ForestSearch package. Column names are mapped internally to the standard names (Y, Event, Treat) expected by weightedsurv.

Subgroup definitions are automatically extracted from the forestsearch object if available:

• fs$grp.consistency$out_sg$sg.harm_label - Human-readable labels

• fs$sg.harm - Technical factor names (fallback)

HR display options controlled by show.cox and show.cox.bc:

• Both TRUE (default): Shows unadjusted HR from weightedsurv AND bias-corrected HR annotation

• show.cox = TRUE, show.cox.bc = FALSE: Shows only unadjusted HR

• show.cox = FALSE, show.cox.bc = TRUE: Shows only bias-corrected HR

• Both FALSE: Shows neither HR estimate

Value

Invisibly returns a list with subgroup data frames and counting data

Description

Plot Spline Treatment Effect Function

Usage

plot_spline_treatment_effect(dgm_result, add_points = TRUE)

Arguments

Value

No return value, called for side effects (produces a plot).

Description

Generates a comprehensive forest plot showing:

• ITT (Intent-to-Treat) population estimate

• Reference subgroups (e.g., by biomarker levels)

• Post-hoc identified subgroups with bias-corrected estimates

• Cross-validation agreement metrics as annotations

Usage

plot_subgroup_results_forestplot(
  fs_results,
  df_analysis,
  subgroup_list = NULL,
  outcome.name,
  event.name,
  treat.name,
  E.name = "Experimental",
  C.name = "Control",
  est.scale = "hr",
  xlog = TRUE,
  title_text = NULL,
  arrow_text = c("Favors Experimental", "Favors Control"),
  footnote_text = c("Eg 80% of training found SG: 70% of B (+) also B in CV testing"),
  xlim = c(0.25, 1.5),
  ticks_at = c(0.25, 0.7, 1, 1.5),
  show_cv_metrics = TRUE,
  cv_source = c("auto", "kfold", "oob", "both"),
  posthoc_colors = c("powderblue", "beige"),
  reference_colors = c("yellow", "powderblue"),
  ci_column_spaces = 20,
  conf.level = 0.95,
  theme = NULL
)

Arguments

Details

Creates a publication-ready forest plot displaying identified subgroups with hazard ratios, bias-corrected estimates, and cross-validation metrics. This wrapper integrates ForestSearch results with the forestploter package.

ForestSearch Labeling Convention

ForestSearch identifies subgroups based on hazard ratio thresholds:

• sg.harm: Contains the definition of the "harm" or "questionable" subgroup (H)

• treat.recommend == 0: Patient is IN the harm subgroup (H)

• treat.recommend == 1: Patient is in the COMPLEMENT subgroup (Hc, typically benefit)

For est.scale = "hr" (searching for harm):

• H (treat.recommend=0): Subgroup defined by sg.harm with elevated HR (harm/questionable)

• Hc (treat.recommend=1): Complement of sg.harm (potential benefit)

For est.scale = "1/hr" (searching for benefit):

• Roles are reversed: H becomes the benefit group

Value

A list containing:

plot

The forestploter grob object (can be rendered with plot())

data

The data frame used for the forest plot

row_types

Character vector of row types for styling reference

cv_metrics

Cross-validation metrics text (if available)

See Also

forestsearch for running the subgroup analysis forestsearch_bootstrap_dofuture for bootstrap bias correction forestsearch_Kfold for cross-validation create_forest_theme for customizing plot appearance render_forestplot for rendering the plot

Description

Dispatches to plot_sg_results for Kaplan-Meier curves, hazard-ratio forest plots, or combined panels.

Usage

## S3 method for class 'forestsearch'
plot(
  x,
  type = c("combined", "km", "forest", "summary"),
  outcome.name = "Y",
  event.name = "Event",
  treat.name = "Treat",
  ...
)

Arguments

Value

Invisibly returns the plot result from plot_sg_results.

See Also

plot_sg_results for full control over appearance, plot_sg_weighted_km for weighted KM curves, plot_subgroup_results_forestplot for publication-ready forest plots.

Description

Print method for cox_ahr_cde objects

Usage

## S3 method for class 'cox_ahr_cde'
print(x, ...)

Arguments

Value

Invisibly returns the input object.

Description

Displays a concise summary of ForestSearch results including the identified subgroup definition, consistency metrics, algorithm details, and computation time.

Usage

## S3 method for class 'forestsearch'
print(x, ...)

Arguments

Value

Invisibly returns x.

See Also

summary.forestsearch for detailed output, plot.forestsearch for visualization.

Description

Print Method for K-Fold CV Results

Usage

## S3 method for class 'fs_kfold'
print(x, ...)

Arguments

Value

Invisibly returns x.

Description

Print Method for Repeated K-Fold CV Results

Usage

## S3 method for class 'fs_tenfold'
print(x, ...)

Arguments

Value

Invisibly returns x.

Description

Print Method for gbsg_dgm Objects

Usage

## S3 method for class 'gbsg_dgm'
print(x, ...)

Arguments

Value

Invisibly returns x.

Examples

dgm <- setup_gbsg_dgm(model = "alt", verbose = FALSE)
print(dgm)

Description

Renders a forest plot from plot_subgroup_results_forestplot().

Usage

render_forestplot(x, newpage = TRUE)

Arguments

Details

To control plot sizing, create a custom theme using create_forest_theme() and pass it to plot_subgroup_results_forestplot():

my_theme <- create_forest_theme(base_size = 14, row_padding = c(6, 4))

result <- plot_subgroup_results_forestplot(..., theme = my_theme)

render_forestplot(result)

Value

Invisibly returns the grob object.

Description

Converts a data frame of pre-computed reference simulation results (e.g., digitized from a published LaTeX table) into a styled gt table. This is useful for displaying published benchmark results alongside new simulation output within vignettes or reports.

Usage

render_reference_table(
  ref_df,
  title = "Reference Simulation Results",
  subtitle = NULL,
  bold_threshold = 0.05
)

Arguments

Value

A gt table object.

Examples

ref <- data.frame(
  Scenario = "M1 Null: N=700",
  Metric   = "any(H)",
  FS       = 0.02,
  FSlg     = 0.03,
  GRF      = 0.25
)
render_reference_table(ref, title = "Reference Results")

Description

General replacement for the legacy run_simulation_analysis() that was coupled to simulate_from_gbsg_dgm() and GBSG-specific column names. This version calls simulate_from_dgm and accepts explicit column-name parameters, making it applicable to any DGM built with generate_aft_dgm_flex.

Usage

run_simulation_analysis(
  sim_id,
  dgm,
  n_sample,
  analysis_time = Inf,
  cens_adjust = 0,
  max_follow = NULL,
  muC_adj = NULL,
  confounders_base = c("v1", "v2", "v3", "v4", "v5", "v6", "v7"),
  n_add_noise = 0L,
  outcome_name = "y_sim",
  event_name = "event_sim",
  treat_name = "treat_sim",
  harm_col = "flag_harm",
  run_fs = TRUE,
  run_fs_grf = TRUE,
  run_grf = TRUE,
  fs_params = list(),
  grf_params = list(),
  cox_formula = NULL,
  cox_formula_adj = NULL,
  n_sims_total = NULL,
  seed_base = 8316951L,
  verbose = FALSE,
  verbose_n = NULL,
  debug = FALSE
)

Arguments

Value

A data.table with one row per analysis method, containing subgroup size, HR, AHR, CDE, and classification metrics.

See Also

simulate_from_dgm, generate_aft_dgm_flex, setup_gbsg_dgm

Description

Identifies the optimal subgroup according to the specified criterion

Usage

select_best_subgroup(values, sg.criterion, dmin.grf, n.max)

Arguments

Value

Data frame row with best subgroup or NULL if none found

Examples

vals <- data.frame(diff = c(8.5, 6.2, 3.1), Nsg = c(120, 95, 80))
select_best_subgroup(values = vals, sg.criterion = "mDiff",
                     dmin.grf = 6, n.max = 500)

Description

Creates a GBSG-based data generating mechanism that is fully compatible with simulate_from_dgm. This is the replacement for create_gbsg_dgm(): it accepts exactly the same arguments and produces the same numeric output, but returns an object of class "aft_dgm_flex" instead of "gbsg_dgm".

Usage

setup_gbsg_dgm(
  model = c("alt", "null"),
  k_treat = 1,
  k_inter = 1,
  k_z3 = 1,
  z1_quantile = 0.25,
  n_super = 5000L,
  cens_type = c("weibull", "uniform"),
  use_rand_params = FALSE,
  seed = 8316951L,
  verbose = FALSE
)

Arguments

Details

Internally the function calls create_gbsg_dgm() and then:

1. Adds a df_super field with column names aligned to simulate_from_dgm() conventions (lin_pred_1, lin_pred_0, lin_pred_cens_1, lin_pred_cens_0, flag_harm).

2. Adds a model_params$tau field (= model_params$sigma) and a model_params$censoring sub-list.

3. Sets class to c("aft_dgm_flex", "gbsg_dgm", "list").

The original df_super_rand field is kept so that compute_dgm_cde() and print.gbsg_dgm continue to work.

Value

An object of class c("aft_dgm_flex", "gbsg_dgm", "list") with all fields from create_gbsg_dgm() plus:

df_super

Super-population data frame with simulate_from_dgm()-compatible column names.

model_params$tau

Copy of model_params$sigma.

model_params$censoring

Sub-list with type, mu, tau for the censoring model.

See Also

create_gbsg_dgm, simulate_from_dgm, compute_dgm_cde

Examples

dgm <- setup_gbsg_dgm(model = "alt", k_inter = 2, verbose = FALSE)
dgm <- compute_dgm_cde(dgm)
print(dgm)
sim <- simulate_from_dgm(dgm, n = 400, seed = 1)

Description

Returns formatted summary tables for subgroups using the gt package, with search metadata and customizable decimal precision. Produces two tables: a treatment effect estimates table and an identified subgroups table, each with fully customizable titles and subtitles.

Usage

sg_tables(
  fs,
  which_df = "est",
  est_title = "Treatment Effect Estimates",
  est_caption = "Training data estimates",
  sg_title = "Identified Subgroups",
  sg_subtitle = NULL,
  potentialOutcome.name = NULL,
  hr_1a = NA,
  hr_0a = NA,
  ndecimals = 3,
  include_search_info = TRUE,
  font_size = 12
)

Arguments

Value

List with gt tables for estimates, subgroups, and optionally search info.

Description

Generates simulated survival data from a previously created AFT data generating mechanism (DGM). Samples from the super population and generates survival times with specified censoring.

Usage

simulate_from_dgm(
  dgm,
  n = NULL,
  rand_ratio = 1,
  entry_var = NULL,
  max_entry = 24,
  analysis_time = 48,
  cens_adjust = 0,
  draw_treatment = TRUE,
  seed = NULL,
  strata_rand = NULL,
  hrz_crit = NULL,
  keep_rand = FALSE,
  time_eos = NULL
)

Arguments

Details

Time-scale consistency

All time parameters (analysis_time, max_entry, time_eos) must be expressed in the same units as outcome_var supplied to generate_aft_dgm_flex(). A common error is building the DGM on days (e.g. rfstime) and then passing analysis_time in months, which causes follow-up windows far shorter than the DGM event-time scale and produces universal administrative censoring (event_sim = 0 for all subjects).

Verify with: exp(dgm$model_params$mu) — the implied median event time should be plausible given your analysis_time.

n = NULL path

When n = NULL the entire super population is used as-is, with no staggered entry and no administrative censoring (follow_up = Inf). Treatment assignments and linear predictors already stored in dgm$df_super are retained unchanged.

Censoring adjustment

cens_adjust shifts the log-scale location parameter of the censoring distribution:

• cens_adjust = log(2) doubles expected censoring times.

• cens_adjust = log(0.5) halves expected censoring times.

Value

A data.frame with columns:

id

Subject identifier.

treat

Original treatment from super population.

treat_sim

Simulated treatment assignment.

flag_harm

Subgroup indicator (1 = all subgroup conditions met).

z_*

Covariate values.

lin_pred_1, lin_pred_0

Counterfactual log-time linear predictors.

y_sim

Observed survival time (min(T, C)).

event_sim

Event indicator (1 = event, 0 = censored).

t_true

Latent true survival time (pre-censoring).

c_time

Effective censoring time (post admin-censoring).

hrz_flag

(Optional) Individual harm-zone indicator.

rand_order

(Optional) Randomisation sequence index.

See Also

generate_aft_dgm_flex, check_censoring_dgm

Examples

dgm <- setup_gbsg_dgm(model = "null", verbose = FALSE)
sim_data <- simulate_from_dgm(dgm, n = 200, seed = 42)
dim(sim_data)
head(sim_data[, c("y_sim", "event_sim", "treat_sim")])