Package {forestsearch}

forestsearch: Exploratory Subgroup Identification in Survival and GLM Outcomes

Description

Implements statistical methods for exploratory subgroup identification in clinical trials. 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. Supports survival endpoints (Cox proportional hazards), binary outcomes (log odds ratio, log relative risk, risk difference), continuous outcomes (mean difference), and count / rate outcomes (log incidence rate ratio via Poisson, quasi-Poisson, or negative-binomial GLMs with optional person-time offset). 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. Methods are described in Leon et al. (2024) doi:10.1002/sim.10163. Post-selection inference for the identified subgroup follows León and Anderson (2026) doi:10.48550/arXiv.2609.38361.

Author(s)

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

Authors:

• 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

Normalise adjustment terms (mirrors forestsearch:::.fs_adjust_terms)

Description

Normalise adjustment terms (mirrors forestsearch:::.fs_adjust_terms)

Usage

.consistency_adj_terms(adjust_covariates)

Per-subject Cox pieces for the consistency approximation

Description

Fits the treatment-only (or covariate-adjusted) Cox model once on a subgroup and returns the treatment log hazard ratio, the per-subject treatment dfbeta influence contributions, and the split-perturbation SD.

Usage

.consistency_cox_pieces(
  df,
  tte.name = "Y",
  event.name = "Event",
  treat.name = "Treat",
  cox_init = 0,
  adjust_covariates = NULL
)

Arguments

Value

List with beta_hat, dfbeta, sigma_D, n, d, log_scale (TRUE), measure ("HR"); or NULL if degenerate.

Per-subject GLM pieces for the consistency approximation

Description

Fits an unweighted, unadjusted treatment-effect model and returns the treatment coefficient, its per-subject dfbeta influence contributions, and the split-perturbation SD.

Usage

.consistency_glm_pieces(
  df,
  outcome_type,
  effect_measure = NULL,
  treat.name = "Treat",
  outcome.name = "Y",
  offset.name = NULL,
  adjust_covariates = NULL,
  adverse_outcome = TRUE
)

Arguments

Details

Supported: "OR" (logistic), "RR" (log-binomial with Poisson fallback), "RD" (identity-link binomial), "MD" (OLS), "IRR" (Poisson with offset(log(time))). "IRD" is not a single coefficient and returns NULL, and the caller falls back to splitting.

Value

List with beta_hat, dfbeta, sigma_D, n, d, log_scale, measure; or NULL if unsupported/degenerate/non-convergent.

Propensity adjustment is not represented here

This function has no ps_method or ps_adjust_method argument and reads no sw / ps_hat / ips_covar column, so under ps_method != "none" it fits the effect as if unadjusted rather than declining to fit it. It cannot do otherwise: nothing in its signature tells it which case it is in.

The consequence is that the identifier ranks candidates on the IPTW or g-computation adjusted effect built by make_effect_estimator(), while every quantity derived here – \hat\beta(g), db_{g,i}, \sigma_{D,g}, and therefore the whole multiplier-resampling correction – describes the unadjusted functional. forestsearch() also forbids adjust_covariates alongside ps_method != "none", so the model fitted here is fully unadjusted, not differently adjusted.

Earlier versions of this documentation claimed that propensity-adjusted effects "return NULL (the caller falls back to splitting)". That was never true and was not implementable as written. Whether to extend the influence path, add the argument purely so the claim can be honoured, or refuse the combination outright is unresolved.

Consistency proportion via literal splitting (single-stage)

Description

Runs n.splits Bernoulli splits through run_single_consistency_split() and returns the rounded proportion of consistent splits, or NA_real_ when too few valid splits are obtained. Survival and GLM (estimator_fn) paths.

Usage

.consistency_via_splits(
  df.x,
  N.x,
  n.splits,
  hr.consistency,
  cox_init,
  estimator_fn,
  consistency_threshold,
  adjust_covariates,
  pconsistency.digits,
  m = NA,
  details = FALSE
)

Internal Workhorse for Single-Theta Detection Probability

Description

Internal Workhorse for Single-Theta Detection Probability

Usage

.detect_prob_glm_single(beta, sigma2_s, k1, k2, method, n_mc, seed)

Arguments

Value

Numeric probability.

Fit Separate Regressions by Treatment Arm

Description

Fit Separate Regressions by Treatment Arm

Usage

.fit_arms(
  Y,
  W,
  X,
  poly_order,
  covariate_select,
  t_threshold,
  regression = "ols",
  offset = NULL
)

Apply multiplier resampling with forestsearch's defaulting rules.

Description

