1 Introduction

Ce TP introduit les modèles de durée sous R via les outils classiques du package {survival}.

Nous utilisons les données pbc, issues d’un essai clinique sur la cirrhose biliaire primitive.

À la fin du TP, vous devrez être capables de :

  • construire correctement un objet Surv() (durée + censure) ;
  • tracer et interpréter des courbes de survie non paramétriques (Kaplan–Meier) ;
  • utiliser l’estimateur de Nelson–Aalen pour diagnostiquer la forme du hazard ;
  • estimer et interpréter des modèles paramétriques (exponentiel, Weibull) ;
  • estimer un modèle de Cox (PH) et un modèle AFT, et comparer leurs interprétations ;
Plan de la séance (3h)
  1. Comprendre et coder les données de durée (Surv()).
  2. Estimer et tracer une courbe Kaplan–Meier (et KM “à la main”).
  3. Introduire Nelson–Aalen et le test du log-rank.
  4. (Si temps) Modèles paramétriques et modèles avec covariables (Cox / AFT).

2 Rappels Markdown

Bonnes pratiques RMarkdown
  • Un chunk setup avec include = FALSE pour charger les packages.
  • Un set.seed() unique au début pour rendre vos résultats reproductibles.
  • Commenter vos gros blocs de code pour expliquer ce que vous faites (pas seulement “Exercice 1”).

3 Les principaux packages utilisés en analyse de survie sous R

Cette section présente les packages essentiels que vous utiliserez tout au long du TP.
Ils permettent d’estimer, visualiser et interpréter les modèles de durée.

3.1 {survival} — Le package indispensable

Le package de base pour toute analyse de survie en R.
Il fournit les briques fondamentales pour :

  • représenter les données de durée et de censure ;
  • estimer des courbes de survie non paramétriques (Kaplan–Meier, Nelson–Aalen) ;
  • ajuster des modèles paramétriques et semi-paramétriques (exponentiel, Weibull, Cox).

Dans un TP, c’est souvent le seul package vraiment indispensable : les autres (survminer, ggsurvfit…) servent surtout à améliorer les graphiques ou la présentation.

3.1.1 Les fonctions clés de {survival}

3.1.1.1 Surv() — Construire l’objet de survie

La fonction fondamentale : elle encode pour chaque individu le couple \((t_i, \delta_i)\).

S <- Surv(time = df$time, event = df$status == 2)

Principales variantes :

  • Censure à droite (cas le plus courant)

    Surv(time, event)

  • Données start–stop (processus en comptage)

    Surv(time_start, time_stop, event)

  • Censure par intervalle

    Surv(time_left, time_right, type = "interval2")

Dans l’affichage de Surv, un + indique une observation censurée.

3.1.1.2 survfit() — Kaplan–Meier et hazard cumulé

Permet d’estimer :

  • la courbe de survie Kaplan–Meier ;
  • le hazard cumulé (Nelson–Aalen / Fleming–Harrington) ;
  • des survies par groupe.
# Survie globale
km <- survfit(S ~ 1)

# Survie par groupe (ex. sexe)
km_sex <- survfit(S ~ sex, data = df)

Pour le hazard cumulé (Nelson–Aalen) :

na <- survfit(S ~ 1, type = "fh")  # Fleming–Harrington ~ Nelson–Aalen

Les objets retournés contiennent : - time : temps distincts où un événement a lieu ; - n.risk : nombre à risque ; - n.event : nombre d’événements ; - surv : survie estimée \(\hat S(t)\) ; - std.err, lower, upper : erreurs standard et bornes d’IC.

3.1.1.3 coxph() — Modèle de Cox à risques proportionnels (PH)

Ajuste le modèle semi-paramétrique le plus utilisé en analyse de survie.

cox1 <- coxph(S ~ age + sex + bili, data = df)
summary(cox1)

On obtient notamment :

  • les coefficients \(\hat \eta\) ;
  • les hazard ratios via exp(coef(cox1)) ;
  • des tests (Wald, score, vraisemblance) ;
  • la log-vraisemblance du modèle.

