Getting started with forestsearch

knitr::opts_chunk$set(comment = "#>", fig.align = "center")
library(forestsearch)
library(survival)

A confirmatory trial is powered for an effect in its whole population, and the indication it supports describes a population in which efficacy and a favourable benefit–risk balance have been concluded. Regulatory guidance asks that this conclusion be examined across relevant subgroups rather than assumed, and that a sub-population in which it does not hold be left out of the indication; it credits a subgroup finding only when credible — tied to prior plausibility and replication — and warns that an estimate singled out for its extremity is biased by the selection (European Medicines Agency 2019). Approvals have been restricted to a subgroup after a positive result in the whole population (Amatya et al. 2021).

A subgroup of concern is reached by searching the trial’s covariates after the data are seen, or by examining many prespecified subgroups and singling one out. Either way it is chosen for the extremity of its estimate, and its effect is reported through the analysis the trial would have run had that subgroup been the only one: the treatment coefficient of the trial’s working model fitted within it, a Cox log hazard ratio for a time-to-event outcome. So chosen, the subgroup carries an estimate that overstates its effect, and the same selection distorts the estimate for the patients outside it. forestsearch is built for this situation. It does two things.

• Identification. It searches for a subgroup in the direction of a sub-optimal effect — the package’s “harm” direction, a hazard ratio above one for a survival outcome. Every candidate is a conjunction of at most two covariate cuts, so whether a patient belongs is read from the patient’s own measurements. Three identifiers share that form: forest search (León et al. 2024), the causal forest of generalized random forests (GRF; Athey et al. 2019, Cui et al. 2023) and the difference-in-natural-parameters estimator (DINA; Gao and Hastie 2025).
• Inference after selection. Refit-free multiplier resampling over the candidate family removes the optimism of selection from the estimate in the identified region and in its complement, and gives one-sided bounds on both (León and Anderson 2026a, 2026b).

This vignette runs the identification, by forest search, on the German Breast Cancer Study Group (GBSG) trial that ships with survival, at the settings of León and Anderson (2026b), and shows the call that adds the inference. It takes a few seconds on one core.

The complete analysis — the three identifiers, the selection-adjusted bounds for each identified region and its complement, and the declaration calibration, which together are Table 1 of León and Anderson (2026b) run end to end — is the article Interpretable subgroup identification with forestsearch, listed under Articles on the package website.

The data

The analysis set is the 686 patients with complete prognostic data, out of 720 randomized: 246 received hormonal therapy and 440 did not, with 299 recurrences or deaths. Seven covariates enter the search: age, menopausal status, tumour size, grade, the number of positive lymph nodes, and the progesterone and estrogen receptor concentrations in fmol/l. Randomization status is masked in the distributed data; hormonal therapy is analysed as randomized.

df <- survival::gbsg
df <- within(df, {
  id          <- seq_len(nrow(df))              # the fit's subject identifier
  time_months <- rfstime / 30.4375              # recurrence-free survival, days to months
  grade3      <- ifelse(grade == "3", 1, 0)     # grade 3 against grade 1 or 2
})
confounders <- c("age", "meno", "size", "grade3", "nodes", "pgr", "er")
stopifnot(nrow(df) == 686, sum(df$status) == 299)

whole <- coxph(Surv(time_months, status) ~ hormon, data = df)
whole_hr <- c(est = unname(exp(coef(whole))), exp(confint(whole))[1, ])
c(patients = nrow(df), hormonal_therapy = sum(df$hormon == 1),
  none = sum(df$hormon == 0), events = sum(df$status))

#>         patients hormonal_therapy             none           events 
#>              686              246              440              299

sprintf("Whole-trial hazard ratio %.3f (%.3f, %.3f), hormonal therapy vs none",
        whole_hr[["est"]], whole_hr[[2]], whole_hr[[3]])

#> [1] "Whole-trial hazard ratio 0.695 (0.544, 0.888), hormonal therapy vs none"

Every estimate below is read against the trial’s own effect, a hazard ratio of 0.695, and against clinically meaningful effect sizes — not as a test at a hazard ratio of one.

The specification

Candidates are conjunctions of at most two covariate cuts (maxk = 2) with more than 60 patients (n.min) and ten events in each arm (d0.min, d1.min). Continuous covariates take the package’s default cuts, except that the two receptors carry ten candidate cuts each (conf.cont_jcuts) with the value zero — no measured receptor — forced in as a cut of its own (conf_force).

Three thresholds govern the search, on the hazard-ratio scale with harm upward:

• the screening threshold c₁ = hr.threshold, which a candidate’s estimated hazard ratio must clear to be admitted;
• the per-split threshold c₂ = hr.consistency, the claim threshold at which an effect stops being beneficial, which both halves of a random half-split must clear for the split to count as consistent;
• the consistency level p^⋆ = pconsistency.threshold, the proportion of consistent splits a candidate must reach.

León and Anderson (2026b) set c₁ = c₂ = 1.0 and p^⋆ = 0.90. The splitting-consistency rate is not computed by drawing splits: with consistency_method = "resample" (the default) it is evaluated in closed form from one fit of the working model within the candidate, so nothing is drawn and the search is deterministic given the data. The article gives the formula and the fallback used when the closed form cannot be evaluated.

The selection rule is max_𝒜N(ε): among the admitted candidates 𝒜, the largest whose hazard ratio is within a fraction ε = 0.20 of the largest hazard ratio over the family (sg_focus = "effMaxSG", selection_rule = "neighborhood", effect_neighborhood = 0.20). The rule is a modelling choice: it moves the returned region along a frontier of size against effect, and so determines what identification and inference are about.

Forest search

One call. parallel_args = list(plan = "sequential") keeps the run on one core; the package’s default plan is parallel.

t0 <- proc.time()[["elapsed"]]
fit_fs <- forestsearch(
  df.analysis = df,
  outcome.name = "time_months", event.name = "status", treat.name = "hormon", id.name = "id",
  confounders.name = confounders,
  is.RCT = TRUE, est.scale = "hr", seedit = 8316951,
  # the candidate family
  conf.cont_jcuts = list(er = 10, pgr = 10), conf_force = c("er <= 0", "pgr <= 0"),
  n.min = 60, d0.min = 10, d1.min = 10, maxk = 2,
  # the identifier, its thresholds and the selection rule
  subgroup_method = "consistency", consistency_method = "resample",
  hr.threshold = 1.00, hr.consistency = 1.00, pconsistency.threshold = 0.90,
  stop_threshold = NULL,                      # evaluate every candidate; no early stop
  sg_focus = "effMaxSG", selection_rule = "neighborhood", effect_neighborhood = 0.20,
  quiet = TRUE, parallel_args = list(plan = "sequential"))
fit_sec_fs <- proc.time()[["elapsed"]] - t0

in_H <- fit_fs$df.est$treat.recommend[match(df$id, fit_fs$df.est$id)] == 0
stopifnot(!anyNA(in_H))
cat(sprintf("Identified subgroup: %s\n", paste(fit_fs$sg.harm, collapse = " & ")))

#> Identified subgroup: {er <= 0} & {pgr <= 26}

cat(sprintf("Patients in it: %d of %d. Fit: %.1f s\n", sum(in_H), nrow(df), fit_sec_fs))

#> Patients in it: 75 of 686. Fit: 1.9 s

The identified region is {er <= 0} & {pgr <= 26}: no measured estrogen receptor and progesterone receptor at most 26 fmol/l. One argument in the call is not the default and is there for a reason: stop_threshold = NULL has the consistency stage evaluate every screened candidate rather than stop at the first that passes. Everything else that is not the default is the specification itself.

The search in numbers

The fit records how the candidates were whittled down.

counts <- fit_fs$find.grps$filter_counts
funnel <- data.frame(
  stage = c("enumerated (at most two cuts, size-eligible)",
            "fittable (ten events in each arm, Cox fit obtained)",
            "clearing the screening threshold, HR > 1.0",
            "passing the consistency screen, p_cons >= 0.90"),
  candidates = c(counts$n_evaluated, counts$n_passed_cox, counts$n_passed_hr,
                 fit_fs$grp.consistency$n_passed))
knitr::kable(funnel, row.names = FALSE, format.args = list(big.mark = ","))

The candidates that pass the screen, with their size, hazard ratio and consistency rate, are the points on which the selection rule acts.

cand <- as.data.frame(fit_fs$grp.consistency$out_sg$result)
cand$candidate <- ifelse(cand$M.2 == "", cand$M.1, paste(cand$M.1, cand$M.2, sep = " & "))
cand <- cand[order(-cand$N), ]
knitr::kable(cand[, c("candidate", "N", "hr", "Pcons")], digits = 3,
             row.names = FALSE, col.names = c("candidate", "N", "HR", "p_cons"),
             caption = "The candidates passing the consistency screen, largest first.")

band <- (1 - 0.20) * max(cand$hr)
within_band <- cand$hr >= band
sel <- cand[within_band, ][which.max(cand$N[within_band]), ]
cat(sprintf("Largest HR over the family: %.2f (N %d). Band: HR >= %.2f.\n",
            max(cand$hr), cand$N[which.max(cand$hr)], band))

#> Largest HR over the family: 2.54 (N 61). Band: HR >= 2.03.