Thin wrapper around fs_mr_inference() shared by the consistency, DINA, and GRF branches so the call (and mr_inference_args overrides) stays identical across methods. Wrapped in tryCatch so a reconstruction failure yields NULL (the branch's normal output is unaffected) rather than aborting.

Usage

.fs_apply_mr(
  df,
  candidates,
  selected_members,
  spec,
  admission,
  effect_neighborhood,
  reselection_default,
  selection_rule_default = "neighborhood",
  mr_inference_args = list(),
  seedit = NULL
)

Arguments

Forwarded argument set

The wrapper forwards the full fs_mr_inference() argument set, so the DINA and GRF branches are controllable through mr_inference_args exactly as the consistency branch is through its own direct call. Before this, ten arguments were dropped here – return_reselection, field_R_out, field_R_in, field_uniform, field_M_cap, field_complement, field_decompose, field_scale_complement, ij_residual, field_recovery – which left field_decompose and field_recovery unreachable on those engines and the other three inert while their defaults happened to coincide with the intended values: a run's meta recorded them as set while they controlled nothing.

Defaults are unchanged on every branch. Each newly forwarded argument defaults to fs_mr_inference()'s own formal default, read from formals(fs_mr_inference) at call time rather than restated here, so the two cannot drift and a call that asks for nothing is byte-identical to the same call before the change. match.arg()-style defaults are forwarded as their whole candidate vector, which fs_mr_inference() resolves exactly as it would have resolved a missing argument. ci_method keeps this wrapper's own "ij" default, which is deliberately not the consistency branch's "field"; changing it is a separate decision.

Assemble the influence matrix from candidate membership lists

Description

For each candidate (a vector of row indices into df), fit the treatment effect once and place its dfbeta into the candidate's rows (zero elsewhere). Candidates whose fit fails or whose dfbeta length does not match are dropped.

Usage

.fs_mr_assemble(df, candidates, spec)

Arguments

Value

List with B (⁠N x S_kept⁠), beta_hat, sigma_D, sizes, names, keep, log_scale.

Default near-null harm-confirmation threshold for a measure

Description

The value t_confirm takes when left NULL: the null effect on the effect scale, i.e. no effect at all rather than the screening threshold.

Usage

.fs_mr_confirm_null(effect_measure)

Arguments

Value

1 for ratio measures (effect scale), 0 for differences.

Build candidate-membership lists from a (v1,d1,c1,v2,d2,c2) table

Description

Build candidate-membership lists from a (v1,d1,c1,v2,d2,c2) table

Usage

.fs_mr_family_from_table(df, tab, op_right = ">=", n_min = 1L, digits = 17L)

Arguments

Value

Named list of integer row-index vectors (one per candidate).

Infinitesimal-jackknife variance of the de-biased (bagged) estimate

Description

Implements Leon et al. (2024) Eq. (VInfJ)/(VInfJ_bc) in multiplier form. The multiplier matrix Xi supplies the centered bootstrap multiplicities K^*_{bi} - \bar K^*_i, and r is the per-draw residual r_b = \hat\beta(\widehat H) - \eta^*_b(\widehat H^*_b) - \eta^*_b(\widehat H) - \hat\beta^*(\widehat H), i.e. ⁠(selection_bias + fixed_bias) - D_{H*_b}(b) - D_{H}(b)⁠, evaluated only on the draws in ok that produced a re-selected winner.

Usage

.fs_mr_ij_var(Xi, r, ok)

Arguments

Value

List with tilde_V (raw IJ), hat_V (Wager 2014 bias-corrected), and B_ok (number of usable draws).

Row-index membership for one candidate conjunction.

Description

Row-index membership for one candidate conjunction.

Usage

.fs_mr_members_from_conj(df, cj)

Arguments

Value

Integer row indices of subjects in the (harm) subgroup.

Mean-zero unit-variance multiplier vectors (one column per draw)

Description

Mean-zero unit-variance multiplier vectors (one column per draw)

Usage

.fs_mr_multipliers(n, draws, type)

Per-candidate influence pieces (Cox / GLM dispatch)

Description

Thin dispatch onto the existing forestsearch internals so MR uses the identical treatment dfbeta, beta_hat, and robust sigma_D the resample consistency engine uses.

Usage

.fs_mr_pieces(df_sub, spec)

Arguments

Value

Output of .consistency_glm_pieces() / .consistency_cox_pieces(), or NULL.

Map sg_focus to MR's re-selection rule, faithfully per engine.

Description

MR must re-select under the same rule the search used. The only rule that differs by engine is the "hr"/"eff" focus: the consistency search ranks by consistency rate, whose MR analog is "maxcons"; DINA/GRF rank by effect, whose analog is "maxeff". Every other focus maps identically (maxSG/minSG pass through; hrMaxSG/hrMinSG -> the size-within-effect- neighborhood rules effMaxSG/effMinSG). This makes sg_focus the single source of truth, so callers never need to specify reselection by hand.

Usage

.fs_mr_reselection_from_focus(sg_focus, engine = c("effect", "consistency"))

Arguments

Resolve a de-biased SE from the IJ variance, with graceful fallback

Description

Prefers the bias-corrected IJ variance; if it is non-positive (too few draws), falls back to the raw IJ variance, then to the subgroup robust SE.

Usage

.fs_mr_se_from_ij(ij, se_fallback)

Arguments

Value

List with se, var, and source ("ij", "ij_raw", or "wald_fallback").

Apply a selection rule among passing candidates on one draw

Description

For the size-within-band rules (effMaxSG/effMinSG) the inclusion band is built to match the search's selection_rule:

• "neighborhood" - effect within nbhd of the max, on the NATURAL effect scale and multiplicative (eff >= (1 - nbhd) * max(eff)), matching sort_subgroups() exactly (G3);

• "pareto" - the 2-D non-dominated set in (effect, size), via the same dominance core the search uses (.pareto_dominated_xy()) (G2);

• "both" - the intersection of the two. beta is on the working scale (log for ratio measures); log_scale controls the conversion to the natural effect for the band.

Usage

.fs_mr_select(
  beta,
  zcons,
  sizes,
  passers,
  rule,
  nbhd,
  selection_rule = "neighborhood",
  log_scale = TRUE
)

Build the Cox/GLM spec for MR from forestsearch()'s resolved columns.

Description

DINA/GRF operate on the user-facing analysis frame df (not the standardized df.fs), so survival uses the user column names here – unlike the consistency hook, which uses the internal Y/Event/Treat names on df.fs.

Usage

.fs_mr_spec(
  outcome_type,
  effect_measure,
  treat.name,
  outcome.name,
  event.name,
  offset.name,
  adjust_covariates,
  adverse_outcome,
  df
)

Can the GLM resampling approximation represent this effect measure?

Description

The resampling approximation requires the treatment effect to be a single model coefficient with a well-defined influence function. That holds for OR/RR (logistic / log-binomial), RD (identity-link binomial), MD (OLS), and IRR (Poisson with offset). It does not hold for IRD (a delta-method rate difference) or propensity-adjusted (IPTW / G-computation) effects, for which the caller falls back to literal splitting.

Usage

.glm_resample_supported(spec)

Arguments

Value

Logical scalar.

Re-format Numeric Cells in FSsg_tab at a Lower Precision

Description

Parses formatted-string cells in FSsg_tab back to numerics and re-formats at the requested digits. Designed for display-time precision control in summarize_bootstrap_results.

Usage

.reformat_FSsg_tab_digits(tab, digits)

Arguments

Details

Cells are identified by column name pattern, not content sniffing. Three cell formats are recognized:

1. Single value with optional sign and percent suffix (Mean(C)/(T), Rate(C)/(T), Diff): "-56.4321", "+12.5", "42.1%".

2. Three-value CI (effect-measure column, "\{measure\} (95% CI)"): "82.7831 (14.8067, 150.7585)".

3. Three-value CI without space (bias-corrected "\{measure\}*" column): "72.4500 (-61.5400,206.4400)".

Columns NOT touched: Subgroup (text label), n and n1 (count-with-percent format like "87 (8.0%)" — different semantics).

Robustness contract: cells that don't match any expected pattern are returned unchanged. Construction-time precision must be at least digits for the reformat to be informative; if construction was at lower precision, this function will zero-pad to digits but cannot recover information.

Value

Reformatted tab, same shape.

Select Covariates for the Crump Test

Description

Returns a list with X (selected columns) and selected_names.

Usage

.select_covariates(Y, W, X, method, threshold)

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.

Examples

## Not run: 
# CV_sgs() is called internally by forestsearch_KfoldOut().
# See forestsearch_Kfold() and forestsearch_KfoldOut() for the entry points.

## End(Not run)

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.

Examples

## Not run: 
labels <- c("age <= 50", "tumor size <= 20", "nodes <= 3")
FS_labels("q1.1", labels)
FS_labels("q1.0", labels)
FS_labels("q2", labels)
FS_labels(c("q1", "q3"), labels)

## End(Not run)

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.

Examples

## Not run: 
library(survival)
df <- data.frame(
  tte          = gbsg$rfstime / 30.4375,
  event        = gbsg$status,
  treat        = gbsg$hormon,
  treat.recommend = as.integer(gbsg$er > 0)
)
SG_tab_estimates(df, SG_flag = "ITT",
                 outcome.name = "tte",
                 event.name   = "event",
                 treat.name   = "treat")

## End(Not run)

Subgroup Summary Table for GLM Outcomes

Description

GLM counterpart to SG_tab_estimates. Produces a summary data frame with per-arm outcome rates and the GLM effect estimate for each subgroup.

Usage

SG_tab_estimates_glm(
  df,
  SG_flag,
  outcome.name,
  treat.name,
  estimator_fn = NULL,
  effect_measure = "RD",
  outcome_type = "binary",
  effect_a_1 = NA,
  effect_a_0 = NA,
  sg1_name = "Recommend",
  sg0_name = "Questionable",
  est.scale = "hr",
  digits = 4
)

Arguments

Value

Data frame. One row when SG_flag = "ITT"; two rows (sg0 then sg1) when SG_flag = "treat.recommend". Columns: Subgroup, n, n1, Rate(C) or Mean(C) (binary vs continuous/count), Rate(T) or Mean(T), Diff, and the effect-measure header "\{effect_measure\} (95% CI)". When effect_a_* is supplied, an additional "\{effect_measure\}*" column carries the bias-corrected estimate.

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

Examples

## Not run: 
# Default (no trimming) — backward compatible
plot_results <- SGplot_estimates(
  results_alt,
  label_training = "Non-Region A, ITT",
  label_itt = "Overall, ITT",
  label_testing = "Region A, ITT",
  label_sg = "Region A, identified subgroup"
)
print(plot_results$plot_estimates)

# With 1% symmetric trimming
plot_results <- SGplot_estimates(
  results_alt,
  label_training = "Non-AP, ITT",
  label_itt = "Overall, ITT",
  label_testing = "AP, ITT",
  label_sg = "AP, identified subgroup",
  trim_fraction = 0.01
)
print(plot_results$plot_estimates)
# Inspect trimming diagnostics
print(plot_results$trim_info)

# With explicit y-axis limits
plot_results <- SGplot_estimates(
  results_alt,
  trim_fraction = 0.01,
  ylim = c(0.2, 3.0)
)

## End(Not run)

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
)

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

