--- title: "Random Forest Survival Analysis with ggRandomForests" author: "John Ehrlinger" date: today format: html: fig-format: png fig-dpi: 96 toc: true toc-depth: 3 html-math-method: mathjax code-fold: false bibliography: ggRandomForests.bib vignette: > %\VignetteIndexEntry{Random Forest Survival Analysis with ggRandomForests} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- ::: {.callout-warning} ## Work in progress This vignette is under active development. Code examples and narrative may change before the next release. ::: ```{r setup, include=FALSE} # vignette dir under R CMD check, package root under an interactive knit .fo <- c("_fig_optim.R", file.path("vignettes", "_fig_optim.R")) if (any(file.exists(.fo))) source(.fo[file.exists(.fo)][1]) knitr::opts_chunk$set( fig.align = "center", fig.width = 7, fig.height = 5, message = FALSE, warning = FALSE, comment = "#>" ) old_opts <- options(mc.cores = 1, rf.cores = 1) ``` # Introduction Random forests [@Breiman:2001] are a non-parametric ensemble method that requires no distributional assumptions on the relation between covariates and the response. Random survival forests (RSF) [@Ishwaran:2007a; @Ishwaran:2008] extend the method to right-censored, time-to-event data by growing trees with a log-rank splitting rule and aggregating Kaplan--Meier estimates within terminal nodes. The forest returns three ensemble quantities for each subject. First, the survival function $\hat{S}(t)$: the probability of surviving past time $t$, bounded between 0 and 1. Second, the cumulative hazard function (CHF): $\hat{H}(t) = -\log \hat{S}(t)$, an unbounded, monotone-increasing summary of accumulated risk. Third, mortality: the expected number of events a subject would accumulate if followed indefinitely under their covariate profile, computed by summing the CHF over the observed time grid. Mortality is not a probability; it is an unbounded relative-risk score. Read it as a single number for ranking patients from low to high risk, not as a survival chance. The **[randomForestSRC](https://www.randomforestsrc.org)** package [@Ishwaran:RFSRC:software:2026] provides a unified implementation for survival, regression, and classification forests. **ggRandomForests** extracts tidy data objects from `rfsrc` fits and renders them with **ggplot2** [@Wickham:2009], so you can see how a survival forest is built, which predictors carry the risk, and how survival shifts across each one. This vignette demonstrates a complete random survival forest workflow on the primary biliary cirrhosis (PBC) data set [@fleming:1991]: 1. **Data preparation and exploration**: cleaning, EDA, Kaplan--Meier curves 2. **Growing the forest**: fitting an RSF, checking convergence and OOB error 3. **Variable selection**: VIMP and minimal depth via `max.subtree()` 4. **Dependence plots**: variable dependence and partial dependence via `gg_variable()` and `gg_partial_rfsrc()` 5. **Variable interactions**: conditioning plots and partial dependence surfaces ```{r packages} library(ggplot2) library(dplyr) library(tidyr) library(randomForestSRC) library(survival) if (requireNamespace("ggRandomForests", quietly = TRUE)) { library(ggRandomForests) } else if (requireNamespace("pkgload", quietly = TRUE)) { pkgload::load_all(export_all = FALSE, helpers = FALSE, attach_testthat = FALSE) } else { stop("Install ggRandomForests (or pkgload for dev builds) to render this vignette.") } theme_set(theme_bw()) # Plotting constants event_marks <- c(1, 4) event_labels <- c("Censored", "Death") event_colors <- c("steelblue", "firebrick") ``` # Data: Primary Biliary Cirrhosis (PBC) The PBC study consists of 424 patients referred to Mayo Clinic between 1974 and 1984, of whom 312 were randomized into a trial of D-penicillamine (DPCA) versus placebo. The data is described in @fleming:1991 Chapter 0.2, with a proportional hazards model developed in Chapter 4.4. We use the copy bundled with **randomForestSRC**. ```{r data-load} data("pbc", package = "randomForestSRC") ``` ## Data cleaning We convert `days` to `years`, recode `treatment` as a factor, and coerce low-cardinality numeric columns (five or fewer unique values, including binary 0/1 indicators) to factors. We avoid converting binary columns to `logical` because `randomForestSRC::partial.rfsrc()` does not handle logical predictors correctly in survival forests. ```{r data-clean} pbc <- pbc |> mutate( years = days / 365.25, age = age / 365.25, treatment = factor( ifelse(treatment == 1, "DPCA", ifelse(treatment == 2, "Placebo", NA)), levels = c("DPCA", "Placebo") ) ) |> select(-days) # Low-cardinality numerics (including binary 0/1) to factor. # NOTE: do NOT convert to logical — partial.rfsrc() fails with logical # predictors in survival forests (randomForestSRC <= 3.5.1). # Exclude the response columns (status, years) from conversion. resp_cols <- c("status", "years") for (nm in setdiff(names(pbc), resp_cols)) { v <- pbc[[nm]] if (is.numeric(v) && !is.factor(v) && length(unique(v[!is.na(v)])) <= 5) { pbc[[nm]] <- factor(v) } } ``` ```{r variable-labels} # Human-readable labels for plot axes st_labs <- c( status = "Death Event", treatment = "Treatment", age = "Age (years)", sex = "Female", ascites = "Ascites", hepato = "Hepatomegaly", spiders = "Spiders", edema = "Edema (0, 0.5, 1)", bili = "Serum Bilirubin (mg/dl)", chol = "Serum Cholesterol (mg/dl)", albumin = "Albumin (gm/dl)", copper = "Urine Copper (ug/day)", alk.phos = "Alkaline Phosphatase (U/liter)", ast = "SGOT (U/ml)", trig = "Triglycerides (mg/dl)", platelet = "Platelets (per cubic ml/1000)", prothrombin = "Prothrombin Time (sec)", stage = "Histologic Stage", years = "Follow-up Time (years)" ) ``` ## Exploratory data analysis Good practice before modeling: scan categorical variables as stacked histograms over follow-up time, and continuous variables as scatter plots colored by event status. ```{r eda-categorical, fig.height=4, fig.cap="EDA for categorical variables. Bar color indicates factor level; white = missing."} cls <- sapply(pbc, class) cnt_idx <- which(cls %in% c("numeric", "integer")) fct_idx <- setdiff(seq_along(pbc), cnt_idx) fct_idx <- union(fct_idx, which(names(pbc) == "years")) dta_cat <- suppressWarnings( pbc[, fct_idx] |> pivot_longer(-years, names_to = "variable", values_to = "value", values_transform = list(value = as.character)) ) ggplot(dta_cat, aes(x = years, fill = value)) + geom_histogram(color = "black", binwidth = 1) + labs(y = "", x = st_labs["years"]) + scale_fill_brewer(palette = "RdBu", na.value = "white") + facet_wrap(~variable, scales = "free_y", nrow = 2) + theme(legend.position = "none") ``` ```{r eda-continuous, fig.height=4, fig.cap="EDA for continuous variables. Points colored by death event; rug marks indicate missing values."} cnt_idx <- union(cnt_idx, which(names(pbc) == "status")) dta_cont <- pbc[, cnt_idx] |> pivot_longer(c(-years, -status), names_to = "variable", values_to = "value") ggplot(dta_cont |> filter(!is.na(value)), aes(x = years, y = value, color = factor(status), shape = factor(status))) + geom_point(alpha = 0.4) + geom_rug(data = dta_cont |> filter(is.na(value)), color = "gray50") + labs(y = "", x = st_labs["years"], color = "Death", shape = "Death") + scale_color_manual(values = event_colors) + scale_shape_manual(values = event_marks) + facet_wrap(~variable, scales = "free_y", ncol = 4) + theme(legend.position = c(0.8, 0.2)) ``` Look for patterns of missingness (white bars, rug marks) and extreme values that fall outside the biological range. ## Kaplan--Meier survival by treatment The Kaplan--Meier estimate is the marginal survival curve with no covariate adjustment. It is the no-covariate baseline: once we fit a forest, comparing the KM curves to the forest's OOB predictions lets us judge how much predictive signal the covariates add beyond the raw event rate. We restrict to the 312 trial patients and construct KM curves with `gg_survival`. ```{r km-survival, fig.cap="KM survival estimates by treatment group with 95% confidence bands."} pbc_trial <- pbc |> filter(!is.na(treatment)) pbc_test <- pbc |> filter(is.na(treatment)) gg_km <- gg_survival(interval = "years", censor = "status", by = "treatment", data = pbc_trial, conf.int = 0.95) plot(gg_km) + labs(y = "Survival Probability", x = "Time (years)", color = "Treatment", fill = "Treatment") + theme(legend.position = c(0.2, 0.2)) + coord_cartesian(ylim = c(0, 1.01)) ``` The curves largely overlap, consistent with the original finding that DPCA offered no clear survival benefit over placebo [@fleming:1991]. ```{r km-cumhaz, fig.cap="KM cumulative hazard estimates by treatment group."} plot(gg_km, type = "cum_haz") + labs(y = "Cumulative Hazard", x = "Time (years)", color = "Treatment", fill = "Treatment") + theme(legend.position = c(0.2, 0.8)) + coord_cartesian(ylim = c(-0.02, 1.22)) ``` We can also stratify on continuous variables. Here we reproduce the bilirubin groupings from @fleming:1991 Figure 4.4.2: ```{r km-bili, fig.width=6, fig.cap="KM survival stratified by bilirubin groups."} pbc_bili <- pbc_trial |> mutate(bili_grp = cut(bili, breaks = c(0, 0.8, 1.3, 3.4, 29))) plot(gg_survival(interval = "years", censor = "status", by = "bili_grp", data = pbc_bili), error = "none") + labs(y = "Survival Probability", x = "Time (years)", color = "Bilirubin") ``` Higher bilirubin strongly predicts worse survival, an effect the random forest will rediscover without any prior specification. # Growing a Random Survival Forest Several predictors in the PBC trial data contain missing values (cholesterol, copper, triglycerides, and others). We handle these with a two-step approach: first impute using `impute.rfsrc()`, which uses the random forest proximity structure to fill in missing values, then fit the survival forest on the complete imputed data. This keeps the fitted forest object free of imputation state, which is required for `partial.rfsrc()` to work correctly. ```{r rfsrc-fit} set.seed(42) # Step 1: impute missing values via random forest proximity pbc_imputed <- impute.rfsrc(Surv(years, status) ~ ., data = pbc_trial, ntree = 50, nimpute = 2) # Step 2: grow the survival forest on the complete imputed data rfsrc_pbc <- rfsrc(Surv(years, status) ~ ., data = pbc_imputed, ntree = 100, nsplit = 10, tree.err = TRUE, importance = TRUE) ``` The forest grew `r rfsrc_pbc$ntree` trees, splitting on `r rfsrc_pbc$mtry` randomly selected candidate variables at each node, and stopping at a minimum terminal node size of `r rfsrc_pbc$nodesize`. ## OOB error convergence ```{r error-plot, fig.height=4, fig.cap="OOB prediction error vs. number of trees."} plot(gg_error(rfsrc_pbc)) ``` The error stabilizes well before the `r rfsrc_pbc$ntree` trees grown here, so the forest is large enough for reliable predictions. ## OOB predicted survival ```{r rfsrc-predicted, fig.cap="OOB predicted survival curves. Blue = censored, red = death events."} gg_rsf <- plot(gg_rfsrc(rfsrc_pbc), alpha = 0.2) + scale_color_manual(values = event_colors) + theme(legend.position = "none") + labs(y = "Survival Probability", x = "Time (years)") + coord_cartesian(ylim = c(-0.01, 1.01)) gg_rsf ``` Each curve is one patient's OOB ensemble survival function $\hat{S}(t)$, extended to the last follow-up time. The forest never saw that patient when building its prediction, so these are genuine out-of-sample estimates. Red (death event) curves generally fall faster and reach lower survival probabilities, so the forest does separate risk groups. Comparing this spread to the marginal KM curves from the previous section shows how much the covariates tighten risk stratification. ```{r rfsrc-by-treatment, fig.cap="Median predicted survival by treatment group with 95% confidence bands."} plot(gg_rfsrc(rfsrc_pbc, by = "treatment")) + theme(legend.position = c(0.2, 0.2)) + labs(y = "Survival Probability", x = "Time (years)") + coord_cartesian(ylim = c(-0.01, 1.01)) ``` ## Test set predictions The non-trial patients (`pbc_test`) have substantial missing data. `predict.rfsrc()` handles these transparently at prediction time via `na.action = "na.impute"`. This is distinct from training-time imputation and does not affect the fitted forest object. ```{r predict-test, fig.cap="Predicted survival for non-trial patients (test set)."} rfsrc_pbc_test <- predict(rfsrc_pbc, newdata = pbc_test, na.action = "na.impute", importance = TRUE) plot(gg_rfsrc(rfsrc_pbc_test), alpha = 0.2) + scale_color_manual(values = event_colors) + theme(legend.position = "none") + labs(y = "Survival Probability", x = "Time (years)") + coord_cartesian(ylim = c(-0.01, 1.01)) ``` Because the training curves use OOB estimates, both plots represent out-of-sample predictions and are directly comparable. # Variable Selection Random forest uses all available predictors. To understand which matter most, we examine variable importance (VIMP) and minimal depth. ## Variable importance (VIMP) VIMP measures the increase in prediction error when a variable's values are randomly permuted across the out-of-bag sample. For survival forests the prediction error is the Harrell C-statistic complement (1 − C), so permuting an important variable hurts C and yields a large positive VIMP. A negative VIMP means the variable is adding noise: the forest would predict better without it. See the [randomForestSRC](https://www.randomforestsrc.org) documentation for details on the VIMP calculation for censored outcomes [@Ishwaran:2007a]. ```{r vimp-plot, fig.width=6, fig.cap="Variable importance ranking. Blue = positive VIMP, red = negative."} plot(gg_vimp(rfsrc_pbc), labels = st_labs) + theme(legend.position = c(0.8, 0.2)) + labs(fill = "VIMP > 0") ``` Bilirubin ranks highest, followed by edema, ascites, albumin, copper, age, and prothrombin time, which closely matches the variables selected in the @fleming:1991 proportional hazards model. ## Minimal depth Minimal depth [@Ishwaran:2010] ranks variables by how close to the root node they first split, on average. Variables that partition large portions of the population early are considered most important. ```{r varsel} md_pbc <- max.subtree(rfsrc_pbc) ``` The `max.subtree()` function computes minimal depth for each variable. The threshold is `r round(md_pbc$threshold, 2)`, selecting `r length(md_pbc$topvars)` variables: `r paste(md_pbc$topvars, collapse = ", ")`. Both selection methods agree on the key continuous predictors: `bili`, `albumin`, `copper`, `prothrombin`, and `age`. Both also rank the categorical `edema` and `ascites` highly; we carry `edema` (selected by the @fleming:1991 model) forward for the remainder of the analysis. ```{r select-vars} xvar <- c("bili", "albumin", "copper", "prothrombin", "age") xvar_cat <- "edema" xvar_all <- c(xvar, xvar_cat) ``` # Variable Dependence ## Variable dependence plots Variable dependence shows each patient's predicted survival probability $\hat{S}(t_0)$ at a fixed time horizon $t_0$, plotted against a predictor of interest. Because $\hat{S}(t_0)$ is a probability, the y-axis is always bounded to [0, 1]. Points are colored by event status; a loess smooth traces the overall trend while the individual points show how much variance remains after the forest accounts for all other covariates. ```{r vardep-bili, fig.height=4, fig.cap="Variable dependence of survival at 1 and 3 years on bilirubin."} gg_v <- gg_variable(rfsrc_pbc, time = c(1, 3), time_labels = c("1 Year", "3 Years")) plot(gg_v, xvar = "bili", alpha = 0.4) + labs(y = "Survival", x = st_labs["bili"]) + theme(legend.position = "none") + scale_color_manual(values = event_colors, labels = event_labels) + scale_shape_manual(values = event_marks, labels = event_labels) + coord_cartesian(ylim = c(-0.01, 1.01)) ``` Survival drops sharply with increasing bilirubin, and the 3-year curve drops further than the 1-year curve, suggesting a non-proportional hazards effect. ```{r vardep-panel, fig.height=5, fig.cap="Variable dependence at 1 and 3 years for continuous predictors."} plot(gg_v, xvar = xvar[-1], panel = TRUE, alpha = 0.4) + labs(y = "Survival") + theme(legend.position = "none") + scale_color_manual(values = event_colors, labels = event_labels) + scale_shape_manual(values = event_marks, labels = event_labels) + coord_cartesian(ylim = c(-0.05, 1.05)) ``` The plots confirm survival increases with albumin and decreases with copper, prothrombin time, and age. The divergence between time curves for `copper` further supports a non-proportional hazards mechanism. ```{r vardep-categorical, fig.height=4, fig.cap="Variable dependence on edema (categorical)."} plot(gg_v, xvar = xvar_cat, alpha = 0.4) + labs(y = "Survival") + theme(legend.position = "none") + scale_color_manual(values = event_colors, labels = event_labels) + scale_shape_manual(values = event_marks, labels = event_labels) + coord_cartesian(ylim = c(-0.01, 1.02)) ``` Patients with edema = 1 (edema despite diuretics) have markedly lower predicted survival. ## Partial dependence Partial dependence integrates out the effects of other covariates, giving a risk-adjusted view of how each predictor influences the response [@Friedman:2000]. We use `gg_partial_rfsrc()` which calls `randomForestSRC::partial.rfsrc()` directly and returns a `gg_partial_rfsrc` object. For survival forests, the result includes a `time` column so `plot(pd)` produces one curve per predictor value over time, faceted by variable name. ```{r partial-dep, error=TRUE, fig.height=5, fig.cap="Partial dependence of survival at approximately 1 and 3 years on continuous predictors."} # partial.rfsrc() requires times that match the model's time.interest grid; # gg_partial_rfsrc() snaps the requested values to the nearest observed times. ti <- rfsrc_pbc$time.interest t1yr <- ti[which.min(abs(ti - 1))] t3yr <- ti[which.min(abs(ti - 3))] pd <- gg_partial_rfsrc(rfsrc_pbc, xvar.names = xvar, partial.time = c(t1yr, t3yr)) # Quick S3 plot — survival forests produce time-series curves per predictor value plot(pd) ``` For a publication-ready layout with custom color scale, access `pd$continuous` directly: ```{r partial-dep-custom, error=TRUE, fig.height=5, fig.cap="Partial dependence (custom styling)."} ggplot(pd$continuous, aes(x = x, y = yhat, color = factor(round(time, 2)), group = factor(time))) + geom_line(linewidth = 1) + facet_wrap(~name, scales = "free_x") + labs(y = "Predicted Survival", x = "", color = "Year") + scale_color_manual(values = setNames(c("steelblue", "firebrick"), as.character(round(c(t1yr, t3yr), 2)))) + theme_bw() ``` The partial dependence curves at approximately 1 and 3 years confirm the variable dependence findings and support the log-transforms used in the @fleming:1991 model for `bili`, `albumin`, and `prothrombin`. The divergence between time curves for `bili` and `copper` supports non-proportional hazards. ## How far can these curves be trusted? Partial dependence curves are smooth enough to read like a fitted model, so it helps to remember what the forest had to work with. The PBC trial has `r nrow(pbc_trial)` patients and `r sum(pbc_trial$status)` deaths, and in a survival forest the events carry the information. Censored patients tell the forest who was still at risk, but every split and every terminal node estimate is built from the deaths. Three things follow when events are sparse. First, the tails are thin. Patients in the top tenth of bilirubin account for `r sum(pbc_trial$status[pbc_trial$bili >= quantile(pbc_trial$bili, 0.9)])` deaths, and the bottom tenth for `r sum(pbc_trial$status[pbc_trial$bili <= quantile(pbc_trial$bili, 0.1)])`. The low end of the bilirubin curve rests on a handful of events, and a cohort with fewer deaths overall would be thin at both ends. Second, trees cannot extrapolate. The terminal node that covers the far end of a predictor averages everyone beyond the last split point, so extreme values are pulled toward the middle of the range and the curves flatten at both ends. A flat stretch at the edge of a partial dependence plot is often an artefact of the method rather than a plateau in the risk, and a threshold or J-shaped effect at the extremes will be understated. Third, each tree is a step function. A tree's prediction changes only at a split, and the terminal node size (`nodesize`, `r rfsrc_pbc$nodesize` here) controls how far a tree grows: a smaller value allows deeper trees whose leaves hold fewer patients, and so fewer events. Averaging over trees with different split points smooths the steps, so the averaged curve has no fixed step width. A narrow bump or dip can still be the work of a few patients, though, so before reading anything into one, refit the forest with a few different seeds, or on a few random subsamples of patients, and recompute the partial dependence on a common grid. A feature that comes and goes between refits is noise. With few events per node, a larger `nodesize` steadies the estimates, and `randomForestSRC::tune.nodesize()` searches for one on validation error (a holdout set by default at this sample size, out-of-bag error on larger data). In practice we read the shape where the data are dense, and we check the edges against the variable dependence plots above, where each point is a patient and gaps in the data are visible. A rug of the predictor (`geom_rug()` on the patient-level data) does the same job on a partial dependence plot. When the shape matters, a Cox model makes a useful cross-check, provided the comparison is like for like. Fit it with the same covariates as the forest and a spline on the predictor of interest. Then marginalise it the same way: set the predictor to each grid value for every patient, predict survival at the same horizon, and average over the same patients. If the two curves agree in the middle and part ways only in the tails, suspect the forest's boundary behaviour first. A Cox fit on the single predictor answers a different, unadjusted question, so disagreement with it says nothing about the forest's tails. The same caution applies to a forest grown on a single predictor. With one variable, every split is on that variable, so the random selection of candidate variables that decorrelates the trees has nothing to choose from, and the forest reduces to bagged survival trees. There is also nothing to integrate over, so partial dependence and the forest's fitted curve are the same thing. The result is an unadjusted association, much like Kaplan--Meier curves split by that predictor. It can show the shape of the relationship, but it adjusts for nothing. ## Conditional dependence To investigate interactions, we can condition variable dependence on group membership in another variable. Here we examine the dependence of survival on bilirubin, stratified by edema group. ```{r coplot-edema, fig.width=7, fig.height=4, fig.cap="Variable dependence on bilirubin, conditional on edema group."} gg_v1 <- gg_variable(rfsrc_pbc, time = 1) gg_v1$edema <- paste("edema =", gg_v1$edema) plot(gg_v1, xvar = "bili", alpha = 0.5) + labs(y = "Survival at 1 Year", x = st_labs["bili"]) + theme(legend.position = "none") + scale_color_manual(values = event_colors, labels = event_labels) + scale_shape_manual(values = event_marks, labels = event_labels) + facet_grid(~edema) + coord_cartesian(ylim = c(-0.01, 1.01)) ``` The decreasing trend with bilirubin holds across all edema groups, but survival is uniformly lower in the edema = 1 panel, confirming an additive effect. We can also condition on a continuous variable by binning into quantile groups: ```{r coplot-albumin, fig.width=7, fig.height=5, fig.cap="Variable dependence on bilirubin, conditional on albumin groups."} albumin_cts <- quantile_pts(gg_v1$albumin, groups = 6, intervals = TRUE) gg_v1$albumin_grp <- cut(gg_v1$albumin, breaks = albumin_cts) levels(gg_v1$albumin_grp) <- paste("albumin =", levels(gg_v1$albumin_grp)) plot(gg_v1, xvar = "bili", alpha = 0.5) + labs(y = "Survival at 1 Year", x = st_labs["bili"]) + theme(legend.position = "none") + scale_color_manual(values = event_colors, labels = event_labels) + scale_shape_manual(values = event_marks, labels = event_labels) + facet_wrap(~albumin_grp) + coord_cartesian(ylim = c(-0.01, 1.01)) ``` The effect of bilirubin attenuates at higher albumin levels, suggesting an interaction between these two liver function markers. # Partial Dependence Surfaces For a richer view of the interaction between bilirubin and albumin, we construct a partial dependence surface. We compute partial dependence on a grid of 5 albumin values, each evaluated at 25 bilirubin points. ```{r surface-data, error=TRUE, cache=FALSE} # Create grid of albumin values alb_grid <- quantile_pts(pbc_trial$albumin, groups = 5) # For each albumin value, compute partial dependence on bili at ~1 year surface_list <- lapply(alb_grid, function(alb_val) { newx <- pbc_trial[, rfsrc_pbc$xvar.names, drop = FALSE] newx$albumin <- alb_val pd_alb <- gg_partial_rfsrc(rfsrc_pbc, xvar.names = "bili", newx = newx, partial.time = t1yr) df <- pd_alb$continuous df$albumin <- alb_val df }) surface_df <- bind_rows(surface_list) ``` ```{r pd-surface, error=TRUE, fig.height=5, fig.cap="Partial dependence surface: survival at 1 year as a function of bilirubin and albumin. Fill color is the predicted survival probability."} if (!exists("surface_df")) { message("surface_df not available --- skipping surface (see surface-data chunk error above).") } else { ggplot(surface_df, aes(x = x, y = albumin, fill = yhat)) + geom_tile() + scale_fill_viridis_c(name = "Survival") + labs(x = "Bilirubin", y = "Albumin") + theme_bw() } ``` The surface shows that survival is highest when bilirubin is low and albumin is high (upper-left corner), and drops steeply as bilirubin increases or albumin decreases. The curvature of the surface, particularly the steep gradient at low albumin and high bilirubin, confirms the interaction detected in the conditional plots. ## Brier Score and CRPS VIMP and minimal depth tell us which predictors matter, but they say nothing about how well the forest's survival predictions are calibrated. The Brier score answers that question. At each point on the event-time grid it measures the mean squared difference between the predicted survival probability $\hat{S}(t)$ and the binary outcome (survived past $t$ or not). Because right-censored subjects have unknown true outcomes, the calculation uses inverse-probability-of-censoring weighting (IPCW): each subject's squared error is upweighted by the inverse probability of being uncensored at that time, so the estimate remains unbiased even under heavy censoring [@graf:1999]. Two reference values help interpret the curve. A Brier score of 0 is perfect prediction. A score around 0.25 is the uninformative reference: a model that assigns every subject a survival probability of 0.5, regardless of their covariates or the time horizon, scores near 0.25 by construction. That is a reference, not a ceiling; a forest can score above it. A forest below 0.25 at a given time beats that constant guess there. One above it is doing worse than the guess at that time, which tells you the forest is miscalibrated there but not why. `gg_brier()` wraps `randomForestSRC::get.brier.survival()` and returns a tidy data frame ready to plot. ```{r brier-overall, fig.width=7, fig.height=3.5, fig.cap="Time-resolved Brier score for the PBC survival forest."} gg_bs <- gg_brier(rfsrc_pbc) plot(gg_bs) ``` Read the curve left to right. Early on it sits near zero; almost everyone is still alive, so predicting survival is easy. It then climbs for the rest of follow-up and is highest at the last event time, about 11.5 years, at just under 0.2, so it never reaches the 0.25 reference. Deaths accumulate, so the outcome grows less certain, and the censoring weights grow as the at-risk pool thins, so the late estimates rest on fewer subjects. Setting `envelope = TRUE` adds a 15--85% ribbon around that line, showing how spread out the individual subjects' Brier contributions are at each time. A narrow ribbon means the forest is assigning similar risks across subjects; a wide one means the predicted risks are pulling apart. ```{r brier-envelope, fig.width=7, fig.height=3.5, fig.cap="Brier score with 15--85% per-subject envelope."} plot(gg_bs, envelope = TRUE) ``` The running CRPS (continuous ranked probability score) integrates the Brier score over time and divides by the time elapsed since the first event, so each point on the curve is the average of the Brier curve up to that time. It collapses the whole trajectory into one running number, which is handy when you want a single calibration figure rather than a curve to read. Because it is an average of Brier scores it reads on the same scale: near 0 is good, and 0.25 is the constant-predictor reference from above. The value at the right edge of the plot is the average over the whole follow-up. ```{r brier-crps, fig.width=7, fig.height=3.5, fig.cap="Running CRPS for the PBC survival forest."} plot(gg_bs, type = "crps") ``` `gg_brier()` also keeps the integrated CRPS that `randomForestSRC::get.brier.survival()` returns, as the `crps_integrated` attribute. Be careful with it. That number is the raw area under the Brier curve, *not* divided by elapsed time, so it carries the units of the time axis (years here) and grows with the length of follow-up. It is not on the 0 to 0.25 scale, and you cannot compare it across studies with different follow-up. For a single number to report, use the `crps_std` attribute instead. It is upstream's `crps.std`, the raw integral divided by the largest event time, `max(gg_bs$time)`, and it is the value `print()` and `summary()` report: ```{r brier-crps-scalar} attr(gg_bs, "crps_integrated") # area under the Brier curve, in years attr(gg_bs, "crps_std") # divided by max(gg_bs$time) tail(gg_bs$crps, 1) # right edge of the running CRPS plot ``` The last two differ slightly. `crps_std` divides by the largest event time, while the running curve divides by the time elapsed since the first one. # Conclusion We have walked a full random survival forest analysis with **randomForestSRC** and **ggRandomForests**, and the results hang together: - `gg_survival()` drew Kaplan--Meier curves that showed no treatment effect in the PBC trial, the same conclusion reached in @fleming:1991. - `gg_error()` showed the OOB error settling quickly as trees were added. - VIMP (`gg_vimp()`) and minimal depth (`max.subtree()`) landed on the same five continuous predictors (bilirubin, albumin, copper, prothrombin, and age), plus edema and ascites, largely where the proportional hazards model landed. - `gg_variable()` exposed non-proportional hazards for bilirubin and copper, where the gap between time horizons widens. - Partial dependence from `gg_partial_rfsrc()` gave the risk-adjusted version of those curves and backed the log-transforms used in the parametric model. - Conditioning plots and the partial dependence surface drew out the bilirubin--albumin interaction. - `gg_brier()` measured how accurate the predictions were, both across time and as a single CRPS summary. Each `gg_*()` function returns a tidy object (often a data frame, sometimes a small list of data frames); the plotting comes after. Lean on the package's `plot()` methods when the default figure works, and drop down to `ggplot2` directly when it does not. ```{r} #| label: restore-options #| include: false #| cache: false options(old_opts) ``` # References