cat(sprintf("max N(eps) takes the largest candidate within the band: N %d, HR %.2f.\n",
            sel$N, sel$hr))

#> max N(eps) takes the largest candidate within the band: N 75, HR 2.22.

cat(sprintf("max N alone would take N %d (HR %.2f); max beta the N %d.\n",
            max(cand$N), cand$hr[which.max(cand$N)], cand$N[which.max(cand$hr)]))

#> max N alone would take N 82 (HR 1.95); max beta the N 61.

stopifnot(sel$N == sum(in_H))   # the rule, re-applied, picks the fit's own selection

All of the passing candidates restrict the estrogen receptor, and the rules differ only in how far that region is trimmed. The max_𝒜N(ε) rule takes the 75 patients with no measured estrogen receptor and progesterone receptor at most 26 fmol/l; max_𝒜N would take the 82 with no measured estrogen receptor; the effect-maximising rules the 61 with, in addition, tumour size at most 35 mm.

The region and its complement, before correction

The standard analysis within the identified region and within its complement: the two-sample Cox fit of treatment, with its 95% Wald interval.

hr_ci <- function(d) {
  f <- coxph(Surv(time_months, status) ~ hormon, data = d)
  ci <- exp(confint(f))
  sprintf("%.3f (%.3f, %.3f)", exp(coef(f)), ci[1, 1], ci[1, 2])
}
knitr::kable(
  data.frame(set = c("identified region", "complement"),
             N = c(sum(in_H), sum(!in_H)),
             events = c(sum(df$status[in_H]), sum(df$status[!in_H])),
             unadjusted = c(hr_ci(df[in_H, ]), hr_ci(df[!in_H, ]))),
  col.names = c("", "N", "events", "unadjusted HR (95% CI)"), row.names = FALSE,
  caption = "The identified region and its complement: unadjusted estimates, before any correction for the selection.")

These are not the estimates to report. The region was chosen for the extremity of its estimate, so its hazard ratio overstates its effect, and the same selection distorts the complement’s. They are shown as what the inference below corrects.

Inference after selection

Adding mr_inference = TRUE to the same call has the fit correct both estimates and bound them:

• the correction — refit-free multiplier resampling over the candidate family removes the optimism of selection from the estimate in the identified region Ĥ and in its complement Ĥ^c (León and Anderson 2026a);
• one-sided bounds — the field lower bound on the effect in Ĥ, the field-s upper bound on the effect in Ĥ^c, relevant to a benefit claim for the patients outside the subgroup, and the Bonferroni pair, which joins the two in one claim at a joint level (León and Anderson 2026b).

The call is shown, not run, to keep this vignette short. The article runs it, for all three identifiers, prints its fit times and reads the results.

fit_fs_inf <- forestsearch(
  df.analysis = df,
  outcome.name = "time_months", event.name = "status", treat.name = "hormon", id.name = "id",
  confounders.name = confounders,
  is.RCT = TRUE, est.scale = "hr", seedit = 8316951,
  conf.cont_jcuts = list(er = 10, pgr = 10), conf_force = c("er <= 0", "pgr <= 0"),
  n.min = 60, d0.min = 10, d1.min = 10, maxk = 2,
  subgroup_method = "consistency", consistency_method = "resample",
  hr.threshold = 1.00, hr.consistency = 1.00, pconsistency.threshold = 0.90,
  stop_threshold = NULL,
  sg_focus = "effMaxSG", selection_rule = "neighborhood", effect_neighborhood = 0.20,
  # selection-adjusted inference
  mr_inference = TRUE,
  mr_inference_args = list(ci_method = "field", draws = 5000L, include_complement = TRUE,
                           field_scale_complement = "selected"),
  quiet = TRUE, parallel_args = list(plan = "sequential"))
mr_estimates_table(fit_fs_inf, est.scale = "hr", digits = 3)

Each estimate and bound is read by its location against the trial’s own effect and against clinically meaningful effect sizes, not as a test at the null. Whether and how to act on a finding is a broader judgement of benefit and risk, which this inference is meant to inform.

Where to go next

• The article Interpretable subgroup identification with forestsearch, listed under Articles on the package website, continues this analysis: the selection-adjusted bounds for this region and its complement, the two model-based identifiers (DINA and GRF), Table 1 of León and Anderson (2026b), and the declaration calibration at a protected level.
• The post-selection article reads the bounds in more depth on the same trial: the re-selection frequency p̂ as a diagnostic, membership agreement, and what each construction does and does not say.
• The Guo and He article sets the bounds beside the de-biased estimate and one-sided bound of Guo and He (2021) on a nested continuum of subgroups, Section S4 of León and Anderson (2026b).
• The ACTG 175 analyses are an article on the package website: forest search on a binary and a continuous endpoint, with the bounds and the declaration calibration (León and Anderson 2026b, Section S3), and the three identifiers on the binary endpoint with the full bootstrap beside multiplier resampling (León and Anderson 2026a, Table S3).
• The package website carries the original workflow of León et al. (2024) — the full argument listing, the bootstrap bias correction and the cross-validation — as the article Exploratory subgroup identification with forestsearch (the 2024 workflow), and the methodology and simulation articles.
• ?forestsearch documents every argument, with examples on this trial and on the binary and continuous ACTG 175 endpoints, and a section on the options not tested against the manuscript analyses.

