--- title: "Interpretable subgroup identification with forestsearch" subtitle: "The German Breast Cancer Study Group analysis of León and Anderson (2026b), Table 1, run end to end" author: "Larry Leon" date: today format: html: embed-resources: true toc: true toc-depth: 3 code-fold: false code-overflow: wrap fontsize: 15px vignette: > %\VignetteIndexEntry{Interpretable subgroup identification with forestsearch} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- ```{r setup} #| message: false #| warning: false knitr::opts_chunk$set(comment = "#>", fig.align = "center") library(forestsearch) library(survival) ``` ```{=html} ``` 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. **Identification — what the package searches for, and how** - **Interpretable candidates.** Every candidate subgroup is a conjunction of at most two covariate cuts, so whether a patient belongs is read from the patient's own measurements. - **Three identifiers, one 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) each score a family of candidates by an estimated treatment effect, admit those whose standard-analysis estimate $\hat\beta(g)$ clears a threshold common to all three, and return one admitted candidate as a selection rule directs — or none, if none is admitted. Forest search enumerates the candidates exhaustively; DINA and GRF propose them from a surface fitted to the outcome, and every proposal is re-ranked on $\hat\beta(g)$ against the same threshold. - **The orientation.** The search runs in the direction of a sub-optimal effect — the package's "harm" direction, a hazard ratio above one for a survival outcome — the direction in which the trial's conclusion of benefit would not hold for the patients concerned. - **The consistency screen, in closed form.** Forest search additionally requires a candidate's effect to replicate across random half-splits. That screen is complementary-pairs stability selection applied to a smooth estimating-function coefficient, so the replication rate follows in closed form from one fit of the working model within the candidate: nothing is drawn, and the search is deterministic given the data. - **The selection rule is a modelling choice.** The rule moves the returned region along a frontier of size against effect — the largest admitted region, the largest within a band of the strongest effect, the strongest effect itself — and so determines what identification and inference are about. - **The targets.** The region returned is $\hat H$; the targets are $\beta(\hat H)$ and, for the patients outside it, $\beta(\hat H^{c})$: the working model's coefficient in each population, its least-false value whether or not the model is correct, on the scale on which the trial was powered and on which regulators and clinicians read its results. **Inference — what is reported about the identified region, and about the rest of the trial** - **The correction.** Refit-free multiplier resampling over the candidate family removes the optimism of selection from the estimate in $\hat H$ and in $\hat H^{c}$ (León and Anderson 2026a). - **One-sided bounds, and a claim about both groups.** The *field* inverts the law of the residual the correction leaves into a one-sided lower bound on $\beta(\hat H)$; *field-s* gives the one-sided upper bound on $\beta(\hat H^{c})$, relevant to a benefit claim for the patients outside the subgroup; the Bonferroni pair joins the two in a two-subgroup claim at a joint level. Coverage is asymptotically exact under separated selection, characterized under local competition, and conditional on the family proposed when the family comes from the outcomes (León and Anderson 2026b). - **Declaration control.** Whether a subgroup should be declared at all is governed separately: the law of the maximum over the enumerated family gives the family-wise rate of falsely declaring a subgroup when every candidate's effect is at least a clinically specified benefit — a diagnostic of the rate the chosen thresholds imply and, where a capped rate is wanted, the calibrated cutoff that holds it (León and Anderson 2026b). - **Reading.** 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. **This vignette** runs the analysis reported as Table 1 of León and Anderson (2026b), on the German Breast Cancer Study Group (GBSG) trial that ships with **survival**, at the manuscript's settings: the three identifiers, each by one call; the selection-adjusted analysis of each identified region and its complement; and the declaration calibration at a protected level. The computation takes about a minute on one core; the fit times are printed as they are measured. ## 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. ```{r data} 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)) sprintf("Whole-trial hazard ratio %.3f (%.3f, %.3f), hormonal therapy vs none", whole_hr[["est"]], whole_hr[[2]], whole_hr[[3]]) ``` Every bound below is read against the trial's own effect, a hazard ratio of `r sprintf("%.3f", whole_hr[["est"]])`, 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`); forest search additionally requires 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_1$ = `hr.threshold`, which a candidate's estimated hazard ratio must clear to be admitted; - the **per-split threshold** $c_2$ = `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^\star$ = `pconsistency.threshold`, the proportion of consistent splits a candidate must reach. The applications of León and Anderson (2026b) set $c_1 = c_2 = 1.0$ and $p^\star = 0.90$. The splitting-consistency rate is *not* computed by drawing splits. With `consistency_method = "resample"` (the default) the rate is evaluated in closed form, as Section 2.2 of that paper gives it: to first order the two half-sample estimates are a mirror pair about the full-sample fit, so the probability that both clear $c_2$ is $$ p_{\mathrm{cons}}(g) \;=\; 2\,\Phi\!\left(\frac{\hat\beta(g) - \log c_2}{\sigma_D(g)}\right) - 1, $$ with $\sigma_D(g)$ the scale of the half-sample difference, obtained from the fit's influence function. Despite its name, `"resample"` draws nothing. Literal random half-splitting (`consistency_method = "split"`, `fs.splits` refits per candidate) is retained only as the fallback for a candidate whose closed form cannot be evaluated — a non-convergent fit, for instance — and `fs.splits` has no effect otherwise. Here every candidate has at least ten events in each arm and the Cox fits converge, so the closed form is evaluable throughout and the search is deterministic given the data. The rate is rounded to `pconsistency.digits = 2` decimals before it is compared with $p^\star$. The **selection rule** is $\max_{\mathcal{A}} N(\varepsilon)$: among the admitted candidates $\mathcal{A}$, the largest whose hazard ratio is within a fraction $\varepsilon = 0.20$ of the largest hazard ratio over the family (`sg_focus = "effMaxSG"`, `selection_rule = "neighborhood"`, `effect_neighborhood = 0.20`). Section 2.1 of León and Anderson (2026b) sets this rule beside the alternatives ($\max_{\mathcal{A}} N$, $\max_{\mathcal{A}} \hat\beta$, and their smallest-subgroup counterparts), which Figure 1 of that paper draws on this trial's candidate frontier. **Inference** (`mr_inference = TRUE`): from 5,000 multiplier draws over every size-eligible candidate the identifier proposes, the corrected (two-term) estimate of the effect in the identified region $\hat H$ and in its complement $\hat H^{c}$ (León and Anderson 2026a); the nominal 95% one-sided *field* lower bound on $\beta(\hat H)$, the *field-s* upper bound on $\beta(\hat H^{c})$, the studentized complement bound, and the Bonferroni pair, a joint 95% claim on both (León and Anderson 2026b). The field's own Monte Carlo runs at 1,000 outer and 500 inner draws, the package's fixed setting on this path. The only randomness in the analysis is in these draws; they are seeded by `seedit`, so the vignette's numbers are reproducible. ## Forest search One call. Identification, the unadjusted estimates and every selection-adjusted quantity come out of it together. `parallel_args = list(plan = "sequential")` keeps the run on one core; the package's default plan is parallel. ```{r fit-fs} #| message: false #| results: hide 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, # selection-adjusted inference mr_inference = TRUE, mr_inference_args = list(ci_method = "field", draws = 5000L, include_complement = TRUE, field_scale_complement = "selected", field_decompose = TRUE, keep_declaration_field = TRUE, declaration_c0 = c(0.70, 0.75, 0.80, 0.85)), quiet = TRUE, parallel_args = list(plan = "sequential")) fit_sec_fs <- proc.time()[["elapsed"]] - t0 ``` ```{r fs-selected} g_fs <- fit_fs$mr_inference cat(sprintf("Identified subgroup: %s\n", paste(fit_fs$sg.harm, collapse = " & "))) cat(sprintf("Fit including inference: %.1f s\n", fit_sec_fs)) cat(sprintf(" multiplier stage %.1f s; field blocks %.1f + %.1f s\n", g_fs$timing_seconds, g_fs$field$timing_seconds, g_fs$field$complement$timing_seconds)) ``` The identified region is `r sprintf("\x60%s\x60", paste(fit_fs$sg.harm, collapse = " & "))`: no measured estrogen receptor and progesterone receptor at most 26 fmol/l. Two arguments in the call are not the defaults and are there for a reason: `stop_threshold = NULL` has the consistency stage evaluate every screened candidate rather than stop at the first that passes, and `field_decompose = TRUE` returns the supporting quantities shown below. Everything else that is not the default is the specification itself. ### The search in numbers The fit records how the candidates were whittled down. Supplementary Table S13 of León and Anderson (2026b) reports the same funnel. ```{r fs-funnel} 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", "re-selection family for inference"), candidates = c(counts$n_evaluated, counts$n_passed_cox, counts$n_passed_hr, fit_fs$grp.consistency$n_passed, g_fs$n_family)) knitr::kable(funnel, row.names = FALSE, format.args = list(big.mark = ",")) ``` The re-selection family is larger than the screened set: multiplier resampling re-selects over every size-eligible candidate under the screening threshold and the consistency screen, without the per-arm event minimum, so the family contains the one the search ranked (Supplementary Section S5). The candidates that pass the screen, with their size, hazard ratio and consistency rate, are the points of Figure 1 of the paper. The selection rule acts on them. ```{r fs-candidates} 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)) cat(sprintf("max N(eps) takes the largest candidate within the band: N %d, HR %.2f.\n", sel$N, sel$hr)) 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)])) stopifnot(sel$N == g_fs$n_selected) # 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_{\mathcal{A}} N(\varepsilon)$ rule takes the 75 patients with no measured estrogen receptor and progesterone receptor at most 26 fmol/l; $\max_{\mathcal{A}} 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 selection-adjusted analysis of a region The nine quantities of Table 1 are read from `fit$mr_inference`. There is no single accessor for them, so the vignette defines one; `print(fit)` and `summary(fit)` report the bounds and the diagnostics in their own layout. ```{r table1-helper} table1_column <- function(fit, data) { g <- fit$mr_inference; cc <- g$complement fb <- g$field; fc <- fb$complement; js <- fb$joint_s in_H <- fit$df.est$treat.recommend[match(data$id, fit$df.est$id)] == 0 stopifnot(!anyNA(in_H), sum(in_H) == g$n_selected, sum(!in_H) == cc$n) f3 <- function(x) sprintf("%.3f", x) ci <- function(e, l, u) sprintf("%s (%s, %s)", f3(e), f3(l), f3(u)) c(paste(fit$sg.harm, collapse = " & "), sprintf("%d (%d)", sum(in_H), sum(data$status[in_H])), ci(g$naive$est, g$naive$lower, g$naive$upper), f3(g$debiased$est), f3(fb$lower_1s), sprintf("%d (%d)", sum(!in_H), sum(data$status[!in_H])), ci(cc$naive$est, cc$naive$lower, cc$naive$upper), f3(cc$debiased$est), f3(fc$upper_1s_s), sprintf("%s / %s", f3(js$bonf_lower_H), f3(js$bonf_upper_Hc))) } table1_rows <- c("Region", "N (events)", "Unadjusted, 95% CI", "Corrected", "Field, lower bound", "Complement N (events)", "Unadjusted, 95% CI", "Corrected", "Field-s, upper bound", "Bonferroni pair") knitr::kable(data.frame(row = table1_rows, FS = table1_column(fit_fs, df)), col.names = c("", "FS"), row.names = FALSE, caption = "Forest search: the identified region and its complement.") ``` The unadjusted hazard ratio in the region is opposite in direction to the trial's; the correction removes `r sprintf("%.3f", g_fs$selection_bias)` on the log scale and leaves `r sprintf("%.3f", g_fs$debiased$est)`. The field lower bound is the one-sided 95% bound on the region's effect; the field-s upper bound the one-sided 95% bound on the complement's; the Bonferroni pair is both at joint 95%, each side at 0.975. The quantities that support the bounds, as in Supplementary Table S18 of the paper: ```{r fs-supporting} fb <- g_fs$field; fc <- fb$complement; js <- fb$joint_s supporting <- data.frame( quantity = c("mean of the simulated selection root, region (log scale)", "second-order estimate, region", "second-order estimate, complement (field-s)", "scale ratio rho_c (field-s vs the unstudentized complement bound)", "Bonferroni pair: simulated joint coverage", "calibrated split gamma (Bonferroni's own alpha/2 is 0.025)", "correlation of the aligned roots"), value = c(sprintf("%.3f", fb$lambda_mean), sprintf("%.3f", fb$est2), sprintf("%.3f", fc$est2_s), sprintf("%.3f", fc$decomp_fields$scale_ratio_c), sprintf("%.4f", js$bonf_joint_prob), sprintf("%.3f", js$gamma), sprintf("%.3f", js$corr))) knitr::kable(supporting, row.names = FALSE) ``` The second-order estimate subtracts from the corrected estimate the mean of the simulated selection root; the bound subtracts a quantile of the same law and is not centred on it. A scale ratio near one means field-s and the unstudentized complement bound agree on these data. ## The two model-based identifiers DINA and GRF propose candidates from a surface fitted to the outcome and have them re-ranked on the same Cox coefficient $\hat\beta(g)$ against the same screening threshold; neither applies the consistency screen. Supplementary Section S5 of León and Anderson (2026b) sets the three identifiers side by side, with every setting used here. The selection rule and the inference are the same. The calls differ from forest search's only in the identifier's own arguments. ```{r fit-dina} #| message: false #| results: hide t0 <- proc.time()[["elapsed"]] fit_dina <- 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 = "dina", dina_select_statistic = "effect", # DINA's own hr.threshold = 1.00, stop_threshold = NULL, sg_focus = "effMaxSG", selection_rule = "neighborhood", effect_neighborhood = 0.20, mr_inference = TRUE, mr_inference_args = list(ci_method = "field", draws = 5000L, include_complement = TRUE, field_scale_complement = "selected", field_decompose = TRUE), quiet = TRUE, parallel_args = list(plan = "sequential")) fit_sec_dina <- proc.time()[["elapsed"]] - t0 ``` ```{r fit-grf} #| message: false #| results: hide t0 <- proc.time()[["elapsed"]] fit_grf <- 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 = "grf", grf_selection = "frontier", # GRF's own grf_select_statistic = "effect", grf_depth = 2L, dmin.grf = 0.0, # GRF's own hr.threshold = 1.00, stop_threshold = NULL, sg_focus = "effMaxSG", selection_rule = "neighborhood", effect_neighborhood = 0.20, mr_inference = TRUE, mr_inference_args = list(ci_method = "field", draws = 5000L, include_complement = TRUE, field_scale_complement = "selected", field_decompose = TRUE), quiet = TRUE, parallel_args = list(plan = "sequential")) fit_sec_grf <- proc.time()[["elapsed"]] - t0 ``` ```{r dina-grf-selected} g_dina <- fit_dina$mr_inference g_grf <- fit_grf$mr_inference cat(sprintf("DINA: %s (fit %.1f s; re-selection family %d)\n", paste(fit_dina$sg.harm, collapse = " & "), fit_sec_dina, g_dina$n_family)) cat(sprintf("GRF: %s (fit %.1f s; re-selection family %d)\n", paste(fit_grf$sg.harm, collapse = " & "), fit_sec_grf, g_grf$n_family)) ``` GRF returns the `r g_grf$n_selected` patients with no measured estrogen receptor, a region that contains forest search's; DINA returns the `r g_dina$n_selected` who are grade 3 with progesterone receptor at most 20 fmol/l, a different classification. For the two model-based identifiers the candidate family is generated from a surface fitted to the outcome, so their bounds are conditional on the family proposed (Supplementary Section S1.9 of León and Anderson 2026b); forest search's family is enumerated from the covariates alone. ```{r agreement} member <- function(fit) { fit$df.est$treat.recommend[match(df$id, fit$df.est$id)] == 0 } m <- data.frame(FS = member(fit_fs), DINA = member(fit_dina), GRF = member(fit_grf)) partition <- aggregate(cbind(N = 1, events = df$status), by = m, FUN = sum) knitr::kable(partition[order(-partition$N), ], row.names = FALSE, caption = "The three-way partition of the 686 patients by the three selections.") ``` ## Table 1 ```{r table1} table1 <- data.frame(row = table1_rows, FS = table1_column(fit_fs, df), DINA = table1_column(fit_dina, df), GRF = table1_column(fit_grf, df)) caption1 <- sprintf(paste( "Selection-adjusted analysis of the three identified regions, hazard-ratio", "scale. Unadjusted: two-sided 95%% Wald interval; field and field-s: one-sided", "at nominal 0.95; Bonferroni pair: region / complement, joint 0.95.", "Whole trial: %d patients, %d events, hazard ratio %.3f (%.3f, %.3f)."), nrow(df), sum(df$status), whole_hr[["est"]], whole_hr[[2]], whole_hr[[3]]) knitr::kable(table1, col.names = c("", "FS", "DINA", "GRF"), row.names = FALSE, align = "lrrr", caption = caption1) ``` This is Table 1 of León and Anderson (2026b). The calls, the seed and the draw counts are those of the manuscript's analysis, so the published values are reproduced to the digit; a different draw count would move the corrected estimates and the bounds in the third decimal and leave the regions and the unadjusted estimates unchanged. ### Reading the table Each estimate is read by its location against the trial's own effect, `r sprintf("%.3f", whole_hr[["est"]])`; the bounds are reported beside it and are not read as tests at a hazard ratio of one. A corrected estimate above 0.90 is read as suggestive of limited benefit, and where its lower bound lies below 1.0 a benefit cannot be ruled out. ```{r reading-values} #| echo: false r3 <- function(x) sprintf("%.3f", x) ev_region <- vapply(m, function(z) sum(df$status[z]), numeric(1)) ``` - **Forest search.** The unadjusted hazard ratio in its region is `r r3(g_fs$naive$est)`, opposite in direction to the trial's. The correction leaves `r r3(g_fs$debiased$est)`, above 0.90, so suggestive of limited benefit; with a field lower bound of `r r3(g_fs$field$lower_1s)`, below 1.0, a benefit cannot be ruled out. The complement, `r round(100 * g_fs$complement$n / nrow(df))`% of the trial, has a corrected estimate of `r r3(g_fs$complement$debiased$est)`, a larger benefit than the trial's, with a field-s upper bound of `r r3(g_fs$field$complement$upper_1s_s)`. - **DINA and GRF** read the same way: corrected estimates of `r r3(g_dina$debiased$est)` and `r r3(g_grf$debiased$est)` in their regions, with lower bounds below 1.0, and of `r r3(g_dina$complement$debiased$est)` and `r r3(g_grf$complement$debiased$est)` in their complements. - **The two-subgroup claim.** Forest search's Bonferroni pair is `r r3(g_fs$field$joint_s$bonf_lower_H)` for the region and `r r3(g_fs$field$joint_s$bonf_upper_Hc)` for the complement. The complement's benefit, a hazard reduction of at least `r round(100 * (1 - g_fs$field$joint_s$bonf_upper_Hc))`%, survives the joint adjustment; the region excludes only hazard ratios below `r r3(g_fs$field$joint_s$bonf_lower_H)`, which is no claim of harm. A claim that the two subgroups differ is not supported; the claim of benefit outside the region is. - **The substantive reading** is the same under all three partitions: a small region whose unfavourable unadjusted estimate the correction largely removes, and a complement, most of the trial, with a corrected benefit larger than the trial's. With `r min(ev_region)` to `r max(ev_region)` events in the regions, a claim of harm at a conventional level is not realistic at this trial's size; the bounds are reported, and how much statistical assurance to require of them is the analyst's decision in their own setting. ## Declaration at a protected level Every region above was declared at the identifiers' default settings, forest search's by the conventional screen at $p^\star = 0.90$. Section 4 of León and Anderson (2026b) asks a different question: would a screen calibrated to protect a clinically meaningful uniform benefit have declared it? The answer is read for forest search, whose family is enumerated. At a protected level $c_0$ and family-wise level $\alpha$, the calibrated cutoff $\hat\kappa_\alpha(c_0)$ is the $(1 - \alpha)$ quantile of the maximum over the family of each candidate's standardized statistic $$ T(g) \;=\; \frac{\hat\beta(g) - \log c_2}{\sigma_D(g)}, $$ with every candidate shifted to sit at $c_0$, under the law of the fit's own multiplier draws. The calibrated rule admits a candidate when $T(g) \ge \hat\kappa_\alpha(c_0)$; it holds the family-wise probability of declaring any harm subgroup at $\alpha$ whenever every candidate's true effect is $c_0$ or better. `fs_declaration_calibration()` reads it from the draws the fit above already made (`keep_declaration_field = TRUE`, with the grid passed as `declaration_c0`); nothing is redrawn and nothing else in the fit changes. ```{r declaration} c0_grid <- c(0.70, 0.75, 0.80, 0.85) cal <- fs_declaration_calibration(fit_fs, alpha = 0.10, c0 = c0_grid) cal_c2 <- fs_declaration_calibration(fit_fs, alpha = 0.10) # c0 = NULL: the claim threshold tb <- cal$c0$table fld <- g_fs$declaration_field T_H <- cal$T_hat[[fld$family_id[fld$meta$selected_index]]] T_max <- max(cal$T_hat) kappa <- c(tb$kappa_hat, cal_c2$kappa_hat) declaration <- data.frame( c0 = c(sprintf("HR %.2f%s", tb$c0, ifelse(tb$c0 == 0.75, " (primary)", "")), "HR 1.00 (claim threshold)"), kappa = sprintf("%.3f", kappa), implied = sprintf("%.3f", 2 * pnorm(kappa) - 1), fw = sprintf("%.3f", c(tb$fw_size, cal_c2$fw_size)), admits = c(lengths(cal$c0$admitted_calibrated), length(cal_c2$admitted_calibrated))) knitr::kable(declaration, row.names = FALSE, align = "lrrrr", col.names = c("protected level c0", "kappa-hat", "implied level", "FW size", "admits"), caption = paste( "Declaration calibration at alpha = 0.10. kappa-hat: the calibrated", "cutoff; implied level: 2*Phi(kappa-hat) - 1; FW size: the family-wise", "size of the p* = 0.90 screen at c0; admits: candidates the calibrated", "rule admits.")) cat(sprintf("T of the declared region: %.3f; largest T in the family: %.3f\n", T_H, T_max)) cat(sprintf("The p* = 0.90 screen's cutoff, as the package applies it: %.3f\n", cal$z_pstar)) ``` The screen's cutoff is the $p^\star = 0.90$ rule as the package applies it, after the rate is rounded to `pconsistency.digits = 2`, which is why it sits a little below the 1.645 of an unrounded 0.90. At the primary protected level, a hazard ratio of 0.75, the calibrated cutoff is `r sprintf("%.3f", tb$kappa_hat[tb$c0 == 0.75])`, an implied screening level of `r sprintf("%.3f", 2 * pnorm(tb$kappa_hat[tb$c0 == 0.75]) - 1)`; the region's statistic is `r sprintf("%.3f", T_H)` and the largest in the family `r sprintf("%.3f", T_max)`, so the calibrated screen admits `r if (lengths(cal$c0$admitted_calibrated)[tb$c0 == 0.75] == 0) "no candidate" else "a candidate"` at that level, and likewise at every protected level on the grid and at the claim threshold. The FW-size column is the diagnostic of the paper's equation (8): the probability that the $p^\star = 0.90$ screen as executed declares some harm subgroup if every candidate carried a uniform benefit of $c_0$ — a property of this family and covariate design, not the trial's realized error rate; at the claim threshold it is a boundary size. To run the screen at the implied level, set `pconsistency.threshold` to it and `pconsistency.digits` to at least 5, so the level is not rounded away; the package's `pstar_settable` and `pstar_fine` columns in `cal$c0$table` give what the digit grid can express. The supplement's Section S1.8.4 adds a realized-rate simulation on this design, which is not run here. The selection and every bound above are those of the fit as run; the calibration reads the fit and alters none of them. ## Where to go next - The [post-selection article](https://larry-leon.github.io/forestsearch/articles/survival-post-selection.html) 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 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 for this analysis and for the binary and continuous ACTG 175 analyses of León and Anderson (2026b), 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. León, L. F. and Anderson, K. M. (2026b). Interpretable subgroup identification across survival and generalized-linear outcomes. 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*. ```{r session} sessionInfo() ```