forestsearch

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

Releases

2

Needs compilation

No

Needs R

>= 4.1.0

Config/roxygen2/version

8.1.0

Config/testthat/edition

3

SystemRequirements

Quarto command line tools (https://github.com/quarto-dev/quarto-cli)

• Clinical trial design

Since March 2026

• data.table
• doFuture
• dplyr
• foreach
• future
• future.apply
• future.callr
• ggplot2
• glmnet
• grf
• gt
• Matrix
• mvtnorm
• parallel
• patchwork
• policytree
• progressr
• randomForest
• rlang
• stats
• survival
• utils
• weightedsurv

Installing it adds

89 packages, 82.4 MB of source

Of those, need compilation

43

• cubature
• DiagrammeR
• doRNG
• forestploter
• htmltools
• katex
• knitr
• latentcor
• MASS
• quarto
• rmarkdown
• sandwich
• svglite
• testthat
• tidyr

Lines by kind (118,185)

Lines by language

Languages are counted from file extensions.

Files

618

R functions

705

Object systems

S3: 60 methods

Deprecated functions

3

Median

5

Mean

10.63

Highest

303

Median lines per function

30

R functions

705

Calls

761

Unconnected

130

Groups

28

Lines counted

33,573

Lines the tests ran

12,793

Tests passed

No

Measured with

covr 3.6.5 on R 4.6.1

Has tests

Yes

Test frameworks

testthat, tinytest

Test cases

485

Lines of tests

12,265

Test lines per line of code

0.16

• Exploratory Subgroup Identification with ForestSearch
• Post-selection inference for an identified survival subgroup

README

Yes, 814 words

Release notes file

Yes

Help pages

427

Lines of documentation

27,531

Commits

2,713

Open pull requests

0

Tagged releases

1

Latest tagged release

v0.4.0, 10 Oct 2026

In development

Default branch at 0.4.0.9000, ahead of CRAN 0.4.0 (checked 2026-10-25)

Issues closed in

–

Pull requests closed in

0 days (median)

Pull requests merged

67%

Created

7 Oct 2025

Licence file reads as

–

Method, repository practices

Across larry-leon

Changes since the last CRAN release (0.1.0)

Breaking changes

• showten_subgroups is removed from forestsearch() and subgroup.consistency(). Neither function has ..., so a call that passes it now fails with "unused argument". Use show_candidate_summary instead.
• Changed defaults in forestsearch() (old -> new):
• use_lasso: TRUE -> FALSE
• use_grf: TRUE -> FALSE
• dmin.grf: 12 -> 0
• frac.tau: 0.6 -> 0.8
• d0.min: 12 -> 10
• d1.min: 12 -> 10
• max_subgroups_search: 10 -> Inf
• vi.grf.min: -0.2 -> NULL
• stop_threshold: 0.95 -> pconsistency.threshold
• parallel_args: list(plan = "callr", workers = 6, show_message = TRUE) -> list(plan = "multisession", workers = .default_parallel_workers(), show_message = TRUE)
• Changed default in subgroup.consistency(): stop_Kgroups: 10 -> 200.
• Changed default in forestsearch_Kfold(): Kfolds: nrow(fs.est$df.est) -> .fs_cv_nobs(fs.est).
• Changed default in cox_ahr_cde_analysis(): output_dir: tempdir() -> "plots" (files are written only when save_plots = TRUE).
• save_leff_calibration() now requires output_dir. Called without it, the function stops instead of writing to "_data" in the working directory; pass output_dir = "_data" to keep the former behaviour.

New features

• GLM outcomes: binary (odds ratio), continuous (mean difference) and count/rate (incidence rate ratio) endpoints, selected with the outcome_type and effect_measure arguments of forestsearch().
• DINA and GRF as subgroup identifiers alongside the consistency search, selected with forestsearch(subgroup_method = c("consistency", "dina", "grf")); see dina_subgroup() and, for GLM outcomes, grf.subg.harm.glm().
• Post-selection inference: multiplier resampling (fs_mr_inference()), with infinitesimal-jackknife (IJ) standard errors and field bounds.
• Operating-characteristics functions: fs_oc_family_enumerate(), fs_oc_predict(), fs_oc_grid() and fs_oc_invert().
• Threshold rules on the forest-search consistency path (survival HR and binary OR only): when only c1 is given, c2 = 0.8 * c1 on the ratio scale; c2 > c1 is an error.
• Calibrated declaration threshold: fs_declaration_calibration().
• forestsearch() gains pconsistency.digits (default 2): the consistency proportion is rounded to this many digits before it is compared with pconsistency.threshold.
• grf forest training is capped at 2 threads when _R_CHECK_LIMIT_CORES_ is set (as under R CMD check --as-cran); otherwise grf's own default applies unchanged.
• Printed post-selection output and documentation simplified; methodology reference León and Anderson (2026) <doi:10.48550/arXiv.2609.38361>.

Development changes after 0.3.5 (not released on CRAN)

• MR's admission set now uses the screen's rounded consistency threshold (forest search only). The search admits on round(Pcons, pconsistency.digits) >= p*, but MR's consistency floor was built from the exact cutoff qnorm((1 + p*) / 2), so MR re-selected under a different admission rule from the one that produced the selected subgroup, in the under-adjusting direction. The floor is now built from the effective threshold .fs_pcons_eff(p*, digits) on both consistency paths. fs_mr_inference() gains pconsistency.digits (default NULL, the subgroup.consistency() default); forestsearch() passes its own value. Corrected estimates, field bounds and the Bonferroni pair can move for forest search. On the GBSG application (seed 8316951) the corrected HR moves from 1.2806 to 1.2762, the field lower bound from 0.6365 to 0.6336, and the Bonferroni pair from 0.5099 / 0.8397 to 0.5087 / 0.8396. GRF and DINA carry no consistency screen and are unaffected (identical output).
• `fs_declaration_calibration()`: FW-hat and the implied p-star are now computed at the screen's rounded admission rule. The screen admits on round(Pcons, pconsistency.digits) >= p*, which is Pcons >= 0.895 for p* = 0.90 at 2 digits, but fw_size was taken at the exact qnorm((1 + p*) / 2). It is now taken at the effective threshold qnorm((1 + pcons_eff) / 2) (returned as pcons_eff and z_pstar), with digits read from the fit's args_call_all (default 2, source recorded), so `fw_size` values rise relative to previous releases and now depend on pconsistency.digits; admitted_pstar uses the same threshold. On the GBSG application (c0 = 0.75, alpha = 0.10) fw_size(c0) moves from 0.6508 to 0.6672. `kappa_hat` is unchanged. The c0 table's pstar_implied column is replaced by what to set to run kappa_hat as a p-star screen: pstar_settable (smallest settable p* at the fit's digits whose rounded screen is at or above kappa_hat; pstar_achievable = FALSE when none is), pcons_eff_settable, z_eff_settable, z_gap, and digits_fine / pstar_fine / z_gap_fine (smallest digits with a gap under 0.01 in z). Under consistency_method = "split" the screen is not a threshold on T, so fw_size, admitted_pstar and the settable columns are NA / NULL.
• The fs_dgm_feasibility() example now runs (it called generate_glm_dgm() without its required data arguments); fs_plot_bias_coverage() captions use \u escapes, so the R source is ASCII -- the rendered text is unchanged.
• `forestsearch()` gains `pconsistency.digits` (default 2). It was previously reachable only on the internal consistency functions (subgroup.consistency(), evaluate_subgroup_consistency(), evaluate_consistency_twostage()); every forestsearch() run was pinned to 2 digits. The consistency proportion is rounded to this many digits *before* it is compared with pconsistency.threshold, so the setting governs admission, not only display; the roxygen for all four functions now says so. Default behaviour is unchanged: at 2 the GBSG application's consistency table is identical to the previous package's. fs_family_report() lists the argument in its consistency-screen row.
• New `fs_dgm_feasibility()`: a design-time check that the planted region can be declared at all. Given a glm_dgm -- the interface fs_oc_predict() and fs_oc_grid() accept -- it draws n_rep replicates per sample size through the DGM's *own* generator (simulate_from_glm_dgm(), the function the campaign templates call, never a re-implementation) and reports, per n: the share of replicates on which the planted region is undeclarable (size <= n.min, the strict test subgroup.search() itself applies), the share under the per-arm events floor, the share on which the estimand does not exist under the estimability boundary's per-estimand condition, and size and per-arm cell summaries. $feasible is TRUE when every undeclarable share is at or below tolerance (default 0.05). The print method is loud and unconditional. On the ACTG175 OR 0.75 design this returns FALSE at n = 500 and n = 750: the planted region is at or below n.min = 60 in 90% of replicates at n = 500, so what the search declares there is a larger overlapping region, not the planted one. It imposes nothing -- the floors are read and reported, never applied, and forestsearch() never consults it. The function does not change the RNG kind and restores the caller's RNG stream on exit; the caller's matching obligation, stated in the roxygen, is to build and calibrate the DGM *before* any replicate switches the kind, since calibrating afterwards yields a different super-population. Survival DGMs are not accepted, and the error and ?fs_dgm_feasibility say exactly what supporting them would need rather than guessing: setup_gbsg_dgm() objects carry a different generator (simulate_from_dgm()) whose analysis_time and cens_adjust arguments have no GLM counterpart and change the per-arm event counts.
• The estimator boundary now returns `NA` with a reason where the estimand does not exist, instead of a finite divergent number. glm() and coxph() do not fail on a slice with an empty cell or a zero-event arm: they return a large finite coefficient with converged = TRUE, which was then admitted as an estimate (a zero-cell OR returned 19.97 with SE 4809). The existence condition is now checked before the fit, per estimand: OR requires all four cells (control/treated x events/non-events) >= 1; RR, IRR and HR require at least one event in each arm. A candidate failing it returns estimate = NA, se = NA, converged = FALSE and a reason naming the empty cell, and takes the existing fit-failure status (status 5) -- an NA effect cannot rank, so no selection logic changes. For the ratio estimands only, a fit reporting converged = FALSE is likewise non-estimable ("non-convergent fit"). Unaffected, byte-identical: RD, IRD and MD, which exist on any slice -- in particular RD's tier-3 raw-proportions fallback returns converged = FALSE by design and is never caught. This is not an admission floor: n.min, d0.min/d1.min and every DINA/GRF floor are unchanged, nothing is forwarded to DINA or GRF, and no minimum count is imposed anywhere. Non-estimable candidates are counted, not dropped silently: filter_counts$n_nonestimable and filter_counts$nonestimable_reasons carry the tally and its reasons, subgroup.search(details = TRUE) prints them, and fs_family_report() gains an "estimability boundary" stage row. Committed results are unaffected -- no declared subgroup in the committed binary record violates the four-cell condition (3,313 replicates, a3374e6f), and no committed survival bundle stores a non-finite estimate or SE (930 bundles, 568,757 declared rows).
• `dina_frontier()`'s display caps now default to `Inf` and say when they trim. max_per_covariate and max_subgroups defaulted to 3L and 10L, which silently dropped covariates whole and truncated the round-robin rank-2 pass -- observed in the applied vignettes, where the printed frontier was not the frontier. Both now default to Inf, aligning with max_subgroups_search's own Inf doctrine, and a user-supplied finite value that actually removes rows raises a warning naming the kept and available counts (condition class dina_frontier_cap_trim). Affected: callers who relied on the finite defaults -- they now receive the full table; every committed in-repo caller passes both caps explicitly, so no committed rendered table changes. Unaffected: the candidate pool anywhere. The frontier is a display object; dina_subgroup() never consults it, and forestsearch(use_dina = TRUE) still supplies the finite screening caps itself, so the screening pool is byte-identical. Under the default selected_only = TRUE that screening frontier is discarded, so the trim warning is muffled there rather than naming a table nobody sees.
• `dina_args` frontier keys now warn when they are ignored. The seven dina_frontier() keys -- scope, m_diff, n_min, direction, max_per_covariate, max_subgroups, digits -- are inert under subgroup_method = "dina", which selects through dina_subgroup() and never consults a frontier table. Supplying any of them there now raises a single warning() per fit naming every offending key at once, rather than accepting them silently. Affected: `subgroup_method = "dina"` callers who pass frontier keys -- the selection is unchanged, only the messaging. Unaffected: the use_dina screening path, where the keys still act (digits under selected_only = TRUE, all seven under selected_only = FALSE), and every other subgroup_method.
• DINA now refuses the identity-scale estimands `"RD"` and `"IRD"`. DINA scores candidates on the family link scale, so its admission floor was derived as log(hr.threshold) for the ratio families and as the identity threshold for gaussian. For an explicit effect_measure = "RD" that derivation was meaningless: a 7-percentage-point threshold became a log-odds floor of log(0.07) = -2.66, which admitted nearly every candidate, silently. There is no identity-scale number DINA could correctly apply on a link scale, so the derivation now stops instead, naming the ratio and non-DINA alternatives. The guard sits at the m_diff derivation itself, so both routes to it are covered: subgroup_method = "dina" (an error, raised before the DINA model is fit) and use_dina screening with dina_args$selected_only = TRUE (where it surfaces through that path's existing "DINA analysis failed" warning and the run continues with no DINA cuts). Affected: only calls that explicitly request `"RD"` or `"IRD"` together with DINA -- nothing committed does, and the defaults per outcome type ("OR", "MD", "IRR") are all outside the guard. Unaffected: family = "gaussian" with "MD", which is a correct identity-scale floor and cannot reach the refusal, and every ratio estimand on every family.
• The two effect thresholds are now a pair on the forest-search consistency path. For subgroup_method = "consistency" with survival (HR) or binary effect_measure = "OR" -- the estimands whose c1 and c2 sit on one comparable ratio scale -- forestsearch() now enforces c2 <= c1 and derives a c2 the caller did not supply. Affected: callers who set `c1` without `c2` (effect.threshold or hr.threshold alone), which used to leave c2 at the package default 1.0 whatever c1 was and now resolves c2 = 0.80 * c1 on the ratio scale -- 1.25 -> 1.00, 1.00 -> 0.80, 0.90 -> 0.72 -- announced with a message and carried into every bootstrap replicate and cross-validation fold. The default pair (1.25, 1.0) is the rule's fixed point, so a call that set neither threshold resolves exactly what it did before, and a call that passed hr.threshold = 1.25 explicitly is likewise unmoved and only gains the message. Also affected: callers with `c2 > c1`, previously accepted and degenerate -- the consistency-stage entry condition re-screens the admitted family on c2, so candidates the screen let through were discarded at the stage boundary and the run could report no subgroup for a reason nothing in the output named. That configuration is now an error, raised before any model is fit, naming both values and the spellings the caller used. Supplying both spellings of one threshold at disagreeing values is an error too. An explicitly supplied c2 is never overridden, and no selection logic changed. Unaffected: subgroup_method = "dina" and "grf", which return before the consistency stage and never consult c2, and the identity-scale estimands "RD", "IRD" and "MD" together with "IRR", all of which resolve exactly as before in every branch.
• Bootstrap replicates and cross-validation folds now resolve the thresholds the original fit resolved. forestsearch() resolves the screening (c1) and consistency (c2) thresholds through branches guarded by whether the caller supplied them, and for the legacy spellings (hr.threshold, hr.consistency) that was detected with missing(). Because forestsearch_bootstrap_dofuture(), forestsearch_Kfold() and forestsearch_tenfold() replay the captured argument list with every formal supplied, missing() was FALSE inside a replicate, and a replicate could resolve thresholds its parent never used. Affected: identity-scale analyses (`effect_measure` `"RD"`, `"IRD"` or `"MD"`) that left `c1` or `c2` at its default, or supplied only one of the two, and then ran the bootstrap or cross-validation -- those replicates used c2 = 1.0 where the fit used 0.0, and on "MD" also c1 = 1.25 where the fit used 0.0. A per-split effect of 1.0 is a threshold no candidate can meet, so the affected folds and replicates identified nothing. Ratio measures ("OR", "RR", "IRR") and survival were never affected, and the fix leaves every original fit's resolution byte-identical -- only what a replay sees changes.

The notes are cut short here.