3.1.1.4 cox.zph() — Vérifier l’hypothèse de risques proportionnels

Teste si l’effet des covariables est constant dans le temps (hypothèse PH).

zph <- cox.zph(cox1)
zph

Couplé avec :

ggcoxzph(zph)  # via le package {survminer}

pour un diagnostic graphique.

3.1.2 survreg() — Modèles paramétriques et AFT

Permet d’ajuster des modèles paramétriques ou de type Accelerated Failure Time (AFT).

aft_weib <- survreg(S ~ age + sex + bili,
                    data = df, dist = "weibull")
summary(aft_weib)

Distributions usuelles :

  • "exponential"
  • "weibull"
  • "lognormal"
  • "loglogistic"

Les coefficients de survreg() s’interprètent sur l’échelle de \(\log(T)\), et leur exponentielle donne des time ratios.

3.1.3 survdiff() — Test du log-rank

Compare des courbes de survie entre groupes.

test_lr <- survdiff(S ~ sex, data = df)
test_lr

La sortie fournit une statistique de type \(\chi^2\). Le p-value associé se calcule par :

chi2 <- as.numeric(test_lr$chisq)
pval <- 1 - pchisq(chi2, df = 1)  # si deux groupes

3.2 Résumé des fonctions à connaître

Fonction Rôle principal
Surv() Créer l’objet de survie (durée + censure)
survfit() Estimer KM et hazard cumulé (NA/FH)
coxph() Ajuster un modèle de Cox (PH)
cox.zph() Tester la proportionnalité des risques (PH)
survreg() Modèles paramétriques / AFT (exp, Weibull)
survdiff() Test du log-rank entre groupes

3.3 {survminer} — Courbes de survie claires et lisibles

Ce package est spécialisé dans la visualisation des objets issus de {survival}.

3.3.1 Ce qu’il permet de faire

  • Tracer facilement les courbes Kaplan–Meier.
  • Ajouter des tables de risque automatiquement.
  • Afficher les intervalles de confiance.
  • Visualiser les diagnostics PH du modèle de Cox.

3.3.2 Fonctions à connaître

  • ggsurvplot() : graphique Kaplan–Meier prêt à l’emploi ;
  • ggcoxzph() : graphique du test PH (résidus de Schoenfeld).

3.4 {ggsurvfit} — Alternative moderne pour les courbes de survie

Package léger, pensé pour les étudiants.
Il produit des courbes propres, lisibles et intuitives.

3.4.1 Ce qu’il permet de faire

  • Tracer une courbe Kaplan–Meier simple.
  • Ajouter des intervalles de confiance.
  • Ajouter des indicateurs de censure.
  • Ajouter des tables d’effectifs.

3.5 {cmprsk} — Étudier les risques concurrents

Quand plusieurs événements possibles se font concurrence
(ex : décès vs greffe dans les données PBC),
ce package devient intéressant.

3.5.1 Ce qu’il permet de faire

  • Estimer la probabilité cumulée d’un événement par cause (cuminc()).
  • Ajuster un modèle de Fine–Gray (crr()).

3.6 {gtsummary} — Produire des tableaux de résultats

Très utilisé en biostatistique, ce package permet de transformer les sorties de modèle en tableaux propres et lisibles, au format publication.

3.6.1 Ce qu’il permet de faire

  • Présenter proprement un modèle de Cox.
  • Ajouter les hazard ratios, intervalles de confiance et p-values.
  • Formater automatiquement le tableau.

3.7 Résumé des packages