Description

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

Usage

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

Arguments

Value

Character vector of z_-prefixed covariate column names.

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.

Examples

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

## End(Not run)

Analyze subgroup for GLM outcomes

Description

GLM counterpart to analyze_subgroup. Computes per-arm outcome rates, sample sizes, and the treatment effect estimate (RD, OR, etc.) using the estimator closure.

Usage

analyze_subgroup_glm(
  df_sub,
  outcome.name,
  treat.name,
  subgroup_name,
  effect_a = NA,
  estimator_fn = NULL,
  effect_measure = "RD",
  outcome_type = "binary",
  N,
  digits = 4
)

Arguments

Value

One-row data frame of formatted subgroup results. Column names are placeholders (Subgroup, n, n1, rate0, rate1, Diff, effect, optionally effect_a) and are overwritten by the calling function (SG_tab_estimates_glm) with the user-facing labels (e.g., "Rate(C)" vs "Mean(C)", the effect-measure header). Returning a data frame – rather than a named character vector via c(...) – ensures type stability across the ITT and subgroup call paths and prevents downstream rendering failures in consumers that go through as.data.frame() (e.g., gt).

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

Random-benchmark subgroup specification

Description

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

Usage

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

Arguments

Value

An object of class "benchmark_spec".

See Also

run_subgroup_sims()

Examples

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

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,
  estimator_fn = NULL,
  effect_measure = NULL,
  boot_index_mat = NULL,
  mr_in_replicates = FALSE
)

