--- title: "Getting Started with hvtiRbootstrap" author: "John Ehrlinger" date: today format: html: toc: true vignette: > %\VignetteIndexEntry{Getting Started with hvtiRbootstrap} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ## What is this package for? `hvtiRbootstrap` is the R port of the CORR macro library's bootstrap macros, written for the analysts who already run them. The job is the one `%bootreg` does: fit the same selection procedure on each of many bootstrap replicates, record which terms survived each time, and report how often each one appeared. A variable chosen in 95% of replicates is a different kind of finding from one chosen in 55%, even though a single stepwise run on the full data would have reported both the same way. If you have run `%bootreg` and `%SUMBOOT`, the core here is the same three steps: | SAS macro | R function | does | |---|---|---| | `%bootreg` | `boot_select()` | resample, fit, record the terms each model kept | | `%SUMBOOT` | `boot_summary()` | selection frequency and coefficient spread per term | | `%cluster` | `boot_clusters()` | how often at least one of a correlated group was kept | `PROC=` becomes the `fitter` argument: `fit_logistic()`, `fit_linear()` or `fit_cox()`. This vignette walks one small screen from start to finish, then shows the reporting layer and the interval branch that sit alongside it. Everything here runs on simulated data: this vignette, the package's examples and its tests contain no cohort data and no PHI. ```{r library} library(hvtiRbootstrap) ``` ## A small screen We simulate 200 patients with a binary outcome. Two things drive the risk: age and body surface area (`bsa`). Weight (`wt`) is built from `bsa`, so it carries the same information as a close correlate, the way the two do in real data. Ejection fraction (`ef`) and `noise` carry no signal at all. ```{r simulate} set.seed(1) n <- 200 age <- rnorm(n, mean = 65, sd = 10) bsa <- rnorm(n, mean = 1.9, sd = 0.2) wt <- 40 * bsa + rnorm(n, sd = 4) ef <- rnorm(n, mean = 55, sd = 8) noise <- rnorm(n) risk <- plogis(-1 + 0.06 * (age - 65) + 2 * (bsa - 1.9)) sim <- data.frame( dead = rbinom(n, 1, risk), age, bsa, wt, ef, noise ) round(cor(sim$bsa, sim$wt), 2) ``` ## Resample and fit `boot_select()` is `%bootreg`. Give it the data, a formula offering the candidate terms, and a fitter. The arguments you know carry over: `n_rep` is `RESAMPL=`, `sle` and `sls` are `SLE=` and `SLS=` (defaults 0.1 and 0.05), `fraction` is `FRACTION=` and `seed` is `SEED=`. The default `fraction = 1` draws as many rows as the data has, which is what the macro always does. ```{r select} fit <- boot_select( sim, dead ~ age + bsa + wt + ef + noise, fitter = fit_logistic, n_rep = 50, seed = 101 ) fit ``` We ask for only 50 replicates so the vignette builds quickly. A real screen uses the default of 1000, and `n_rep` counts *valid* models, as `RESAMPL=` does: a replicate whose fit fails is drawn again and not counted. The result keeps one row per replicate and one column per candidate term. A term the model did not select is `NA` in that row: ```{r coefficients} head(fit$coefficients) ``` That missingness is the design, not a gap to clean up. Counting the non-missing values down a column gives the number of replicates that chose the term, which is how `%SUMBOOT` counts too. Fill those `NA`s with zero and every selection frequency in the package becomes 100%. ## How often did each term survive? `boot_summary()` is `%SUMBOOT`, held to exact parity with the macro: the same replicate table in gives the same `n`, `pct`, mean, standard deviation, minimum and maximum out. ```{r summary} s <- boot_summary(fit) s ``` ```{r pct-helper, include = FALSE} pct <- function(v) s$pct[s$variable == v] ``` `n` is the count of replicates that selected the term and `pct` is that count over `n_rep`. The two real effects hold up. Age survives in `r pct("age")`% of replicates and `bsa` in `r pct("bsa")`%. `ef` and `noise` carry nothing by construction, yet they still turn up in `r pct("ef")`% and `r pct("noise")`%. That is stepwise keeping a useless term by chance, and a fair reminder that a non-zero `pct` is not evidence of anything on its own. The coefficient columns describe the replicates that *selected* the term, not all of them. A term chosen a handful of times has a mean and spread computed over that handful. ## What about correlated terms? Look at `wt`. It is selected in only `r pct("wt")`% of replicates, which reads as unimportant. But `wt` and `bsa` measure the same thing, and once one of them is in the model the other has little left to add. The replicates split their votes between the pair, so each one's individual frequency understates how often "body size" made it into the model. `boot_clusters()` is `%cluster`. It asks how often *at least one* member of a group was selected: ```{r clusters} cl <- boot_clusters(fit, list(size = c("bsa", "wt"))) cl ``` `n_any` is not the sum of the members' counts. A replicate that selected both `bsa` and `wt` counts once, which is why the cluster reads `r cl$pct_any`% rather than `r pct("bsa") + pct("wt")`%. ## Other model families The fitter is the only thing that changes between a logistic, linear and Cox screen. Each pins its own entry and removal tests to the `PROC=` it replaces: `fit_linear()` uses partial F like `PROC REG`, `fit_logistic()` enters on the score test and removes on Wald like `PROC LOGISTIC`, and `fit_cox()` enters on the likelihood ratio and removes on Wald. That last one is a registered divergence from `PROC PHREG`, which enters on the score test; R has no score test for a Cox model to call. `fit_cox()` needs the `survival` package, so this chunk runs only when it is installed: ```{r cox, eval = requireNamespace("survival", quietly = TRUE)} set.seed(2) sim$time <- rexp(n, rate = 0.1 * exp(0.05 * (sim$age - 65))) sim$event <- rbinom(n, 1, 0.7) cox_fit <- boot_select( sim, survival::Surv(time, event) ~ age + ef + noise, fitter = fit_cox, n_rep = 20, seed = 202 ) boot_summary(cox_fit) ``` A Cox model has no intercept, so there is no `(Intercept)` row here. ## From a screen to a report The reporting layer turns a finished screen into tables. It reads a *long* "bag" rather than the wide matrix `boot_select()` returns, and `boot_bag()` is the conversion. It takes the facts the screen cannot know about itself: which terms are the base model, how many candidates the runner offered, and which dataset was screened. ```{r bag} bag <- boot_bag( fit, base_params = "(Intercept)", requested = 5, manifest = list(sha256 = "synthetic-example") ) ``` `boot_frequencies()` is the per-term table a report shows. It drops the base model, which is in every replicate by construction, and adds the Monte-Carlo error each frequency carries, roughly `sqrt(p * (1 - p) / n)`. Give it a retention cutoff and it marks which terms clear it, and which sit close enough that a rerun with different seeds could move them across. ```{r frequencies} boot_frequencies(bag, threshold = 50) ``` With 50 replicates the errors are several percentage points wide. That is one more reason a real screen runs 1000. `boot_health()` asks a different question: did the screen run at all? A screen that selected nothing, or one whose base parameter never varied across replicates, can still produce a tidy table of numbers that mean nothing. This check catches both. ```{r health} boot_health(bag) ``` ## Banding an estimate instead Everything above is the *selection* branch. The replicates are a vote, and `NA` is how a vote is cast. The package does a second job as well. `boot_predict_ci()` is the port of `%BNMNR` and `%BNPREV`. It resamples to put a band around an estimate, so nothing is selected and there is no `NA` semantics: the replicates are a distribution. You supply a `statistic`, a function that takes one resampled data frame and returns a named numeric vector. Here it is the predicted risk at three ages: ```{r intervals} risk_at_age <- function(d, ...) { m <- glm(dead ~ age, family = binomial, data = d) p <- predict(m, newdata = data.frame(age = c(55, 65, 75)), type = "response") stats::setNames(p, c("age55", "age65", "age75")) } ci <- boot_predict_ci(sim, risk_at_age, n_rep = 50, seed = 7) summary(ci) ``` No argument sets a confidence level, here or anywhere in the package. The macros hardcode the 2.5, 16, 50, 84 and 97.5 percentiles, so both the 95% band (`cll_p95`, `clu_p95`) and the 68% band (`cll_p68`, `clu_p68`) come back in columns named for their coverage. The percentiles use `quantile(type = 4)`, which is SAS's `PCTLDEF=1`. The bands are pointwise: each quantity is summarised on its own, so a curve drawn through `cll_p95` is not a 95% region for the curve. ## What matches SAS, and what does not `boot_summary()` and `boot_clusters()` are held to exact parity with `%SUMBOOT` and `%cluster`, and the interval arithmetic in `boot_predict_ci()` matches `PROC STDIZE PCTLDEF=1`. Resampling and model fitting are not parity-tested and cannot be, because R and SAS draw different samples. So compare an R screen with a SAS screen of the same data on its conclusions, not to the decimal. Where the port departs from a macro on purpose, the README's divergence register lists the departure and the function's help page marks it as a **Divergence**. Most are set up so the default still behaves as SAS does. Two are not: `boot_predict_ci()` defaults to 1000 replicates where `%BNMNR` defaults to 100 (pass `n_rep = 100` to match), and `fit_cox()` enters on the likelihood ratio because there is no score test to call. ## Where to go next - The function help pages, `?boot_select` first. Each names the macro it ports and the SAS parameter each argument replaces. - The [reference index](https://ehrlinger.github.io/hvtiRbootstrap/reference/index.html), grouped by job: resample, summarise, report, pool and band. - `boot_pool_chunks()`, `boot_chunk_files()` and `boot_shortfall()`, when a screen is too long to run in one go and has to be split into restartable chunks. - `boot_provenance()`, `boot_seeds()`, `boot_dropped()` and `boot_concepts()` for the rest of a report: where the screen came from, the candidates dropped before screening, and competing forms of one variable counted together. - The [README](https://github.com/ehrlinger/hvtiRbootstrap#readme) for the full SAS-to-R argument map and the divergence register.