De sept prédicteurs candidats à un modèle de Cox parcimonieux — screening univarié, AIC pas-à-pas, tests du rapport de vraisemblance, vérification d’une interaction, la concordance et un forest plot ggforest
Construisez et sélectionnez un modèle de Cox à risques proportionnels multivariable en R. Faites un screening univarié des covariables candidates, comparez les stratégies de sélection — connaissance clinique, AIC pas-à-pas (MASS::stepAIC) et LASSO pénalisé — comparez des modèles emboîtés avec un test du rapport de vraisemblance, testez une interaction traitement-par-covariable, lisez la concordance (C-index) et tracez le forest plot des hazard ratios ggforest de survminer, sur les données de cancer du poumon NCCTG.
coxph() une covariable à la fois) pour voir quels prédicteurs portent un signal par eux-mêmes — mais ne décidez jamais du modèle final sur le seul screening.MASS::stepAIC(), piloté par les données) ou pénalisé (LASSO via glmnet, pour de nombreux prédicteurs). L’AIC arbitre entre ajustement et complexité (plus bas = mieux) ; le pas-à-pas est piloté par les données et instable — traitez-le comme une suggestion.anova(reduced, full)) et testez les interactions de la même façon — un terme traitement-par-covariable significatif signifie que l’effet diffère selon les patients (un facteur prédictif, et pas seulement pronostique).summary(fit)$concordance (0.5 = hasard, plus haut = meilleure discrimination) ; utilisez rms::validate() pour l’optimisme par bootstrap. Tracez ensuite le forest plot des hazard ratios ggforest() emblématique.Vous avez le tableur d’un clinicien : 168 patients atteints d’un cancer du poumon, le délai jusqu’au décès et sept prédicteurs candidats — âge, sexe, deux scores de statut de performance, perte de poids, apport calorique. Vous savez déjà ajuster un modèle de Cox unique et lire un hazard ratio. La vraie question est celle qui vient ensuite : lesquels de ces sept appartiennent au modèle, et comment décider ?
Mettez les sept et vous risquez un modèle surajusté qui ne se reproduira pas. Gardez-en trop peu et vous manquez un vrai facteur pronostique ou vous confondez l’effet qui vous intéresse. Cette leçon est le workflow pratique entre ces deux erreurs — les étapes qu’un chercheur suit réellement pour passer d’une pile de variables candidates à un seul modèle de Cox multivariable parcimonieux, interprétable et validé.
Nous travaillons sur le jeu de données NCCTG lung-cancer (lung, du package survival) et visualisons avec survminer. Vous allez :
MASS::stepAIC) et une note sur la régression pénalisée (LASSO) ;ggforest() final.Pré-spécifié (réglementaire / confirmatoire) : vous fixez la liste des variables à partir de la connaissance clinique avant de voir les données, de sorte qu’il n’y a aucun biais de sélection à défendre. Exploratoire (construction d’un modèle pronostique) : vous laissez les données aider à choisir, à l’aide de l’AIC ou de la pénalisation — et vous devez alors valider, car une sélection pilotée par les données gonfle la performance apparente. La plupart des analyses sont un hybride : un noyau clinique (traitement, principaux facteurs pronostiques) que vous gardez toujours, plus quelques variables choisies systématiquement et validées.
Nous utilisons lung, le jeu de données du North Central Cancer Treatment Group fourni avec survival (228 patients atteints d’un cancer du poumon avancé). Plusieurs covariables ont des valeurs manquantes, et cela compte pour la construction du modèle (voir plus bas), donc nous construisons d’abord un jeu de données en cas complets sur les sept candidats et la réponse. Le remodelage et la sélection ici utilisent base R uniquement :
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 Nous conservons 168 patients avec 121 décès — contre 228, car meal.cal à elle seule est manquante pour 47 patients. Les colonnes candidates :
time délai de survie (jours) ; status 1 = censuré, 2 = décédé (l’événement).sex 1 = homme, 2 = femme ; age années.ph.ecog score de performance ECOG (0 = actif … plus haut = altéré) ; ph.karno score de Karnofsky évalué par le médecin (0–100, plus haut = mieux) ; pat.karno score de Karnofsky évalué par le patient.wt.loss perte de poids (livres, 6 derniers mois) ; meal.cal calories consommées aux repas.Comme les valeurs manquantes diffèrent d’une covariable à l’autre, ajuster le modèle complet et un modèle réduit sur des sous-ensembles de cas complets différents les rend non comparables — un test du rapport de vraisemblance n’aurait aucun sens. Fixez l’ensemble d’analyse une seule fois (lung_cc ci-dessus) et ajustez chaque modèle candidat dessus. C’est pourquoi nous avons préparé les données d’abord.
Avant de construire quoi que ce soit, examinez chaque candidat isolément : ajustez un modèle de Cox univarié par covariable et tabulez le hazard ratio, son IC à 95 % et la p-value. C’est un balayage, pas une décision — il vous indique quelles variables portent un signal et à peu près sa force. Bouclez sur les noms avec 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.9700À lire : sex (HR 0.59, protecteur), ph.ecog (HR 1.6, néfaste) et les deux scores de Karnofsky (HR juste en dessous de 1 — meilleure performance, hazard plus faible) ressortent comme les signaux les plus forts au screening ; age est limite, et wt.loss et meal.cal ne portent presque aucun signal univarié ici. Une règle courante consiste à conserver tout ce qui a un p de screening < 0.20 (un seuil généreux — une variable peut compter dans le modèle multivariable même si elle paraît faible seule), mais ne laissez jamais le seul screening choisir votre modèle : il ignore comment les variables s’expliquent les unes les autres, et il gonfle les faux positifs. Le screening resserre le champ ; les étapes suivantes choisissent.
Le choix par défaut conforme aux exigences réglementaires : inclure les variables qu’un expert du domaine désigne comme pertinentes, fixées à l’avance. Ici c’est le modèle complet sur les sept candidats — chaque effet ajusté sur les six autres :
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-04Lecture du modèle complet : après ajustement mutuel, sex (HR 0.58, p = 0.006) et ph.ecog (HR 2.08, p = 0.001) restent clairement significatifs ; ph.karno est limite (p = 0.046), wt.loss et pat.karno sont suggestifs (p ≈ 0.06–0.12) ; age et meal.cal n’apportent essentiellement rien (p = 0.36 et p = 0.90). La concordance du modèle est de 0.65 et son AIC est de 1011.5 — des chiffres que nous tenterons d’améliorer en élaguant le poids mort.
stepAIC)Lorsque vous n’avez pas de liste fixe, laissez un critère d’information choisir. L’AIC (Akaike Information Criterion) note un modèle comme l’ajustement pénalisé par le nombre de paramètres — plus bas est mieux — donc il récompense le fait d’expliquer les données mais punit l’ajout de variables qui ne méritent pas leur place. MASS::stepAIC() parcourt le modèle dans les deux directions (ajout et retrait de termes) jusqu’à ce que l’AIC cesse de s’améliorer :
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 À lire : l’AIC pas-à-pas retire age et meal.cal et garde sex, ph.ecog, ph.karno, pat.karno, wt.loss — abaissant l’AIC de 1011.5 à 1008.4. (Mettez trace = TRUE pour observer chaque étape d’ajout/retrait.)
stepAIC() est pratique, mais elle choisit les variables en regardant le résultat, donc ses p-values et ses intervalles de confiance sont optimistes (ils ignorent que les variables ont été sélectionnées), et un petit changement dans les données peut faire basculer quelles variables survivent. Traitez le modèle pas-à-pas comme une suggestion à confronter à la connaissance clinique, pas comme une réponse finale — et si vous le rapportez, validez-le (ci-dessous). Le BIC (k = log(n) dans stepAIC) pénalise la complexité plus lourdement et produit des modèles plus parcimonieux ; l’AIC privilégie la prédiction, le BIC privilégie le « vrai » modèle plus petit.
Avec de nombreux prédicteurs candidats (des dizaines de biomarqueurs, un ratio EPV en dessous de ~10 événements par variable), la sélection pas-à-pas s’effondre et surajuste. La régression pénalisée est l’alternative moderne : le LASSO (glmnet avec family = "cox") rétrécit les coefficients vers zéro et met les plus faibles exactement à zéro, réalisant la sélection et l’estimation en une seule étape, avec la force de la pénalité choisie par validation croisée (cv.glmnet). Il gère élégamment les prédicteurs corrélés (là où le pas-à-pas en choisit un arbitrairement) et ne nécessite aucun test séquentiel. Pour un problème à 7 variables comme le nôtre c’est surdimensionné — mais pour un problème de grande dimension c’est l’outil adéquat. La mécanique (construire une matrice de modèle, cv.glmnet(x, y, family = "cox"), lire les coefficients à lambda.1se) est couverte dans les leçons de sélection de modèle de ML supervisé ; utilisez-le quand les prédicteurs sont plus nombreux que les événements, pas avant.
Le pas-à-pas a donné un modèle à cinq variables ; un raisonnement clinique pourrait préférer juste sex + ph.ecog. Lequel est meilleur ? Lorsqu’un modèle est emboîté dans un autre (les termes du modèle plus petit sont un sous-ensemble de ceux du plus grand), comparez-les avec un test du rapport de vraisemblance (RV) via anova() — il demande si les variables supplémentaires améliorent l’ajustement au-delà du hasard :
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 ' ' 1À lire : ajouter les trois termes Karnofsky/perte de poids donne χ² = 7.92 sur 3 df, p = 0.048 — une amélioration tout juste significative de l’ajustement. Le modèle plus grand s’ajuste donc significativement (quoique marginalement) mieux que sex + ph.ecog seul. Le test du rapport de vraisemblance est la façon rigoureuse de décider si un bloc de variables vaut la peine d’être conservé ; c’est aussi exactement ce que la sélection pas-à-pas automatise une variable à la fois.
Les deux comparent des modèles, différemment. Le test du rapport de vraisemblance donne une p-value pour « les termes supplémentaires aident-ils ? » — mais seulement pour des modèles emboîtés. L’AIC donne un nombre que vous pouvez classer entre modèles qu’ils soient emboîtés ou non (le plus bas gagne), avec une pénalité de complexité intégrée. Utilisez le test du rapport de vraisemblance pour justifier un bloc précis de variables ; utilisez l’AIC/BIC pour comparer plusieurs modèles candidats à la fois.
Jusqu’ici chaque effet est pronostique — il déplace le risque de la même façon pour tout le monde. Une interaction pose une question différente : l’effet d’une variable dépend-il d’une autre ? Dans un essai, une interaction traitement × biomarqueur signifierait que le traitement aide certains patients plus que d’autres — un facteur prédictif qui oriente la médecine de précision. Testez-la de la même façon : ajustez le modèle avec l’interaction (a * b se développe en a + b + a:b) et confrontez-le par test du rapport de vraisemblance au modèle sans :
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.3317À lire : l’interaction sex:ph.ecog donne χ² = 0.94 sur 1 df, p = 0.33 — non significative. L’ effet du score de performance ne diffère pas selon le sexe (et inversement), donc il n’y a aucune interaction à rapporter : les deux restent de simples facteurs pronostiques. Gardez le modèle additif. Si cela avait été significatif, vous rapporteriez l’effet au sein de chaque sous-groupe (p. ex. le HR de ph.ecog séparément chez les hommes et les femmes) plutôt qu’un unique HR moyenné.
Un facteur pronostique prédit le résultat indépendamment du traitement (il entre dans le modèle comme effet principal). Un facteur prédictif modifie l’effet du traitement (il apparaît comme une interaction traitement × facteur significative) — la base pour choisir qui traiter. Pré-spécifiez toujours les interactions que vous testez ; les pêcher à l’aveugle gonfle les faux positifs.
La sélection n’est pas la validation. Avant de faire confiance au modèle, demandez-vous à quel point il discrimine — étant donné deux patients, à quelle fréquence classe-t-il celui qui meurt le plus tôt comme le plus à risque ? C’est la concordance de Harrell (C-index), rapportée par summary() : 0.5 = pile ou face, 1.0 = parfait, ~0.7+ = utile. Prenez le modèle final à cinq variables :
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 À lire : C = 0.66 (SE 0.03) — le modèle classe une paire aléatoire de patients dans le bon ordre de survie environ deux fois sur trois, modestement mieux que le hasard. C’est une performance honnête, typique pour des données de survie cliniques avec une poignée de prédicteurs.
Calculé sur les mêmes données utilisées pour ajuster et sélectionner, le C-index est optimiste : il mesure en partie le bruit que le modèle a mémorisé. Pour une estimation honnête, estimez l’optimisme par bootstrap avec rms::validate() — réajustez le modèle sur des rééchantillons bootstrap, mesurez la performance sur les données d’origine et soustrayez l’écart moyen. Le C-index corrigé de l’optimisme est celui que vous rapportez. Avec une sélection lourde (pas-à-pas/LASSO), placez la procédure de sélection entière à l’intérieur du bootstrap, sinon la correction sous-estime l’optimisme.
ggforest)Rassemblez le tout : rapportez les hazard ratios du modèle final en langage clair, puis tracez le forest plot ggforest() emblématique — le HR, l’IC à 95 % et la p-value de chaque covariable sur un seul panneau prêt à publier, avec la ligne de référence à HR = 1 (gauche = protecteur, droite = néfaste) :
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)
Interprétez le modèle final, facteur par facteur :
sex : HR ≈ 0.57 (IC à 95 % 0.39–0.85) — être une femme réduit le hazard de décès d’environ 43 %, l’ effet le plus fort du modèle et clairement significatif.ph.ecog : HR ≈ 2.09 (IC à 95 % 1.34–3.26) — chaque aggravation d’un point du score ECOG double à peu près le hazard. Un fort mauvais facteur pronostique.ph.karno : HR ≈ 1.02 par point — étonnamment au-dessus de 1 une fois ph.ecog dans le modèle (les deux scores de performance sont corrélés, donc les effets ajustés peuvent paraître contre-intuitifs) ; limite (p ≈ 0.07).pat.karno : HR ≈ 0.99 par point — meilleure performance évaluée par le patient, hazard légèrement plus faible (p ≈ 0.11).wt.loss : HR ≈ 0.99 par livre — un effet petit et limite (p ≈ 0.06) ; contre-intuitivement protecteur, une bizarrerie connue de ce jeu de données qu’il vaut mieux signaler que surinterpréter.Les deux termes de Karnofsky et la perte de poids sont limites et corrélés à ph.ecog — un rappel qu’un modèle « sélectionné » est un point de départ pour le jugement clinique, pas le dernier mot.
Un modèle de Cox à risques proportionnels multivariable a été construit sur les données de cancer du poumon NCCTG (n = 168 cas complets, 121 décès). Les prédicteurs candidats ont été soumis à un screening univarié et sélectionnés par AIC pas-à-pas (
MASS::stepAIC, dans les deux directions), retenant le sexe, le score de performance ECOG, les scores de Karnofsky évalués par le médecin et par le patient, et la perte de poids (AIC 1011.5 → 1008.4). Dans le modèle final, le sexe féminin était associé à un hazard de décès plus faible (HR = 0.57, IC à 95 % 0.39–0.85, p = 0.005) et un score ECOG plus élevé à un hazard plus élevé (HR = 2.09, IC à 95 % 1.34–3.26, p = 0.001), en ajustant sur les scores de Karnofsky et la perte de poids. Une interaction sexe-par-ECOG n’était pas significative (RV χ²(1) = 0.94, p = 0.33). La discrimination du modèle était modeste (C de Harrell = 0.66, SE 0.03) ; une correction de l’optimisme par bootstrap est recommandée avant un usage externe.
Les deux reposent sur la log-vraisemblance partielle maximisée du modèle \(\ell\). L’AIC arbitre entre ajustement et complexité :
\[ AIC = -2\,\ell + 2k \]
où \(k\) est le nombre de coefficients estimés. Ajouter une variable inutile augmente \(\ell\) seulement légèrement mais ajoute \(2\) à la pénalité, donc l’AIC monte — c’est pourquoi il favorise la parcimonie. Le BIC remplace la pénalité par \(k\log(n)\), plus sévère pour les grands \(n\).
Le test du rapport de vraisemblance compare deux modèles emboîtés de log-vraisemblances \(\ell_0\) (réduit) et \(\ell_1\) (complet) :
\[ G = 2\,(\ell_1 - \ell_0) \;\sim\; \chi^2_{\,df} \]
avec \(df\) = le nombre de paramètres supplémentaires. Un grand \(G\) (petite p-value) indique que les termes supplémentaires améliorent l’ajustement au-delà du hasard. L’AIC et le test du rapport de vraisemblance s’accordent dans l’esprit ; l’AIC intègre simplement la pénalité et fonctionne aussi pour les modèles non emboîtés.
glmnet, nombreux prédicteurs / EPV faible).anova(reduced, full)) sur des modèles emboîtés.a * b) et confrontez-la par test du rapport de vraisemblance.summary(fit)$concordance), corrigée de l’optimisme avec rms::validate().Construisez le modèle vous-même — changez la liste des candidats, essayez stepAIC avec le BIC (k = log(nrow(lung_cc))), ou testez une autre interaction. Le bac à sable démarre au premier Run.
Demandez à Prova « aide-moi à construire un modèle de Cox multivariable à partir de mes prédicteurs candidats — fais un screening, lance la sélection AIC pas-à-pas, teste une interaction, et donne-moi la concordance et un ggforest plot » — elle répond avec du code survival + survminer que vous pouvez exécuter sur vos propres données. The runtime is the judge. Demander à Prova →
Les valeurs manquantes changent silencieusement votre échantillon — et cassent la comparaison de modèles. coxph() utilise les cas complets, donc ajuster le modèle complet et un modèle réduit sur des données ayant des patrons de valeurs manquantes différents signifie qu’ils ont été ajustés sur des patients différents — un test du rapport de vraisemblance ou une comparaison d’AIC entre eux est invalide. Construisez un jeu d’analyse en cas complets fixe d’abord (lung[complete.cases(lung[, keep]), ]) et ajustez chaque modèle candidat dessus, ou imputez les valeurs manquantes en amont. Vérifiez fit$n pour voir combien de patients le modèle a réellement utilisés.
La sélection pas-à-pas donne des p-values trop confiantes. Comme les variables ont été choisies en regardant le résultat, les intervalles de confiance du modèle final sont trop étroits et sa performance trop flatteuse. Ne rapportez pas les p-values du pas-à-pas comme si le modèle était pré-spécifié — confrontez la sélection à la connaissance clinique, et validez le C-index avec une correction de l’optimisme par bootstrap (rms::validate).
Les prédicteurs corrélés font paraître les HR ajustés étranges. ph.karno et pat.karno (et ph.ecog) mesurent des choses qui se recoupent, donc une fois l’un dans le modèle un autre peut changer de signe ou perdre sa significativité — c’est de la colinéarité, pas une erreur de données. Quand deux variables sont quasi identiques, gardez la plus cliniquement pertinente (ou combinez-les) plutôt que les deux ; le LASSO gère cela automatiquement.
Ajustez le modèle complet avec coxph(), puis passez-le à MASS::stepAIC() : step_model <- stepAIC(full, direction = "both", trace = FALSE). Il ajoute et retire des termes pour minimiser l’AIC et renvoie le modèle sélectionné. Utilisez direction = "backward" ou "forward" pour une sélection à sens unique, et k = log(n) pour le BIC (un modèle plus parcimonieux). Rappelez-vous que la sélection pas-à-pas est pilotée par les données — ses p-values sont optimistes, donc validez le résultat.
L’AIC (Akaike Information Criterion) note un modèle comme son ajustement pénalisé par le nombre de paramètres : AIC plus bas = mieux. Il vous permet de classer des modèles concurrents — même non emboîtés — tout en décourageant les variables inutiles. Obtenez-le avec AIC(fit) ; stepAIC() automatise la recherche. Le BIC (BIC(fit)) pénalise plus la complexité et favorise les modèles plus petits. Utilisez l’AIC/BIC pour comparer des modèles entiers, et un test du rapport de vraisemblance pour justifier un bloc précis de variables.
Les deux, dans l’ordre. Le screening univarié (une covariable à la fois) est un balayage rapide pour voir quels prédicteurs portent un signal — utilisez un seuil généreux (p. ex. p < 0.20) pour décider quoi conserver, car une variable peut compter dans le modèle complet même si elle paraît faible seule. Choisissez ensuite le modèle multivariable par connaissance clinique, AIC pas-à-pas ou LASSO. Ne choisissez jamais votre modèle final sur les seules p-values univariées : elles ignorent comment les variables s’expliquent les unes les autres.
Si un modèle est emboîté dans l’autre (ses termes sont un sous-ensemble), utilisez un test du rapport de vraisemblance : anova(reduced, full) renvoie un χ² et une p-value indiquant si les variables supplémentaires améliorent l’ajustement. Pour des modèles non emboîtés, comparez AIC(model1) vs AIC(model2) (le plus bas gagne) — l’AIC fonctionne aussi pour les modèles emboîtés. Pour comparer la performance prédictive, comparez la concordance (summary(fit)$concordance), idéalement corrigée de l’optimisme.
Un facteur pronostique prédit le résultat indépendamment du traitement — il entre dans le modèle de Cox comme effet principal (p. ex. le score ECOG). Un facteur prédictif modifie l’effet du traitement — il apparaît comme une interaction traitement × facteur significative, ce qui signifie que le traitement aide certains patients plus que d’autres (la base de la médecine de précision). Testez-le en ajustant l’interaction (treatment * factor) et en la confrontant par test du rapport de vraisemblance au modèle sans ; pré-spécifiez les interactions que vous testez.
stepAIC avec la pénalité BIC (k = log(nrow(lung_cc))). Remplissez le blanc. Le BIC garde-t-il plus ou moins de variables que l’AIC, et lesquelles survivent ?anova() : le plus grand ajoute 2 variables et donne χ² = 11.4 sur 2 df, p = 0.003. En clair, que concluez-vous, et garderiez-vous le plus grand modèle ?Remplissez le blanc avec log(nrow(lung_cc)) — cela transforme la pénalité de stepAIC d’AIC (\(2k\)) en BIC (\(k\log
n\)). Comme le BIC pénalise plus lourdement chaque variable supplémentaire, il tend à garder moins de variables que l’AIC. Comparez les termes retenus avec ceux du modèle AIC (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.Pour la question 2 : les deux variables supplémentaires donnent χ² = 11.4, df = 2, p = 0.003 < 0.05, donc elles améliorent significativement l’ajustement du modèle au-delà du hasard. Gardez le plus grand modèle — le test du rapport de vraisemblance dit que le bloc ajouté porte une information réelle. (Confirmez tout de même que les variables supplémentaires sont cliniquement sensées et pas du simple surajustement de bruit.)
L’AIC pas-à-pas sélectionne un modèle à cinq variables dont le C-index intra-échantillon est de 0.74, mais la correction de l’optimisme par bootstrap (rms::validate) le ramène à 0.66. Qu’est-ce que cela vous apprend ?
Le modèle est modérément optimiste — sa performance honnête est C ≈ 0.66, pas 0.74. L’écart de 0.08 est l’ optimisme : le C-index intra-échantillon a en partie mesuré le bruit que la procédure pas-à-pas a ajusté. Le 0.66 corrigé est ce que vous rapportez et ce que le modèle livrerait sur de nouveaux patients. Un écart de cette taille est courant avec une sélection pilotée par les données ; c’est la raison pour laquelle vous validez plutôt que de faire confiance au chiffre intra-échantillon.
Vous savez désormais construire et sélectionner un modèle de Cox multivariable en R comme un workflow, pas une seule commande : screener les candidats univariément, sélectionner par connaissance clinique ou AIC pas-à-pas (MASS::stepAIC) ou LASSO, comparer des modèles emboîtés avec un test du rapport de vraisemblance (anova), tester les interactions pour séparer les facteurs prédictifs des pronostiques, et valider avec la concordance (C-index), corrigée de l’optimisme via rms::validate. Sur les données de cancer du poumon, cela a réduit sept candidats à un modèle parcimonieux à cinq variables — sexe féminin protecteur, plus mauvais score ECOG néfaste — tracé en un forest plot ggforest() prêt à publier. La discipline compte plus que n’importe quelle fonction : un modèle parcimonieux, justifié et validé tiendra là où un modèle tout-inclus ne tiendra pas. La dernière case à cocher est celle que la sélection ne peut pas : confirmer que les hazard ratios sont constants au fil du temps.
ggforest) ; la fondation sur laquelle cette leçon s’appuie. · Tester l’hypothèse des risques proportionnels — la vérification suivante essentielle : vos HR sélectionnés sont-ils constants au fil du temps ? · Test du log-rank — la comparaison univariée de groupes que Cox généralise. · Qu’est-ce que l’analyse de survie ? — la censure, le hazard et les fonctions de survie expliqués.Prouvez que vous savez le faire. Maîtrisez toute la série Analyse de survie en R — suivez votre parcours, construisez des projets et obtenez un certificat.
Passez à Pro — Prova illimité sur vos propres données et un certificat vérifiable qui atteste la compétence.
dès 15 $/mois facturé annuellement
✓ Vous êtes Pro — continuez. The runtime is the judge.
Prêt à passer à la vitesse supérieure ?
Chaque résultat de cette page a été produit par le code montré, exécuté au moment du build contre un environnement R figé — modifiez n’importe quel bloc et faites Run pour le reproduire vous-même.
@online{2026,
author = {},
title = {Sélection de variables Cox en R : modèle multivariable},
date = {2026-06-26},
url = {https://www.datanovia.com/learn/biostatistics/survival-analysis/cox-model-building-selection},
langid = {fr}
}