Package À quoi sert-il ?
survival Modèles de survie (KM, Cox, paramétriques).
survminer Graphiques KM riches et lisibles.
ggsurvfit Courbes KM simples, pédagogiques.
cmprsk Risques concurrents (Fine–Gray).
gtsummary Tableaux professionnels pour modèles de survie.
Mémo en 10 lignes
  • Objet de survie : S <- Surv(time, event)
  • Courbe KM globale : survfit(S ~ 1)
  • KM par groupe : survfit(S ~ groupe)
  • Graphique : ggsurvplot()
  • Hazard cumulé : survfit(S ~ 1, type = “fh”) + ggsurvplot(…, fun = “cumhaz”)
  • Modèle de Cox : coxph(S ~ X1 + X2, data = df)
  • PH ? cox.zph() (+ ggcoxzph())
  • Modèle AFT Weibull : survreg(S ~ X1 + X2, dist = “weibull”)
  • HR = exp(coef(cox))
  • TR = exp(coef(aft))

4 Mise en place

library(survival)
library(survminer)
library(dplyr)
library(ggplot2)

data(pbc)
df <- pbc

5 Ex1 – Exploration des données, codage et censure

Objectif de cours : comprendre le couple \((t_i, \delta_i)\), le rôle de la censure, et comment cela se reflète dans le code Surv().

5.1 Aperçu général des données

head(df)
str(df)
summary(df$time)
table(df$status)
En analyse de survie, chaque individu est résumé par le couple \((t_i,\delta_i)\), où :
  • \(t_i\) = durée d’observation,
  • \(\delta_i = 1\) si l’événement est observé, \(\delta_i = 0\) si l’observation est censurée.

La vraisemblance s’écrit alors :

\[ L(\theta) = \prod_{\delta_i = 1} f(t_i \mid \theta) \times \prod_{\delta_i = 0} S(t_i \mid \theta). \] Si event est mal codé → toute l’analyse est biaisée.

Questions

  1. Quelle est la variable de durée dans les données PBC ?

  2. Quelle variable indique l’événement (décès) ?

  3. Calculez la proportion de censures : 

    mean(df$status == 0)

  4. Combien de modalités distinctes a la variable status ? 

    table(df$status)

    → Interprétez chaque modalité : laquelle correspond à la censure ? Au décès ? À une greffe ?

  5. La variable time contient-elle des valeurs aberrantes ? 

    summary(df$time)

    → Y a-t-il des durées très courtes (< 5 jours) ou très longues (> 4000 jours) ?
    Comment cela pourrait-il affecter une estimation paramétrique ?

  6. La censure semble-t-elle administrative ou liée à un mécanisme plus complexe ? 

    plot(df$time, df$status == 0)

    → Censures concentrées en fin de suivi → censure administrative.
    → Sinon → possible censure informative.

  7. Quelles seraient les conséquences d’un codage incorrect de l’événement ?

    • un décès codé comme censure ;
    • une censure codée comme décès.

    Expliquez l’impact sur la courbe KM et sur un modèle de Cox.

  8. Pourquoi delta_i doit-il être binaire (0/1 ou TRUE/FALSE) ?

    Que se passerait-il si on conservait les valeurs brutes 0 / 1 / 2 dans Surv() ?

  9. Combien d’individus sont observés avant 1000 jours ? 

    sum(df$time < 1000)

    → Expliquez pourquoi cela influence la précision de la KM au début du suivi.

  10. Quelle est la médiane brute des durées (sans censure) ? 

    median(df$time)

    → Pourquoi cette médiane n’est-elle pas une médiane de survie ?

  11. Peut-on détecter rapidement si la censure dépend d’une covariable ? 

    tapply(df$status == 0, df$sex, mean)

    → Si les différences sont fortes, que cela implique-t-il pour une analyse naïve ?

5.2 Objet Surv() et première courbe Kaplan–Meier

5.2.1 Création de l’objet Surv()

S <- Surv(time = df$time,
          event = df$status == 2)  # 2 = décès
S
Rôle de Surv()
L’objet Surv() encode pour chaque individu le couple \((t_i, \delta_i)\) utilisé dans la vraisemblance : \[ L(\theta) = \prod_{\delta_i = 1} f(t_i \mid \theta) \times \prod_{\delta_i = 0} S(t_i \mid \theta). \]
  • \(\delta_i = 1\) (événement observé) → contribution via la densité \(f(t_i)\) ;
  • \(\delta_i = 0\) (censuré) → contribution via la survie \(S(t_i)\).