Arguments

Value

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

boot_id

Integer. Bootstrap iteration number (1 to nb_boots)

H_biasadj_1

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

H_biasadj_2

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

Hc_biasadj_1

Bias-corrected estimate for H^c using method 1

Hc_biasadj_2

Bias-corrected estimate for H^c using method 2

max_sg_est

Numeric. Maximum subgroup hazard ratio found

L

Integer. Number of candidate factors evaluated

max_count

Integer. Maximum number of factor combinations

events_H_0

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

events_H_1

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

events_Hc_0

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

events_Hc_1

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

events_Hstar_0

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

events_Hstar_1

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

events_Hcstar_0

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

events_Hcstar_1

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

tmins_search

Numeric. Minutes spent on subgroup search in this iteration

tmins_iteration

Numeric. Total minutes for this bootstrap iteration

Pcons

Numeric. Consistency p-value for top subgroup

any_found

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

hr_sg

Numeric. Hazard ratio for top subgroup

N_sg

Integer. Sample size of top subgroup

E_sg

Integer. Number of events in top subgroup

K_sg

Integer. Number of factors defining top subgroup

g_sg

Numeric. Subgroup group ID

m_sg

Numeric. Subgroup index

M.1

Character. First factor label

M.2

Character. Second factor label

M.3

Character. Third factor label

M.4

Character. Fourth factor label

M.5

Character. Fifth factor label

M.6

Character. Sixth factor label

M.7

Character. Seventh factor label

grf_cuts_b

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

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

Bias Correction Methods

Two bias correction approaches are implemented:

1. Method 1 (Simple Optimism):

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

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

2. Method 2 (Double Bootstrap):

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

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

where:

• H_obs: Original subgroup HR on original data

• H_star: Original subgroup HR on bootstrap data

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

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

Computational Details

• Uses doFuture backend for parallel execution (configured externally)

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

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

• Sequential execution within each bootstrap prevents nested parallelization

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

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

Bootstrap Configuration

Each bootstrap iteration modifies ForestSearch arguments to:

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

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

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

Performance Considerations

• Typical runtime: 1-5 seconds per bootstrap iteration

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

• Memory usage scales with dataset size and number of workers

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

Error Handling

The function gracefully handles three failure modes:

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

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

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

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

Note

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

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

• Packages listed in BOOTSTRAP_REQUIRED_PACKAGES to be installed

• Proper parallel backend setup via setup_parallel_SGcons

See Also

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

Examples

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

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

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

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

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

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

# 6. Check results
summary(boot_results)

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

## End(Not run)

Bootstrap Ystar Count Matrix

Description

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

Usage

bootstrap_ystar(df, nb_boots, seed = 8316951L)

Arguments

Value

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

Implementation note (v0.2.0)

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

Reproducibility note

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

Examples

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

## End(Not run)

Build Polynomial Basis Matrix

Description

Build Polynomial Basis Matrix

Usage

build_basis(X, poly_order = 1L, include_intercept = TRUE)

Arguments

Value

Numeric matrix (N x K).

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,
  subgroup_notation = c("harm", "benefit")
)

Arguments

Details

For each scenario the function computes:

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

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

• sens(Hc): Mean specificity.

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

• ppv(Hc): Mean negative predictive value.

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

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

Value

A gt table object.

See Also

format_oc_results, summarize_simulation_results

Examples

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

build_classification_table(scenarios, n_sims = 100)

## End(Not run)

Build Cox Model Formula

Description

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

Usage

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

Arguments

Value

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

Examples

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

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,
  subgroup_notation = c("harm", "benefit"),
  trim_threshold = 1000,
  trim_fraction = 0.01
)

Arguments

Details

Uses the paper's notation conventions:

• theta-dagger: Marginal (causal) HR truth

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

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

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

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

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

• Avg: Mean of the estimates across estimable simulations.

• SD: Empirical standard deviation.

• Min / Max: Range.

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

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

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

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

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

Value

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

See Also

build_classification_table, format_oc_results, get_dgm_hr

Examples

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

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

## End(Not run)

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 L_eff from Simulation Results

Description

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

Usage

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

Arguments

Details

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

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

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

Value

A list of class "leff_calibration" with components:

C

Numeric. Fitted constant.

alpha

Numeric. Fitted power-law exponent.

n_min

Integer. Minimum subgroup size.

r_squared

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

data

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

fit

The lm object.

Examples

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

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,
  tol_rel = 2.5,
  n_eval = 2000,
  verbose = TRUE,
  ...
)

Arguments

Details

Two calibration targets are supported:

"rate"

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

"km_median"

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

How the objective function works

At each candidate cens_adjust value, the objective function:

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

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

3. Returns sim_metric - ref_metric.

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

Monotonicity

The objective is monotone in cens_adjust for both targets:

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

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

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

Stochastic noise

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

Value

A named list with elements:

cens_adjust

Calibrated cens_adjust value.

target

Calibration target used.

ref_value

Reference metric value from dgm$df_super.

sim_value

Achieved metric value in simulated data at the calibrated cens_adjust.

residual

Absolute difference between sim_value and ref_value.

iterations

Number of uniroot iterations.

diagnostic

Output of check_censoring_dgm at the calibrated value (invisibly).

See Also

simulate_from_dgm, check_censoring_dgm, generate_aft_dgm_flex

Examples

## Not run: 
library(survival)

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

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

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

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

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

## End(Not run)

Calibrate GLM Interaction for a Target Subgroup Effect Size

Description

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

Usage

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

Arguments

Details

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

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

Value

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

See Also

generate_glm_dgm, simulate_from_glm_dgm

Examples

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

## End(Not run)

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,
  tol_rel = 2.5,
  use_ahr = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

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

Value

Numeric value of k_inter that achieves the target HR