References

Amatya, A. K., Fiero, M. H., Bloomquist, E. W., Sinha, A. K., Lemery, S. J., Singh, H., Ibrahim, A., Donoghue, M., Fashoyin-Aje, L. A., de Claro, R. A., Gormley, N. J., Amiri-Kordestani, L., Sridhara, R., Theoret, M. R., Kluetz, P. G., Pazdur, R., Beaver, J. A. and Tang, S. (2021). Subgroup analyses in oncology trials: Regulatory considerations and case examples. Clinical Cancer Research 27(21), 5753–5756.

Athey, S., Tibshirani, J. and Wager, S. (2019). Generalized random forests. The Annals of Statistics 47(2), 1148–1178.

Cui, Y., Kosorok, M. R., Sverdrup, E., Wager, S. and Zhu, R. (2023). Estimating heterogeneous treatment effects with right-censored data via causal survival forests. Journal of the Royal Statistical Society Series B 85(2), 179–211.

European Medicines Agency (2019). Guideline on the investigation of subgroups in confirmatory clinical trials. EMA/CHMP/539146/2013, Committee for Medicinal Products for Human Use.

Gao, Z. and Hastie, T. (2025). Estimating heterogeneous treatment effects for general responses. Biometrics 81(4), ujaf162.

León, L. F. and Anderson, K. M. (2026a). Inference for standard trial estimands after data-driven subgroup discovery. arXiv:2609.38361. https://doi.org/10.48550/arXiv.2609.38361

León, L. F. and Anderson, K. M. (2026b). Interpretable subgroup identification across survival and generalized-linear outcomes. arXiv:2610.28379. https://doi.org/10.48550/arXiv.2610.28379

León, L. F., Jemielita, T., Guo, Z., Marceau West, R. and Anderson, K. M. (2024). Exploratory subgroup identification in the heterogeneous Cox model: A relatively simple procedure. Statistics in Medicine. https://doi.org/10.1002/sim.10163

sessionInfo()

#> R version 4.6.2 (2026-10-29)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Pop!_OS 24.04 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/Los_Angeles
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] future_1.80.0      survival_3.8-12    forestsearch_0.5.0
#> 
#> loaded via a namespace (and not attached):
#>  [1] gt_1.4.0             generics_0.1.4       xml2_1.6.0          
#>  [4] shape_1.4.6.1        lattice_0.23-1       listenv_1.1.0       
#>  [7] digest_0.6.39        magrittr_2.0.5       grf_2.6.1           
#> [10] evaluate_1.0.5       grid_4.6.2           RColorBrewer_1.1-3  
#> [13] iterators_1.0.14     policytree_1.2.5     mvtnorm_1.4-2       
#> [16] fastmap_1.2.0        foreach_1.5.2        jsonlite_2.0.0      
#> [19] Matrix_1.7-6         glmnet_5.1           scales_1.4.0        
#> [22] weightedsurv_0.1.0   codetools_0.2-20     cli_3.6.6           
#> [25] rlang_1.3.0          parallelly_1.48.0    future.apply_1.21.0 
#> [28] splines_4.6.2        yaml_2.3.12          otel_0.2.0          
#> [31] tools_4.6.2          parallel_4.6.2       doFuture_1.4.0      
#> [34] dplyr_1.2.1          ggplot2_4.0.3        globals_0.20.0      
#> [37] vctrs_0.7.3          R6_2.6.1             lifecycle_1.0.5     
#> [40] randomForest_4.7-1.2 fs_2.1.0             pkgconfig_2.0.3     
#> [43] progressr_1.1.0      pillar_1.11.1        gtable_0.3.6        
#> [46] data.table_1.18.6.1  glue_1.8.1           Rcpp_1.1.2          
#> [49] xfun_0.62            tibble_3.3.1         tidyselect_1.2.1    
#> [52] rstudioapi_0.19.0    knitr_1.52           farver_2.1.2        
#> [55] patchwork_1.3.2      htmltools_0.5.9      rmarkdown_2.32      
#> [58] compiler_4.6.2       S7_0.2.2