Dans l’affichage de Surv, un + à droite d’un temps indique une observation censurée.
Questions de compréhension
  1. En regardant S[1:10], repérez :
    • un individu pour lequel l’événement est observé ;
    • un individu pour lequel l’observation est censurée (avec +).
    Pour chacun, notez \((t_i, \delta_i)\).
  2. Pourquoi est-il important de coder event en 0/1 ou TRUE/FALSE plutôt qu’en 0/1/2 ? Indice : que se passerait-il si vous écriviez 
    S_bad <- Surv(time = df$time, event = df$status)
    en gardant les 3 modalités de status ? Que voyez-vous dans S_bad[1:10] ?
  3. Dans les données PBC, nous avons choisi “décès” comme événement d’intérêt. Proposez un autre choix raisonnable d’événement (par exemple “décès ou greffe”) et écrivez le code pour l’objet Surv correspondant.
  4. La définition de S dépend uniquement de deux vecteurs : time et event. Donnez deux exemples de situations (dans un autre contexte que PBC) où :
    • la variable de temps serait un temps d’attente ;
    • l’événement d’intérêt serait différent.

5.2.2 Première courbe KM avec survfit() et ggsurvplot()

km <- survfit(S ~ 1)

ggsurvplot(
  km,
  conf.int   = TRUE,
  risk.table = TRUE,
  censor     = TRUE,
  ggtheme    = theme_minimal(),
  xlab       = "Temps (jours)",
  ylab       = "Probabilité de survie"
)
Pourquoi l’estimateur KM fonctionne ?
L’estimateur de Kaplan–Meier maximise une vraisemblance non paramétrique : \[ \hat S(t) = \prod_{t_i \le t} \left(1 - \frac{d_i}{n_i}\right), \] où \(d_i\) est le nombre d’événements au temps \(t_i\) et \(n_i\) le nombre d’individus à risque juste avant \(t_i\).
Lecture rapide de la KM
  1. Médiane de survie :
    À peu près à quel temps la survie passe-t-elle sous 0,5 (médiane de survie si elle existe) ? Si la courbe ne passe jamais sous 0,5, que peut-on dire de la médiane de survie ?
  2. Rôle des censures :
    Comment les censures (petits traits ou symboles) modifient-elles la forme de la courbe par rapport à une courbe sans censure ? Donnez un exemple de segment de la courbe où l’on observe des censures mais aucun événement.
  3. Table de risque :
    À partir de quel temps la table de risque indique-t-elle qu’il reste très peu d’individus à risque ? Que conclure sur la fiabilité de la forme de la courbe à droite (fin de suivi) ?
  4. Lecture ponctuelle de S(t) :
    Lisez graphiquement \(\hat S(1000)\) et \(\hat S(2000)\). Comment interpréter la différence entre ces deux valeurs en termes de probabilité de survie ?
  5. Lien avec la vraisemblance :
    En vous appuyant sur la formule \[ \hat S(t) = \prod_{t_i \le t}\left(1 - \frac{d_i}{n_i}\right), \] expliquez en une phrase ce que représente le facteur \(1 - d_i / n_i\) au temps \(t_i\). Pourquoi peut-on l’interpréter comme une probabilité conditionnelle de « survivre juste après \(t_i\) » ?
  6. Effet de la troncature de l’axe des temps :
    Modifiez le graphique pour limiter l’axe des abscisses à 3000 jours (argument xlim = c(0, 3000) dans ggsurvplot). La conclusion sur la médiane de survie et la forme globale de la courbe change-t-elle ? Comment le choix de l’échelle de temps peut-il influencer la perception visuelle des résultats ?
  7. Sans intervalle de confiance :
    Reproduisez le graphique en mettant conf.int = FALSE. Est-ce plus « lisible » ? Quelles informations importantes perdez-vous en supprimant les intervalles de confiance ?