Examples

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

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

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

## End(Not run)

Calibrate k_treat for Target Overall Hazard Ratio

Description

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

Usage

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

Arguments

Details

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

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

Value

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

See Also

calibrate_k_inter, calibrate_cens_adjust, generate_aft_dgm_flex

Examples

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

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

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

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

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

## End(Not run)

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

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

## End(Not run)

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.

Format a numeric threshold for a cut expression

Description

Renders a (already-rounded) numeric threshold as a clean string with no trailing zeros or floating-point noise, for use in "var <= value" cut expressions. Integer-valued thresholds print without a decimal point.

Usage

collapse_cuts_fmt(v, digits = 0L)

Arguments

Value

Character scalar.

Round half away from zero

Description

Base R round() uses round-half-to-even ("banker's rounding"), so round(0.5) == 0 and round(74.5) == 74. Cut coarsening rounds representative thresholds half up (away from zero) so that, e.g., a cluster centroid of 74.5 maps to 75, matching the intuitive "nearest integer" rule.

Usage

collapse_cuts_round(x, digits = 0L)

Arguments

Value

Numeric vector rounded half away from zero.

Collapse near-redundant continuous candidate cuts

Description

Continuous covariates can generate many candidate threshold cuts that are practically redundant – e.g. "age <= 35.7" and "age <= 35", or "wtkg <= 75.2" and "wtkg <= 74.5" – because the candidate pool unions quantile cuts, GRF / DINA splits, and calibrated-DGM cuts at different precisions. This helper merges such near-duplicates to a single representative threshold, leaving categorical and indicator cuts untouched.

Usage

collapse_redundant_cuts(
  cuts,
  df,
  confounders.name = NULL,
  c_band = 1,
  safety_tol = 0.05,
  digits = 0L,
  cont.cutoff = 4,
  details = FALSE,
  protect = character(0)
)

Arguments

Details

For each (variable, operator) group of continuous-variable cuts, distinct thresholds are single-linkage clustered using a per-variable band band = c_band * sd(x) / sqrt(n) (the standard error of the variable: thresholds finer than ~1 SE are not statistically resolvable). Each cluster collapses to one representative, the cluster mean rounded half-up to digits places (see collapse_cuts_round()). Rounding alone can also merge singletons that round to the same value (e.g. 75.2 and 74.5 both round to 75 at digits = 0); de-duplication then drops the duplicate.

Safety check: a cluster is collapsed only if replacing each member threshold by the representative changes subgroup membership for no more than safety_tol subjects (a fraction of n when safety_tol < 1, an absolute count when safety_tol >= 1). Clusters that would move more than that are kept unchanged, so coarsening never silently redefines a candidate subgroup by more than the tolerance.

Cut expressions that are categorical, indicator-valued (fewer than cont.cutoff unique values), bare variable names, or equality tests (==) pass through unchanged. Cuts listed in protect are also passed through verbatim: they are model-identified screening candidates (e.g. GRF/DINA cuts) rather than points on a default continuous grid, so coarsening them would defeat the purpose of the screening engine.

Value

Character vector of de-duplicated cut expressions (length <= input).

Collect Results from foreach with Error Handling

Description

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

Usage

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

Arguments

Details

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

Design rationale for the reporting strategy:

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

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

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

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

Value

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

n_failed

Integer. Number of tasks that returned errors.

n_total

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

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

Preserving per-replicate diagnostics

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

See Also

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

Examples

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

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

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.

Examples

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

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

## End(Not run)

Compare Subgroup Membership Across Two Analyses

Description

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

Usage

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

Arguments

Value

A list with:

table

A gt cross-tabulation table.

crosstab

The raw cross-tabulation matrix.

concordance

Agreement, kappa, counts.

membership

Per-subject assignments from both analyses.

Examples

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

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

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

## End(Not run)

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

Examples

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

## End(Not run)

Compare forestsearch Runs Across Selection-Rule Combinations

Description

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

Usage

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

Arguments

Details

For each combination the wrapper:

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

2. Optionally computes the frontier CIs via compute_frontier_cis.

3. Builds the Pareto-frontier plot via plot_pareto_frontier.

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

Value

A list with class "forestsearch_comparison" and components:

combos

data.frame with columns sg_focus, selection_rule, label.

fs

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

ci_tab

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

plots

Named list of ggplot objects, one per combination.

plot_grid

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

plot_combined

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

plot_combined_subsets

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

combined_skip_reason

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

console

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

diagnostics

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

errors

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

Tuple semantics

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

Console capture

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

See Also

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

Examples

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

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

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

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

## End(Not run)

Compare two extreme-subgroups simulation studies

Description

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

Usage

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

Arguments

Details

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

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

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

Value

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

See Also

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

Examples

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

## End(Not run)

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 Marginal Causal Effect from GLM Potential Outcomes

Description

Computes the marginal (average) treatment effect from individual-level potential outcomes stored in GLM simulation data. This is the GLM analogue of compute_ahr for survival data.

Usage

compute_aor(df, subset_indicator = NULL, effect_measure = "OR")

Arguments

Details

For binary outcomes (p0, p1):

• OR: odds(mean(p1)) / odds(mean(p0))

• RD: mean(p1) - mean(p0)

• RR: mean(p1) / mean(p0)

For continuous or count outcomes (mu0, mu1):

• MD: mean(mu1) - mean(mu0)

• IRR: mean(mu1) / mean(mu0)

• IRD: mean(mu1) - mean(mu0)

Value

Numeric effect value, or NA_real_ if required columns are missing.

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 CDE Analogue from GLM Potential Outcomes

Description

Computes the Controlled Direct Effect analogue for GLM outcomes. For binary (logit link), this averages the individual-level odds first, then takes the ratio — analogous to mean(exp(theta_1)) / mean(exp(theta_0)) in survival:

Usage

compute_cde_glm(df, subset_indicator = NULL, effect_measure = "OR")

Arguments

Details

CDE(S) = mean(p1[S]/(1-p1[S])) / mean(p0[S]/(1-p0[S]))

This differs from the AOR (odds(mean(p1)) / odds(mean(p0))) due to Jensen's inequality whenever individual probabilities are heterogeneous. For identity-link (MD) and log-link (IRR) outcomes, CDE equals AOR and this function returns NA_real_.

Value

Numeric CDE value (binary OR only), or NA_real_.

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 for the full subgroup (d = total subgroup events); each 50/50 split-half (d/2 events) has variance 8/d

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

Value

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

Examples

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

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

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

## End(Not run)

Compute Detection Probability for GLM Outcomes

Description

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

Usage

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

Arguments

Details

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

1. Screening: W_1 + W_2 \ge 2k_1

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

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

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

Value

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

See Also

compute_detection_probability for the survival-specific version.

Examples

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

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

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

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

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

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

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

## End(Not run)

Compute Naive and Split-Derived 95% CIs for Pareto Frontier Members

Description

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

Usage

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

Arguments

Details

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

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

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

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

Value

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

m

Candidate id (matches pareto_frontier$m).

estimate

Refit effect estimate on the natural scale.

naive_lcl, naive_ucl

Naive Wald CI bounds.

split_sd

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

split_lcl, split_ucl

Split-derived CI bounds.

fsbc_estimate

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

fsbc_lcl, fsbc_ucl

FSBC-mimic CI bounds.

n_valid_splits

Number of splits that produced valid per-half estimates.

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

Details — Half-jackknife SE

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

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

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

Details — FSBC-mimic interpretation

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

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

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

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

See Also

pareto_frontier_table, frontier_member_flags, plot_pareto_frontier.

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

Examples

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

## End(Not run)

Compute Pareto Frontier on (Effect, N)

Description

Returns the subset of candidate subgroups that are not dominated on the two-objective space (effect size, sample size), where both objectives are maximized. A candidate i is dominated iff another candidate j has hr_j \ge hr_i and N_j \ge N_i, with at least one inequality strict.

Usage

compute_pareto_frontier(result_dt, effect_log_scale = FALSE)

Arguments

Details

The frontier is intended as a post-hoc reporting artifact, not a selection criterion. The selected subgroup (under sg_focus) may or may not appear on the frontier – in particular, "hrMinSG" may select an N-dominated point by design (preferring small subgroups).

Value

A data.table of non-dominated rows, sorted by hr descending. Returns an empty 0-row data.table if input is empty.

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

Wald confidence intervals for DINA coefficients.

Description

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

Usage

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

Arguments

Value

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

Resampling approximation to the splitting consistency rate

Description

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

Usage

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

Arguments

Details

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

Value

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

See Also

run_single_consistency_split(), consistency_resample_compare()

Examples

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

## End(Not run)

Validate the consistency approximation against literal splitting (Cox)

Description

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

Usage

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

Arguments

Value

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

See Also

consistency_resample(), run_single_consistency_split()

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.

Examples

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

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 = "plots",
  verbose = TRUE
)

Arguments

Value

List of class "cox_ahr_cde" containing:

cox_fit

Results from cox_cs_fit function.

ahr_results

AHR calculations for different subgroup definitions.

cde_results

CDE calculations if theta variables available.

optimal_cutpoint

Optimal biomarker cutpoint, or NULL when hr_threshold is NULL.

subgroup_stats

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

data

List with z_values, loghr_po, and subgroup assignments.

Examples

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

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

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

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

## End(Not run)

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

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

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

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

## End(Not run)

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

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

## End(Not run)

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

Description

Called in analyze_subgroup() <– SG_tab_estimates

Usage

cox_summary_legacy(Y, E, Treat, Strata)

Arguments

Details

Calculates hazard ratio and confidence interval for a subgroup using Cox regression.

Value

Character string with formatted HR and CI.

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,
  effect_label = "HR"
)

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

Examples

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

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

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

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

## End(Not run)

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

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

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

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

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

render_forestplot(result)

## End(Not run)

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

Examples

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

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

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

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

## End(Not run)

Compute GLM Effect Estimate for a Single Subgroup

Description

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

Usage

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

Arguments

Value

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

See Also

grf.subg.harm.glm

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

Examples

cfg <- create_grf_config(frac.tau = 0.6, n.min = 60, dmin.grf = 6,
                         RCT = TRUE, sg.criterion = "mDiff",
                         maxdepth = 2, seedit = 42L)
str(cfg)

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

Examples

## Not run: 
ss_table <- create_sample_size_table(
  theta_values = c(1.5, 1.75, 2.0, 2.5),
  prop_cens_values = c(0.2, 0.3, 0.4),
  target_power = 0.80
)
print(ss_table)

## End(Not run)

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)

Examples

## Not run: 
# Basic usage
create_summary_table(
  data = trial_data,
  treat_var = "treatment",
  vars_continuous = c("age", "bmi"),
  vars_categorical = c("sex", "stage")
)

# Customized appearance
create_summary_table(
  data = trial_data,
  treat_var = "treatment",
  vars_continuous = c("age", "bmi"),
  vars_categorical = c("sex", "stage"),
  font_size = 11,
  header_font_size = 13,
  use_alternating_rows = TRUE,
  highlight_pval = 0.05,
  compact_mode = TRUE
)

## End(Not run)

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.

Generate J-quantile cut expressions for a continuous variable

Description

For a continuous variable, returns J cut expressions of the form "X <= qj(X, k, J + 1)" for k = 1, ..., J. These J cut points are placed at the (k/(J+1))-th empirical quantiles of X and partition its range into J + 1 non-overlapping intervals

