From seven candidate predictors to one parsimonious Cox model — univariate screening, stepwise AIC, likelihood-ratio tests, an interaction check, the concordance, and a ggforest forest plot
Build and select a multivariable Cox proportional hazards model in R. Screen candidate covariates univariately, compare selection strategies — clinical knowledge, stepwise AIC (MASS::stepAIC), and penalized LASSO — compare nested models with a likelihood-ratio test, test a treatment-by-covariate interaction, read the concordance (C-index), and draw the survminer ggforest hazard-ratio plot, on the NCCTG lung-cancer data.
coxph() one covariate at a time) to see which predictors carry signal on their own — but never decide the final model on screening alone.MASS::stepAIC(), data-driven), or penalized (LASSO via glmnet, for many predictors). AIC trades fit against complexity (lower = better); stepwise is data-driven and unstable — treat it as a suggestion.anova(reduced, full)) and test interactions the same way — a significant treatment-by-covariate term means the effect differs across patients (a predictive, not just prognostic, factor).summary(fit)$concordance (0.5 = chance, higher = better discrimination); use rms::validate() for bootstrap optimism. Then draw the signature ggforest() hazard-ratio plot.You have a clinician’s spreadsheet: 168 lung-cancer patients, time to death, and seven candidate predictors — age, sex, two performance-status scores, weight loss, calorie intake. You already know how to fit a single Cox model and read a hazard ratio. The real question is the one that comes next: which of those seven belong in the model, and how do you decide?
Throw all seven in and you risk an overfitted model that won’t replicate. Keep too few and you miss a real prognostic factor or confound the effect you care about. This lesson is the practical workflow between those two mistakes — the steps a researcher actually follows to go from a pile of candidate variables to one parsimonious, interpretable, validated multivariable Cox model.
We work on the NCCTG lung-cancer dataset (lung, from the survival package) and visualize with survminer. You will:
MASS::stepAIC), and a note on penalized regression (LASSO);ggforest() hazard-ratio plot.Pre-specified (regulatory / confirmatory): you fix the variable list from clinical knowledge before seeing the data, so there is no selection bias to defend. Exploratory (prognostic-model building): you let the data help choose, using AIC or penalization — and then you must validate, because data-driven selection inflates the apparent performance. Most analyses are a hybrid: a clinical core (treatment, key prognostics) you always keep, plus a few variables chosen systematically and validated.
We use lung, the North Central Cancer Treatment Group dataset that ships with survival (228 advanced lung-cancer patients). Several covariates have missing values, and that matters for model building (more on this below), so we first build a complete-case dataset on the seven candidates and the response. Reshaping and selection here use base R only:
library(survival)
# Candidate predictors + the response, dropping rows with any NA among them (complete cases)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
c(rows = nrow(lung_cc), events = sum(lung_cc$status == 2)) rows events
168 121 We keep 168 patients with 121 deaths — down from 228, because meal.cal alone is missing for 47 patients. The candidate columns:
time survival time (days); status 1 = censored, 2 = dead (the event).sex 1 = male, 2 = female; age years.ph.ecog ECOG performance score (0 = active … higher = impaired); ph.karno physician-rated Karnofsky score (0–100, higher = better); pat.karno patient-rated Karnofsky score.wt.loss weight loss (lbs, last 6 months); meal.cal calories consumed at meals.Because missing values differ across covariates, fitting the full model and a reduced model on different complete-case subsets makes them non-comparable — a likelihood-ratio test would be meaningless. Fix the analysis set once (lung_cc above) and fit every candidate model on it. This is why we prepared the data first.
Before building anything, look at each candidate on its own: fit a univariate Cox model per covariate and tabulate the hazard ratio, its 95% CI, and the p-value. This is a scan, not a decision — it tells you which variables carry signal and roughly how strong. Loop over the names with base R:
library(survival)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
covariates <- c("age", "sex", "ph.ecog", "ph.karno", "pat.karno", "wt.loss", "meal.cal")
screen <- lapply(covariates, function(x) {
m <- coxph(as.formula(paste("Surv(time, status) ~", x)), data = lung_cc)
s <- summary(m)
data.frame(
covariate = x,
HR = round(s$conf.int[1], 2),
lower = round(s$conf.int[3], 2),
upper = round(s$conf.int[4], 2),
p.value = signif(s$coef[5], 2)
)
})
screen <- do.call(rbind, screen)
screen[order(screen$p.value), ] covariate HR lower upper p.value
3 ph.ecog 1.60 1.23 2.07 0.0004
5 pat.karno 0.98 0.97 0.99 0.0022
2 sex 0.62 0.42 0.91 0.0140
1 age 1.02 1.00 1.04 0.0630
4 ph.karno 0.99 0.97 1.00 0.0660
7 meal.cal 1.00 1.00 1.00 0.6100
6 wt.loss 1.00 0.99 1.01 0.9700Read it: sex (HR 0.59, protective), ph.ecog (HR 1.6, harmful), and the two Karnofsky scores (HR just under 1 — higher performance, lower hazard) screen as the strongest signals; age is borderline, and wt.loss and meal.cal carry almost no univariate signal here. A common rule is to carry forward anything with screening p < 0.20 (a generous threshold — a variable can matter in the multivariable model even if it looks weak alone), but never let screening alone pick your model: it ignores how variables explain each other, and it inflates false positives. Screening narrows the field; the next steps choose.
The regulator-friendly default: include the variables a domain expert says belong, fixed in advance. Here that is the full model on all seven candidates — each effect adjusted for the other six:
library(survival)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
full <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno +
pat.karno + wt.loss + meal.cal, data = lung_cc)
summary(full)Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog + ph.karno +
pat.karno + wt.loss + meal.cal, data = lung_cc)
n= 168, number of events= 121
coef exp(coef) se(coef) z Pr(>|z|)
age 1.065e-02 1.011e+00 1.161e-02 0.917 0.35906
sex -5.509e-01 5.765e-01 2.008e-01 -2.743 0.00609 **
ph.ecog 7.342e-01 2.084e+00 2.233e-01 3.288 0.00101 **
ph.karno 2.246e-02 1.023e+00 1.124e-02 1.998 0.04574 *
pat.karno -1.242e-02 9.877e-01 8.054e-03 -1.542 0.12316
wt.loss -1.433e-02 9.858e-01 7.771e-03 -1.844 0.06518 .
meal.cal 3.329e-05 1.000e+00 2.595e-04 0.128 0.89791
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0107 0.9894 0.9880 1.0340
sex 0.5765 1.7347 0.3889 0.8545
ph.ecog 2.0838 0.4799 1.3452 3.2277
ph.karno 1.0227 0.9778 1.0004 1.0455
pat.karno 0.9877 1.0125 0.9722 1.0034
wt.loss 0.9858 1.0144 0.9709 1.0009
meal.cal 1.0000 1.0000 0.9995 1.0005
Concordance= 0.651 (se = 0.029 )
Likelihood ratio test= 28.33 on 7 df, p=2e-04
Wald test = 27.58 on 7 df, p=3e-04
Score (logrank) test = 28.41 on 7 df, p=2e-04Read the full model: after mutual adjustment, sex (HR 0.58, p = 0.006) and ph.ecog (HR 2.08, p = 0.001) stay clearly significant; ph.karno is borderline (p = 0.046), wt.loss and pat.karno are suggestive (p ≈ 0.06–0.12); age and meal.cal add essentially nothing (p = 0.36 and p = 0.90). The model’s concordance is 0.65 and its AIC is 1011.5 — numbers we will try to improve by trimming the dead weight.
stepAIC)When you don’t have a fixed list, let an information criterion choose. AIC (Akaike Information Criterion) scores a model as fit penalized by the number of parameters — lower is better — so it rewards explaining the data but punishes adding variables that don’t earn their place. MASS::stepAIC() walks the model in both directions (add and drop terms) until AIC stops improving:
library(survival)
library(MASS)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
full <- coxph(Surv(time, status) ~ age + sex + ph.ecog + ph.karno +
pat.karno + wt.loss + meal.cal, data = lung_cc)
step_model <- stepAIC(full, direction = "both", trace = FALSE)
step_modelCall:
coxph(formula = Surv(time, status) ~ sex + ph.ecog + ph.karno +
pat.karno + wt.loss, data = lung_cc)
coef exp(coef) se(coef) z p
sex -0.556256 0.573352 0.198616 -2.801 0.00510
ph.ecog 0.738091 2.091939 0.225670 3.271 0.00107
ph.karno 0.020412 1.020621 0.011083 1.842 0.06551
pat.karno -0.012682 0.987398 0.007926 -1.600 0.10960
wt.loss -0.014601 0.985505 0.007702 -1.896 0.05799
Likelihood ratio test=27.47 on 5 df, p=4.616e-05
n= 168, number of events= 121 Read it: stepwise AIC drops age and meal.cal and keeps sex, ph.ecog, ph.karno, pat.karno, wt.loss — lowering AIC from 1011.5 to 1008.4. (Set trace = TRUE to watch each add/drop step.)
stepAIC() is convenient, but it chooses variables by looking at the outcome, so its p-values and confidence intervals are optimistic (they ignore that the variables were selected), and a small change in the data can flip which variables survive. Treat the stepwise model as a suggestion to sanity-check against clinical knowledge, not a final answer — and if you report it, validate it (below). BIC (k = log(n) in stepAIC) penalizes complexity more heavily and yields sparser models; AIC favours prediction, BIC favours the “true” smaller model.
With many candidate predictors (dozens of biomarkers, an EPV ratio below ~10 events per variable), stepwise selection breaks down and overfits. Penalized regression is the modern alternative: LASSO (glmnet with family = "cox") shrinks coefficients toward zero and sets the weakest exactly to zero, doing selection and estimation in one step, with the penalty strength chosen by cross-validation (cv.glmnet). It handles correlated predictors gracefully (where stepwise picks one arbitrarily) and needs no sequential testing. For a 7-variable problem like ours it is overkill — but for a high-dimensional one it is the right tool. The mechanics (build a model matrix, cv.glmnet(x, y, family = "cox"), read coefficients at lambda.1se) are covered in the supervised-ML model-selection lessons; reach for it when predictors outnumber events, not before.
Stepwise gave a five-variable model; clinical reasoning might prefer just sex + ph.ecog. Which is better? When one model is nested inside another (the smaller model’s terms are a subset of the larger’s), compare them with a likelihood-ratio (LR) test via anova() — it asks whether the extra variables improve fit by more than chance:
library(survival)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
reduced <- coxph(Surv(time, status) ~ sex + ph.ecog, data = lung_cc)
fuller <- coxph(Surv(time, status) ~ sex + ph.ecog + ph.karno +
pat.karno + wt.loss, data = lung_cc)
anova(reduced, fuller)Analysis of Deviance Table
Cox model: response is Surv(time, status)
Model 1: ~ sex + ph.ecog
Model 2: ~ sex + ph.ecog + ph.karno + pat.karno + wt.loss
loglik Chisq Df Pr(>|Chi|)
1 -503.14
2 -499.18 7.9242 3 0.04761 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Read it: adding the three Karnofsky/weight-loss terms gives χ² = 7.92 on 3 df, p = 0.048 — a just significant improvement in fit. So the larger model fits significantly (if marginally) better than sex + ph.ecog alone. The LR test is the principled way to decide whether a block of variables is worth keeping; it is also exactly what stepwise selection automates one variable at a time.
Both compare models, differently. The LR test gives a p-value for “do the extra terms help?” — but only for nested models. AIC gives a number you can rank across models whether or not they are nested (lower wins), with a built-in complexity penalty. Use the LR test to justify a specific block of variables; use AIC/BIC to compare several candidate models at once.
So far every effect is prognostic — it shifts risk the same way for everyone. An interaction asks a different question: does one variable’s effect depend on another? In a trial, a treatment × biomarker interaction would mean the treatment helps some patients more than others — a predictive factor that drives precision medicine. Test it the same way: fit the model with the interaction (a * b expands to a + b + a:b) and LR-test it against the model without:
library(survival)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
no_int <- coxph(Surv(time, status) ~ sex + ph.ecog, data = lung_cc)
with_int <- coxph(Surv(time, status) ~ sex * ph.ecog, data = lung_cc)
anova(no_int, with_int)Analysis of Deviance Table
Cox model: response is Surv(time, status)
Model 1: ~ sex + ph.ecog
Model 2: ~ sex * ph.ecog
loglik Chisq Df Pr(>|Chi|)
1 -503.14
2 -502.67 0.9421 1 0.3317Read it: the sex:ph.ecog interaction gives χ² = 0.94 on 1 df, p = 0.33 — not significant. The effect of performance score does not differ by sex (and vice versa), so there is no interaction to report: both stay simple prognostic factors. Keep the additive model. Had this been significant, you would report the effect within each subgroup (e.g. the HR for ph.ecog separately in men and women) rather than a single averaged HR.
A prognostic factor predicts outcome regardless of treatment (it goes in the model as a main effect). A predictive factor modifies the treatment effect (it shows up as a significant treatment × factor interaction) — the basis for choosing who to treat. Always pre-specify the interactions you test; fishing for them inflates false positives.
Selection is not validation. Before trusting the model, ask how well it discriminates — given two patients, how often does it rank the one who dies sooner as higher-risk? That is Harrell’s concordance (C-index), reported by summary(): 0.5 = coin flip, 1.0 = perfect, ~0.7+ = useful. Take the final five-variable model:
library(survival)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
final <- coxph(Surv(time, status) ~ sex + ph.ecog + ph.karno +
pat.karno + wt.loss, data = lung_cc)
summary(final)$concordance C se(C)
0.65568638 0.02896871 Read it: C = 0.66 (SE 0.03) — the model ranks a random pair of patients in the correct survival order about two-thirds of the time, modestly better than chance. That is honest, typical performance for clinical survival data with a handful of predictors.
Computed on the same data used to fit and select, the C-index is optimistic: it partly measures noise the model memorized. For an honest estimate, bootstrap the optimism with rms::validate() — refit the model on bootstrap resamples, measure performance on the original data, and subtract the average gap. The optimism-corrected C-index is what you report. With heavy selection (stepwise/LASSO), put the whole selection procedure inside the bootstrap, or the correction understates the optimism.
ggforest)Bring it together: report the final model’s hazard ratios in plain language, then draw the signature ggforest() plot — every covariate’s HR, 95% CI, and p-value on one publication-ready panel, with the reference line at HR = 1 (left = protective, right = harmful):
library(survival)
library(survminer)
keep <- c("time", "status", "age", "sex", "ph.ecog",
"ph.karno", "pat.karno", "wt.loss", "meal.cal")
lung_cc <- lung[complete.cases(lung[, keep]), ]
final <- coxph(Surv(time, status) ~ sex + ph.ecog + ph.karno +
pat.karno + wt.loss, data = lung_cc)
ggforest(final, data = lung_cc)
Interpret the final model, factor by factor:
sex: HR ≈ 0.57 (95% CI 0.39–0.85) — being female cuts the hazard of death by ~43%, the strongest effect in the model and clearly significant.ph.ecog: HR ≈ 2.09 (95% CI 1.34–3.26) — each one-point worsening of the ECOG score roughly doubles the hazard. A strong bad prognostic factor.ph.karno: HR ≈ 1.02 per point — surprisingly above 1 once ph.ecog is in the model (the two performance scores are correlated, so adjusted effects can look counter-intuitive); borderline (p ≈ 0.07).pat.karno: HR ≈ 0.99 per point — higher patient-rated performance, slightly lower hazard (p ≈ 0.11).wt.loss: HR ≈ 0.99 per lb — a small, borderline effect (p ≈ 0.06); counter-intuitively protective, a known quirk of this dataset worth flagging rather than over-interpreting.The two Karnofsky terms and weight loss are borderline and correlated with ph.ecog — a reminder that a “selected” model is a starting point for clinical judgement, not the final word.
A multivariable Cox proportional hazards model was built on the NCCTG lung-cancer data (n = 168 complete cases, 121 deaths). Candidate predictors were screened univariately and selected by stepwise AIC (
MASS::stepAIC, both directions), retaining sex, ECOG performance score, physician- and patient-rated Karnofsky scores, and weight loss (AIC 1011.5 → 1008.4). In the final model, female sex was associated with a lower hazard of death (HR = 0.57, 95% CI 0.39–0.85, p = 0.005) and a higher ECOG score with a higher hazard (HR = 2.09, 95% CI 1.34–3.26, p = 0.001), adjusting for the Karnofsky scores and weight loss. A sex-by-ECOG interaction was not significant (LR χ²(1) = 0.94, p = 0.33). The model’s discrimination was modest (Harrell’s C = 0.66, SE 0.03); bootstrap optimism correction is recommended before external use.
Both are built on the model’s maximized partial log-likelihood \(\ell\). The AIC trades fit against complexity:
\[ AIC = -2\,\ell + 2k \]
where \(k\) is the number of estimated coefficients. Adding a useless variable raises \(\ell\) only slightly but adds \(2\) to the penalty, so AIC rises — which is why it favours parsimony. BIC replaces the penalty with \(k\log(n)\), harsher for large \(n\).
The likelihood-ratio test compares two nested models with log-likelihoods \(\ell_0\) (reduced) and \(\ell_1\) (full):
\[ G = 2\,(\ell_1 - \ell_0) \;\sim\; \chi^2_{\,df} \]
with \(df\) = the number of extra parameters. A large \(G\) (small p-value) says the extra terms improve fit beyond chance. AIC and the LR test agree in spirit; AIC just bakes the penalty in and works for non-nested models too.
glmnet, many predictors / low EPV).anova(reduced, full)) on nested models.a * b) and LR-test it.summary(fit)$concordance), optimism-corrected with rms::validate().Build the model yourself — change the candidate list, try stepAIC with BIC (k = log(nrow(lung_cc))), or test a different interaction. The sandbox boots on first Run.
Ask Prova “help me build a multivariable Cox model from my candidate predictors — screen them, run stepwise AIC selection, test an interaction, and give me the concordance and a ggforest plot” — it answers with survival + survminer code you can run on your own data. The runtime is the judge. Ask Prova →
Missing values silently change your sample — and break model comparison. coxph() uses complete cases, so fitting the full model and a reduced model on data with different missing patterns means they were fit on different patients — a likelihood-ratio test or AIC comparison between them is invalid. Build one fixed complete-case analysis set first (lung[complete.cases(lung[, keep]), ]) and fit every candidate model on it, or impute the missing values up front. Check fit$n to see how many patients the model actually used.
Stepwise selection gives over-confident p-values. Because the variables were chosen by looking at the outcome, the final model’s confidence intervals are too narrow and its performance too rosy. Don’t report stepwise p-values as if the model were pre-specified — sanity-check the selection against clinical knowledge, and validate the C-index with bootstrap optimism correction (rms::validate).
Correlated predictors make adjusted HRs look strange. ph.karno and pat.karno (and ph.ecog) measure overlapping things, so once one is in the model another can flip sign or lose significance — that is collinearity, not a data error. When two variables are near-duplicates, keep the more clinically meaningful one (or combine them) rather than both; LASSO handles this automatically.
Fit the full model with coxph(), then pass it to MASS::stepAIC(): step_model <- stepAIC(full, direction = "both", trace = FALSE). It adds and drops terms to minimize AIC and returns the selected model. Use direction = "backward" or "forward" for one-way selection, and k = log(n) for BIC (a sparser model). Remember stepwise selection is data-driven — its p-values are optimistic, so validate the result.
AIC (Akaike Information Criterion) scores a model as its fit penalized by the number of parameters: lower AIC = better. It lets you rank competing models — even non-nested ones — while discouraging useless variables. Get it with AIC(fit); stepAIC() automates the search. BIC (BIC(fit)) penalizes complexity more and favours smaller models. Use AIC/BIC to compare whole models, and a likelihood-ratio test to justify a specific block of variables.
Both, in order. Univariate screening (one covariate at a time) is a quick scan to see which predictors carry signal — use a generous threshold (e.g. p < 0.20) to decide what to carry forward, because a variable can matter in the full model even if it looks weak alone. Then choose the multivariable model by clinical knowledge, stepwise AIC, or LASSO. Never pick your final model on univariate p-values alone: they ignore how variables explain each other.
If one model is nested in the other (its terms are a subset), use a likelihood-ratio test: anova(reduced, full) returns a χ² and p-value for whether the extra variables improve fit. For non-nested models, compare AIC(model1) vs AIC(model2) (lower wins) — AIC also works for nested models. To compare predictive performance, compare the concordance (summary(fit)$concordance), ideally optimism-corrected.
A prognostic factor predicts outcome regardless of treatment — it enters the Cox model as a main effect (e.g. ECOG score). A predictive factor modifies the treatment effect — it shows up as a significant treatment × factor interaction, meaning the treatment helps some patients more than others (the basis for precision medicine). Test it by fitting the interaction (treatment * factor) and LR-testing it against the model without; pre-specify which interactions you test.
stepAIC with the BIC penalty (k = log(nrow(lung_cc))). Fill the blank. Does BIC keep more or fewer variables than AIC did, and which ones survive?anova(): the larger adds 2 variables and gives χ² = 11.4 on 2 df, p = 0.003. In plain words, what do you conclude, and would you keep the larger model?Fill the blank with log(nrow(lung_cc)) — that turns stepAIC’s penalty from AIC (\(2k\)) into BIC (\(k\log
n\)). Because BIC penalizes each extra variable more heavily, it tends to keep fewer variables than AIC. Compare the retained terms with the AIC model’s (sex, ph.ecog, ph.karno, pat.karno, wt.loss).
bic_model <- stepAIC(full, direction = "both", trace = FALSE,
k = log(nrow(lung_cc)))
bic_model
#> BIC keeps FEWER variables than AIC — typically just the strongest
#> signals (sex and ph.ecog), because its heavier complexity penalty
#> drops the borderline Karnofsky/weight-loss terms the AIC model retained.For question 2: the extra two variables give χ² = 11.4, df = 2, p = 0.003 < 0.05, so they significantly improve the model’s fit beyond chance. Keep the larger model — the LR test says the added block carries real information. (Still confirm the extra variables are clinically sensible and not just over-fitting noise.)
Stepwise AIC selects a five-variable model whose in-sample C-index is 0.74, but bootstrap optimism correction (rms::validate) brings it down to 0.66. What does this tell you?
The model is moderately optimistic — its honest performance is C ≈ 0.66, not 0.74. The 0.08 gap is the optimism: the in-sample C-index partly measured noise the stepwise procedure fit. The corrected 0.66 is what you report and what the model would deliver on new patients. A gap this size is common with data-driven selection; it is the reason you validate rather than trust the in-sample number.
You can now build and select a multivariable Cox model in R as a workflow, not a single command: screen candidates univariately, select by clinical knowledge or stepwise AIC (MASS::stepAIC) or LASSO, compare nested models with a likelihood-ratio test (anova), test interactions to separate predictive from prognostic factors, and validate with the concordance (C-index), optimism-corrected via rms::validate. On the lung-cancer data, that took seven candidates down to a parsimonious five-variable model — female sex protective, worse ECOG score harmful — drawn as a publication-ready ggforest() plot. The discipline matters more than any one function: a model that is parsimonious, justified, and validated will hold up where an everything-in model won’t. The last box to tick is the one selection can’t: confirm the hazard ratios are constant over time.
@online{2026,
author = {},
title = {Cox {Model} {Variable} {Selection} in {R:} {Build} a
{Multivariable} {Model}},
date = {2026-06-26},
url = {https://www.datanovia.com/learn/biostatistics/survival-analysis/cox-model-building-selection},
langid = {en}
}