6 Ex2 – Construire une courbe de survie « à la main »

Objectif de cours : reconstruire pas à pas l’estimateur de Kaplan–Meier,
comprendre le rôle des quantités \(t_i, d_i, n_i\),
et visualiser comment elles déterminent toute la courbe de survie.

6.1 Sous-échantillon et définition de l’événement

On commence par extraire un sous-échantillon de taille modérée pour garder des tables lisibles.

Questions :

  1. Dans le codage ci-dessus, que représentent exactement les variables event et is_cens ?

6.2 Temps de décès et construction de la table (t_i, d_i, n_i)

On construit maintenant une table de survie où chaque ligne correspond à un temps où au moins un événement (décès) se produit.

On ne garde que les temps où il y a au moins un décès, car ce sont les seuls temps où la survie va “descendre”.

Questions :

  1. Vérifiez que n_i décroît quand le temps augmente. Pourquoi est-ce logique ?
  2. Expliquez pourquoi on doit prendre en compte à la fois les décès et les censures pour calculer n_i.
  3. Que se passerait-il si l’on construisait n_i en ne retirant que les décès et en ignorant les censures ?

6.3 Hazard discret h_i et survie S_i

Le hazard discret au temps \(t_i\) est défini par :

\[ h_i = \frac{d_i}{n_i}. \]

La survie à la date \(t_i\) est obtenue comme produit-limite :

\[ S_i = \prod_{j : t_j \le t_i} (1 - h_j). \]

Questions :

  1. Pourquoi peut-on interpréter \(1 - h_i\) comme une probabilité conditionnelle de “survivre juste après t_i” ?
  2. Pourquoi la suite \(S_i\) est-elle décroissante en \(i\) ?

6.4 Comparaison avec la courbe KM de survfit()

On estime maintenant la courbe de Kaplan–Meier standard sur le même sous-échantillon et on la compare à la courbe reconstruite “à la main”.

S_sub  <- Surv(df_km$time, df_km$event)

km_sub <- survfit(Surv(time, event) ~ 1, data = df_km)

6.4.1 Courbe KM standard