[\min(X), c_1),\ [c_1, c_2),\ \ldots,\ [c_J, \max(X)]

where c_k is the (k/(J+1))-th quantile of X. Note the cut expressions themselves are nested half-spaces (each is a subset of the next); ForestSearch consumes them as binary candidate factors and combines them via intersection during the search.

Usage

cut_var_jq(x, J)

Arguments

Details

Expressions are emitted in deferred form (with literal qj(...) calls inside the string) so that they are correctly recomputed when the same expression is processed against a different data subset (e.g., a bootstrap replicate). They are subsequently resolved to literal numerics by process_conf_force_expr().

Value

Character vector of J 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.

Examples

## Not run: 
# Compare CV results from different configurations
cv_comparison <- cv_compare_results(
  cv_list = list(
    "maxk=1" = cv_maxk1,
    "maxk=2" = cv_maxk2,
    "maxk=3" = cv_maxk3
  )
)

## End(Not run)

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

Examples

## Not run: 
# After running forestsearch_tenfold
tenfold_results <- forestsearch_tenfold(
  fs.est = fs_result,
  sims = 100,
  Kfolds = 10
)

# Create combined metrics table
cv_tables <- cv_metrics_tables(tenfold_results)
cv_tables

# Create separate tables
cv_tables <- cv_metrics_tables(tenfold_results, table_style = "separate")
cv_tables$agreement_table
cv_tables$finding_table

# Minimal one-row summary
cv_metrics_tables(tenfold_results, table_style = "minimal")

## End(Not run)

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

Examples

## Not run: 
# Run K-fold CV
cv_results <- forestsearch_Kfold(fs.est = fs_result, Kfolds = 10)

# Get detailed output
kfold_out <- forestsearch_KfoldOut(cv_results, outall = TRUE)

# Create summary tables
cv_tables <- cv_summary_tables(kfold_out)
cv_tables$combined_table
cv_tables$metrics_table

## End(Not run)

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.

Examples

## Not run: 
cv_text <- cv_summary_text(tenfold_results)
# Returns: "CV found = 95%, Agree(+,-) = 88%, 92%"

## End(Not run)

Effective Information for Binary Outcomes (Logistic)

Description

Computes d_{\text{eff}} = n_H \bar{p} (1 - \bar{p}) from the Fisher information of the logistic regression treatment coefficient.

Usage

d_eff_binary(n_sg, p_event)

Arguments

Details

For the model \text{logit}(p_i) = \beta_0 + \beta V_i with balanced arms (n_1 = n_0 = n_H/2):

\text{Var}(\hat\beta) = \frac{4}{n_H \bar{p}(1-\bar{p})}

Value

Numeric. Effective information.

Examples

d_eff_binary(n_sg = 100, p_event = 0.30)
# 21

Effective Information for Continuous Outcomes (Gaussian)

Description

Computes d_{\text{eff}} = n_H / \sigma_Y^2 from the Fisher information of the Gaussian linear model treatment coefficient.

Usage

d_eff_continuous(n_sg, sigma_y)

Arguments

Details

For the model Y_i = \beta_0 + \beta V_i + \epsilon_i, \epsilon_i \sim N(0, \sigma^2), with balanced arms:

\text{Var}(\hat\beta) = \frac{4\sigma^2}{n_H}

Value

Numeric. Effective information.

Examples

d_eff_continuous(n_sg = 100, sigma_y = 1.5)
# 44.4

Effective Information for Count Outcomes (Poisson)

Description

Computes d_{\text{eff}} = D / \phi where D is the total expected events and \phi is the dispersion parameter.

Usage

d_eff_count(
  n_sg = NULL,
  event_rate = NULL,
  total_events = NULL,
  dispersion = 1
)

Arguments

Details

For Poisson regression with log-offset:

\text{Var}(\hat\beta) = \phi \left(\frac{1}{D_0} + \frac{1}{D_1}\right) \approx \frac{4\phi}{D}

Under proportional hazards, D = d (Cox events), so d_{\text{eff}}^{\text{Poisson}} = d_{\text{eff}}^{\text{Cox}} (Laird and Olivier, 1981).

Value

Numeric. Effective information.

Examples

# From total events
d_eff_count(total_events = 55)

# From sample size and rate
d_eff_count(n_sg = 100, event_rate = 0.8)

# Quasi-Poisson with overdispersion
d_eff_count(total_events = 55, dispersion = 1.5)

Effective Information for Survival (Cox PH)

Description

Computes d_{\text{eff}} = n_H (1 - p_c) where p_c is the censoring proportion.

Usage

d_eff_survival(n_sg, prop_cens)

Arguments

Value

Numeric. Effective information (expected events).

Examples

d_eff_survival(n_sg = 60, prop_cens = 0.45)
# 33

Default GRF parameters (general)

Description

Default GRF parameters (general)

Usage

default_grf_params_gen()

Default ForestSearch parameters (general)

Description

Returns the default fs_params list used by run_simulation_analysis(). The parallel_args entry defaults to list(plan = "sequential") because run_simulation_analysis() is designed to be called inside a foreach() %dofuture% loop (one replicate per worker). Running the inner forestsearch() multisession in that context produces nested parallelism: each outer worker tries to spawn its own pool of inner workers, which parallelly rejects with a 300\ limit error. Users who call run_simulation_analysis() once at the top level (i.e., not inside %dofuture%) can opt back into multisession by passing fs_params = list(parallel_args = list(plan = "multisession", workers = N)).

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 Leon et al. (2024, Section 2.1) approximation, each 50/50 split-half log(HR) estimator follows N(log(theta), 8/d) where d = n_sg * (1 - prop_cens) is the expected TOTAL number of events in the subgroup. Each split holds d/2 events, so its variance is 4/(d/2) = 8/d.

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:

Examples

## Not run: 
example_data <- data.frame(
  id = 1:100,
  age = rnorm(100, 50, 10),
  grade = sample(1:3, 100, replace = TRUE),
  status = sample(c("Active", "Inactive"), 100, replace = TRUE),
  score = runif(100, 0, 100)
)
result <- detect_variable_types(example_data,
                                 max_unique_for_cat = 10,
                                 exclude_vars = "id")
result$continuous_vars  # c("age", "score")
result$cat_vars         # c("grade", "status")
names(result$data_subset)  # c("age", "grade", "status", "score")

## End(Not run)

Estimate heterogeneous treatment effects via DINA (data-frame interface)

Description

dina() is the high-level user-facing wrapper around dina_fit(). It accepts a data frame plus column-name arguments, matching the conventions used elsewhere in the forestsearch package, and extracts the covariate matrix, response, and treatment indicator before dispatching to dina_fit().

Usage

dina(
  df,
  outcome,
  treatment,
  covariates,
  family = c("gaussian", "binomial", "poisson", "cox"),
  status = NULL,
  ...
)

Arguments

Details

All other arguments are forwarded to dina_fit() via ..., including propensity_method, baseline_method, cross_fitting, n_folds, cens_type, cens_params, n_grid, eps, and seed. See dina_fit() for their full documentation.

Covariate columns must be numeric (or integer). Convert factors to dummy variables (e.g., with stats::model.matrix()) before calling if needed. Rows with any NA value in the referenced columns are rejected with an error; use stats::na.omit() or impute beforehand.

Value

An object of class "dina", as returned by dina_fit(). The ⁠$call⁠ component is overwritten with the original dina() call so that print() and summary() show the high-level invocation rather than the internal matrix-based dispatch.

See Also

dina_fit() for the underlying matrix-based estimator; dina_bagged() for the bagged version with infinitesimal-jackknife variance.

Examples

set.seed(1)
n <- 400
df_demo <- data.frame(
  y         = NA,
  w         = stats::rbinom(n, 1, 0.5),
  x1        = stats::runif(n, -1, 1),
  x2        = stats::runif(n, -1, 1),
  x3        = stats::runif(n, -1, 1)
)
tau_x      <- 0.5 + 0.8 * df_demo$x1 - 0.3 * df_demo$x2
df_demo$y  <- 1 + df_demo$x1 + df_demo$w * tau_x + stats::rnorm(n)

fit <- dina(
  df         = df_demo,
  outcome    = "y",
  treatment  = "w",
  covariates = c("x1", "x2", "x3"),
  family     = "gaussian",
  seed       = 1L
)
coef(fit)
confint(fit)

Bagged DINA with IJ variance (data-frame interface)

Description

dina_bagged() is the high-level user-facing wrapper around dina_fit_bagged(). It accepts a data frame plus column-name arguments and dispatches to dina_fit_bagged(), which fits the DINA estimator on n_bags bootstrap replicates and returns the infinitesimal-jackknife variance per Wager, Hastie and Efron (2014).

Usage

dina_bagged(
  df,
  outcome,
  treatment,
  covariates,
  family = c("gaussian", "binomial", "poisson", "cox"),
  status = NULL,
  ...
)

Arguments

Details

All bagging-related arguments (n_bags, cross_fitting_per_bag, n_folds_per_bag, parallel, ij_finite_sample_correction, project_psd, verbose, ...) are forwarded via .... See dina_fit_bagged() for full documentation.

For parallel = "bags", set future::plan() to a parallel backend before calling dina_bagged(). Restore with ⁠future::plan("sequential"); gc()⁠ afterwards.

Value

An object of class c("dina_bagged", "dina"), as returned by dina_fit_bagged(). The ⁠$call⁠ component is overwritten with the original dina_bagged() call.

See Also

dina_fit_bagged() for the underlying matrix-based bagged estimator; dina() for the single-pass DINA with sandwich variance.

Examples

set.seed(1)
n <- 400
df_demo <- data.frame(
  y         = NA,
  w         = stats::rbinom(n, 1, 0.5),
  x1        = stats::runif(n, -1, 1),
  x2        = stats::runif(n, -1, 1),
  x3        = stats::runif(n, -1, 1)
)
tau_x      <- 0.5 + 0.8 * df_demo$x1 - 0.3 * df_demo$x2
df_demo$y  <- 1 + df_demo$x1 + df_demo$w * tau_x + stats::rnorm(n)

fit <- dina_bagged(
  df         = df_demo,
  outcome    = "y",
  treatment  = "w",
  covariates = c("x1", "x2", "x3"),
  family     = "gaussian",
  n_bags     = 30L,
  seed       = 1L
)
coef(fit)
sqrt(diag(vcov(fit)))

Estimate heterogeneous treatment effects via DINA

Description

dina_fit() implements the two-step DINA (DIfference in NAtural parameters) estimator of Gao and Hastie (2025) for heterogeneous treatment effects (HTE) with general response types: Gaussian (mean difference), binomial (log odds ratio), Poisson (log mean ratio), and Cox (log hazard ratio). The treatment effect is assumed to be linear in the covariates,

\tau(x) = \beta_0 + x^\top \beta_{1:d},

following Proposition 1 of the paper. The procedure has two steps: (1) estimate nuisance functions a(x) and \nu(x) on a training fold, and (2) plug those estimators into a GLM (or Cox) likelihood on the held-out fold to estimate \beta. Cross-fitting averages the two estimators obtained by swapping the folds.

Usage

dina_fit(
  X,
  Y,
  W,
  family = c("gaussian", "binomial", "poisson", "cox"),
  propensity_method = "logistic",
  baseline_method = NULL,
  cross_fitting = TRUE,
  n_folds = 2L,
  cens_type = NULL,
  cens_params = list(),
  n_grid = 500L,
  eps = 1e-06,
  seed = NULL
)

Arguments