ggsurvplot(
  km_sub,
  conf.int   = FALSE,
  risk.table = TRUE,
  ggtheme    = theme_minimal()
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
##   Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Ignoring unknown labels:
## • colour : "Strata"

6.4.2 Survie “à la main” vs survfit() sur le même graphique

Questions de synthèse :

  1. Les deux courbes coïncident-elles à l’œil nu ? Si vous voyez de très légères différences, d’où peuvent-elles venir (arrondis, gestion des temps égaux, etc.) ?
  2. Pourquoi dit-on que survfit() implémente exactement le produit-limite de Kaplan–Meier à partir des quantités \(t_i, d_i, n_i\) ?
  3. En fin de suivi, lorsque n_i devient petit, la courbe reconstruite vous semble-t-elle plus instable ? Comment l’interprétez-vous en termes d’information disponible ?
  4. Expliquez en quelques phrases en quoi cet exercice vous aide à mieux comprendre ce que représente une courbe de Kaplan–Meier et pourquoi elle est adaptée aux données censurées.

7 Ex3 – Hazard cumulé de Nelson–Aalen (calcul à la main)

Objectif : comprendre le lien entre hazard discret, hazard cumulé et survie,
et utiliser la forme de \(H(t)\) pour diagnostiquer l’évolution du hazard.

7.1 Rappel : hazard cumulé (définition)

Le hazard cumulé estimé par Nelson–Aalen est défini par :

\[ \hat H(t_i) = \sum_{j : t_j \le t_i} \frac{d_j}{n_j}. \]

On peut l’interpréter comme une mesure cumulée du “risque instantané” accumulé dans le temps.

7.2 Calcul de \(H_i\) et de la survie associée

tab <- tab %>%
  dplyr::mutate(
    H_i  = cumsum(d_i / n_i),
    S_na = exp(-H_i)
  )

tab

Questions :

  1. Pourquoi \(H_i\) est-il toujours croissant ?
  2. Que représente exactement la quantité \(\exp(-H_i)\) ?

7.3 Visualiser le hazard cumulé

ggplot(tab, aes(x = time, y = H_i)) +
  geom_step(linewidth = 1) +
  labs(
    title = "Hazard cumulé de Nelson–Aalen (calcul à la main)",
    x = "Temps",
    y = "H(t)"
  ) +
  theme_minimal()

Questions d’interprétation :

  1. La courbe \(H(t)\) est-elle plutôt linéaire, convexe ou concave ?
  2. Que cela suggère-t-il sur la forme du hazard : constant, croissant, décroissant ?

7.4 Comparaison des survies : KM vs Nelson–Aalen

Nelson–Aalen fournit une approximation de la survie :

\[ \tilde S(t) = \exp[-\hat H(t)]. \]

On compare cette approximation à la survie \(S_i\) obtenue via le produit-limite de Kaplan–Meier.

ggplot(tab, aes(x = time)) +
  geom_step(aes(y = S_i),  linewidth = 1) +
  geom_step(aes(y = S_na), linetype = "dashed") +
  labs(
    title = "Survie KM (plein) vs exp(-H_NA) (pointillés)",
    x = "Temps",
    y = "Survie"
  ) +
  theme_minimal()

Questions :

  1. Dans quels cas \(S_{KM}(t)\) et \(\exp[-H(t)]\) coïncident-ils presque parfaitement ?
    (Indice : approximation \(\log(1 - x) \approx -x\) pour \(x\) petit.)
  2. Pourquoi la version Nelson–Aalen tend-elle à légèrement surestimer la survie en fin de suivi ?

7.5 Résumé visuel : lire un hazard cumulé

Lecture d’un hazard cumulé :
  • Courbe à peu près linéaire → hazard environ constant.
  • Courbe convexe → hazard croissant dans le temps.
  • Courbe concave → hazard décroissant.

8 Ex5 – Modèles paramétriques (exponentiel et Weibull)

Objectif de cours : comprendre le lien entre la forme du hazard, la forme de la survie, et les hypothèses implicites des modèles exponentiel et Weibull.

8.1 Modèle exponentiel

m_exp <- survreg(S ~ 1, data = df, dist = "exponential")
summary(m_exp)

# Dans survreg : scale = 1 / lambda
lambda_hat <- 1 / m_exp$scale
lambda_hat
# Superposition KM vs modèle exponentiel
ggsurvplot(km, conf.int = FALSE) +
  stat_function(fun = function(t) exp(-lambda_hat * t))
Questions – Modèle exponentiel
  1. La forme hazard constant vous semble-t-elle compatible avec la forme empirique de \(H(t)\) ?
  2. Comparez la médiane théorique \(\log(2)/\hat\lambda\) à la médiane KM.

8.2 Modèle Weibull

m_weib <- survreg(S ~ 1, data = df, dist = "weibull")
summary(m_weib)

k   <- 1 / m_weib$scale
lam <- exp(-m_weib$coefficients * k)

k
lam
# Comparaison exponentiel vs Weibull
AIC(m_exp, m_weib)
BIC(m_exp, m_weib)
Questions – Modèle Weibull
  1. Que signifie \(k > 1\) / \(k < 1\) en termes de hazard (croissant/décroissant) ?
  2. Le Weibull améliore-t-il nettement l’ajustement (AIC/BIC, médiane) par rapport à l’exponentiel ?

Liens avec les modèles paramétriques :
Si \(H(t)\) est à peu près linéaire, un modèle exponentiel (hazard constant) est plausible.
Si \(H(t)\) est clairement convexe ou concave, un modèle Weibull (hazard croissant ou décroissant) est souvent plus adapté.

9 Ex6 – Modèles PH (Cox) et AFT (Weibull) avec covariables

Objectif de cours : distinguer les deux grands cadres (PH vs AFT) et apprendre à lire HR vs TR.

Rappel de cours – Modèle exponentiel
Sans covariables, le modèle exponentiel suppose : \[ h(t) = \lambda \quad (\text{hazard constant}), \qquad S(t) = \exp(-\lambda t). \] La médiane de survie est : \[ \tilde t_{\text{exp}} = \frac{\log(2)}{\lambda}. \] Un hazard constant implique que le risque instantané ne dépend pas de l’ancienneté (pas de “vieillissement” dans le temps).

9.1 Modèle de Cox (PH)

cox1 <- coxph(S ~ age + sex + bili, data = df)
summary(cox1)

exp(coef(cox1))       # hazard ratios
exp(confint(cox1))    # intervalles de confiance sur HR

Rappel : modèle PH de Cox
\[ h(t \mid X) = h_0(t)\,\exp(X\beta) \] Un HR > 1 signifie un risque instantané plus élevé d’événement ; HR < 1, un risque plus faible.

Questions – Modèle exponentiel
  1. Comparez la courbe KM et la courbe exponentielle théorique.
    • Sont-elles globalement proches ou très différentes ?
    • Les différences sont-elles concentrées au début, au milieu ou en fin de suivi ?
  2. Reprenez le graphique de l’hazard cumulé de Nelson–Aalen (Ex3). La courbe \(\hat H(t)\) vous semble-t-elle à peu près linéaire (hazard constant) ou clairement courbée ? En quoi cela confirme ou infirme la plausibilité du modèle exponentiel ?
  3. Comparez la médiane théorique du modèle exponentiel med_exp à la médiane de survie lue sur la KM (Ex1 / Ex2).
    • L’ordre de grandeur est-il raisonnable ?
    • Un écart important signale quel type de mauvaise spécification ?
  4. Si, en réalité, le hazard est croissant dans le temps, que se passe-t-il si l’on impose un hazard constant ?
    • La courbe exponentielle aura tendance à être trop pessimiste ou trop optimiste au début / à la fin ?
    • Quel risque pour des prédictions à long terme (survie extrapolée) ?

9.2 Vérification de l’hypothèse PH

zph <- cox.zph(cox1)
zph
ggcoxzph(zph)
Lire les graphiques de ggcoxzph() :
  • Ligne horizontale ≈ effet constant dans le temps → compatible avec PH.
  • Tendance clairement croissante ou décroissante → effet qui change dans le temps → violation possible de PH.
  • On regarde la tendance globale, pas les petites oscillations locales.
Questions – PH
  1. Certaines covariables semblent-elles violer PH (p-value faible) ?
  2. Quelle forme prend leur courbe de résidus de Schoenfeld (croissante, décroissante) ?

9.3 Modèle Weibull AFT

Rappel de cours – Modèle Weibull
Sans covariables, le modèle Weibull suppose : \[ h(t) = \lambda k t^{k-1}, \qquad S(t) = \exp(-\lambda t^k), \] où :
  • \(k\) = paramètre de forme
  • \(\lambda\) = paramètre d’échelle
Interprétation :
  • \(k > 1\) : hazard croissant (risque qui augmente avec le temps)
  • \(k = 1\) : hazard constant (cas exponentiel)
  • \(k < 1\) : hazard décroissant
La médiane de survie est : \[ \tilde t_{\text{Weibull}} = \left(\frac{\log(2)}{\lambda}\right)^{1/k}. \] 
m_weib <- survreg(S ~ 1, data = df, dist = "weibull")
summary(m_weib)

# Dans survreg (dist = "weibull") :
# scale = 1 / k
k_hat   <- 1 / m_weib$scale

# coef = log(lambda^{-1/k}) => lambda = exp(-coef * k)
lam_hat <- exp(-m_weib$coefficients * k_hat)

k_hat
lam_hat

# Médiane théorique Weibull
med_weib <- (log(2) / lam_hat)^(1 / k_hat)
med_weib


p_weib <- ggsurvplot(
  km,
  conf.int   = FALSE,
  risk.table = FALSE,
  ggtheme    = theme_minimal(),
  xlab       = "Temps (jours)",
  ylab       = "S(t)"
)

p_weib$plot +
  stat_function(
    fun = function(t) exp(-lam_hat * t^k_hat),
    linetype = "dashed"
  ) +
  labs(title = "Kaplan–Meier (plein) vs Weibull (pointillé)")
Questions – Modèle Weibull
  1. Le paramètre de forme estimé \(\hat k\) est-il supérieur, égal ou inférieur à 1 ? Reliez ce signe à la forme observée de \(\hat H(t)\) (linéaire, convexe, concave).
  2. Comparez visuellement :
    • KM vs exponentiel,
    • KM vs Weibull.
    Quel modèle colle le mieux à la courbe empirique, et sur quels segments de temps (début / milieu / fin) ?
  3. Comparez la médiane théorique Weibull med_weib à :
    • la médiane exponentielle med_exp ;
    • la médiane de Kaplan–Meier.
    Quel modèle semble le plus cohérent avec la médiane non paramétrique ?
  4. Un modèle plus flexible comme Weibull “s’adapte” mieux à la forme de la survie. Pourquoi ne suffit-il pas de regarder la superposition graphique pour conclure qu’il est “meilleur” ? (Indice : pensez à la pénalisation de la complexité du modèle.)

Rappel : modèle AFT
\[ \log(T) = X\gamma + \varepsilon \quad \Leftrightarrow \quad T = \exp(X\gamma)\, T_0. \] TR > 1 : durée plus longue ; TR < 1 : durée plus courte.

# Comparaison AIC
AIC(cox1)
AIC(aft_weib)
Rappel – AIC / BIC
\[ \text{AIC} = -2 \ell(\hat\theta) + 2k, \qquad \text{BIC} = -2 \ell(\hat\theta) + k \log(n), \] où \(k\) est le nombre de paramètres, \(n\) la taille de l’échantillon. Plus AIC/BIC sont faibles, meilleur est le compromis ajustement / complexité.
Questions – Choix de modèle
  1. Lequel des deux modèles minimise l’AIC ? Le BIC ? La différence est-elle faible (< 2), modérée (≈ 2–5) ou forte (> 5) ?
  2. Comment combiner :
    • les critères AIC/BIC,
    • la comparaison visuelle des courbes de survie,
    • et la plausibilité clinique/économique de la forme du hazard (croissant, décroissant, constant) ?
    pour justifier un choix de modèle auprès d’un non-spécialiste ?
  3. Donnez un exemple de contexte (hors PBC) où un hazard croissant est naturel (ex. panne d’un moteur, faillite d’entreprise), et un contexte où un hazard décroissant serait plus plausible (ex. abandon d’un programme après une phase initiale).
À retenir sur les modèles paramétriques
  • La forme de \(\hat H(t)\) (linéaire / convexe / concave) est un outil clé pour choisir un modèle.
  • L’exponentiel est un cas particulier du Weibull (\(k = 1\)).
  • Un bon modèle paramétrique doit être compatible avec :
    • les diagnostics non paramétriques (KM, Nelson–Aalen),
    • les critères d’information (AIC/BIC),
    • et la logique économique/clinique du phénomène étudié.

10 Conclusion – Checklist de compétences

Auto-checklist – Cochez ce que vous savez faire seul·e :
  • [ ] Créer un objet Surv() correctement codé.
  • [ ] Tracer et commenter une courbe de Kaplan–Meier.
  • [ ] Lire la forme du hazard (constant / croissant / décroissant) sur un hazard cumulé.
  • [ ] Estimer et comparer un modèle exponentiel et Weibull (AIC, médiane théorique, forme du hazard).
  • [ ] Estimer un modèle de Cox, interpréter les HR et vérifier l’hypothèse PH.
  • [ ] Estimer un modèle AFT (Weibull), interpréter les time ratios et comparer à Cox.