Modèles de survie en R

Durée de grèves et cycle économique

Pierre Beaucoral

2025-11-28

Code 
knitr::opts_chunk$set(fig.width = 7, fig.height = 4, fig.align = "center")

library(survival)
library(AER)      # StrikeDuration
library(dplyr)
library(ggplot2)

if (requireNamespace("survminer", quietly = TRUE)) {
  library(survminer)
}

# -------------------------
# Données brutes de grève
# -------------------------
data("StrikeDuration")
strikes_raw <- StrikeDuration

# Construction d'un jeu de données pour le TD
# - on crée une censure artificielle à 200 jours
# - on construit 2 versions catégorielles du choc d'activité
strikes <- strikes_raw |>
  mutate(
    # variable de durée
    time   = pmin(duration, 200),
    status = ifelse(duration > 200, 0, 1),  # 1 = fin de grève observée avant 200, 0 = censuré à 200

    # choc d'activité non anticipé (déjà continu) :
    #   uoutput > 0 -> cycle favorable ; < 0 -> cycle défavorable
    uoutput_sign = case_when(
      uoutput < 0  ~ "Cycle défavorable",
      uoutput == 0 ~ "Cycle neutre",
      uoutput > 0  ~ "Cycle favorable"
    ),

    # version en 3 classes par terciles (plutôt pour les KM)
    uoutput_q = cut(
      uoutput,
      breaks = quantile(uoutput, probs = c(0, .33, .66, 1), na.rm = TRUE),
      labels = c("Choc faible", "Choc moyen", "Choc fort"),
      include.lowest = TRUE
    )
  )

# Objet de survie pour les grèves (avec censure artificielle à 200 jours)
S_strikes <- Surv(time = strikes$time, event = strikes$status)

1 1. Notions de base & censure

1.1 Données de durée de vie

  • On s’intéresse à une durée \(T\) (positive) :
    • durée d’une grève
    • durée de chômage
    • durée de séjour à l’hôpital, etc.
  • Pour chaque unité (grève) \(i = 1,\dots,n\) :
    • \(t_i\) : durée observée
    • \(\delta_i\) : indicatrice d’événement
      • \(\delta_i = 1\) si l’événement est observé (fin de grève avant la date limite)
      • \(\delta_i = 0\) si la durée est censurée (on sait juste que la grève dure au moins jusqu’à \(t_i\))

En R, on encode \((t_i, \delta_i)\) avec Surv().

1.2 Fonctions de survie & de risque

  • Fonction de survie : \[S(t) = \Pr(T > t)\]
    Probabilité que la grève soit encore en cours à la date \(t\).

  • Fonction de répartition : \[F(t) = \Pr(T \le t) = 1 - S(t)\]

  • Fonction de risque (hazard) : \[h(t) = \lim_{\Delta t \to 0} \frac{\Pr(t \le T < t + \Delta t \mid T \ge t)}{\Delta t}\]

    → “probabilité instantanée” que la grève se termine juste après \(t\), sachant qu’elle durait encore à \(t\).

1.3 Lien \(S(t)\) – \(h(t)\)

  • En continu, on a : \[S(t) = \exp\left(-\int_0^t h(u)\,du\right)\]
  • Cas simple : risque constant \(h(t) = \lambda\)
    → \(S(t) = \exp(-\lambda t)\) (modèle exponentiel)

1.4 Censure à droite

  • On ne connaît pas toujours la date exacte de fin :
    • fin de la période d’observation alors que la grève continue
    • données administratives tronquées
  • On observe alors :
    • \(t_i = \min(T_i, C_i)\)
    • \(\delta_i = \mathbb{1}(T_i \le C_i)\)
  • Hypothèse clé : censure non informative
    → le mécanisme de censure est indépendant de la durée sous-jacente, conditionnellement aux covariables.

1.5 Censure à droite – Rappel R

Code 
# Objet de survie pour les grèves (avec censure artificielle à 200 jours)
S_strikes <- Surv(time = strikes$time, event = strikes$status)
head(S_strikes)
[1]  7  9 13 14 26 29

2 2. Kaplan–Meier & log-rank dans R

2.1 Estimateur Kaplan–Meier

  • But : estimer \(S(t)\) sans paramètre (non paramétrique)
  • On empile les durées ordonnées où surviennent des événements :
    • à chaque temps \(t_j\), on observe \(d_j\) événements parmi \(n_j\) grèves encore en cours.
  • Estimateur : \[\hat S(t) = \prod_{t_j \le t} \left(1 - \frac{d_j}{n_j}\right)\]

2.2 Kaplan–Meier en R

Code 
km_all <- survfit(S_strikes ~ 1, data = strikes)
km_all
Call: survfit(formula = S_strikes ~ 1, data = strikes)

      n events median 0.95LCL 0.95UCL
[1,] 62     61     27      21      41
Code 
plot(km_all, xlab = "Durée de grève (jours)", ylab = "Probabilité de grève encore en cours",
     main = "Courbe de survie Kaplan–Meier\nToutes les grèves")

Ou avec {survminer} :

Code 
if (requireNamespace("survminer", quietly = TRUE)) {
  ggsurvplot(km_all, conf.int = TRUE,
             xlab = "Durée de grève (jours)",
             ylab = "Probabilité de grève en cours")
}

2.3 Kaplan–Meier par groupes

On veut comparer la durée des grèves selon le choc d’activité :

Code 
km_cycle <- survfit(S_strikes ~ uoutput_q, data = strikes)
summary(km_cycle)
Call: survfit(formula = S_strikes ~ uoutput_q, data = strikes)

                uoutput_q=Choc faible 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    2     25       1     0.96  0.0392      0.88618        1.000
    3     24       1     0.92  0.0543      0.81957        1.000
    5     23       1     0.88  0.0650      0.76141        1.000
   12     22       2     0.80  0.0800      0.65761        0.973
   15     20       1     0.76  0.0854      0.60974        0.947
   17     19       1     0.72  0.0898      0.56386        0.919
   19     18       1     0.68  0.0933      0.51967        0.890
   21     17       2     0.60  0.0980      0.43566        0.826
   27     15       1     0.56  0.0993      0.39563        0.793
   28     14       1     0.52  0.0999      0.35681        0.758
   38     13       1     0.48  0.0999      0.31919        0.722
   42     12       1     0.44  0.0993      0.28275        0.685
   49     11       1     0.40  0.0980      0.24749        0.646
   61     10       1     0.36  0.0960      0.21346        0.607
   72      9       1     0.32  0.0933      0.18071        0.567
   98      8       1     0.28  0.0898      0.14934        0.525
   99      7       1     0.24  0.0854      0.11947        0.482
  104      6       1     0.20  0.0800      0.09132        0.438
  114      5       1     0.16  0.0733      0.06517        0.393
  117      4       1     0.12  0.0650      0.04151        0.347
  152      3       1     0.08  0.0543      0.02117        0.302
  153      2       1     0.04  0.0392      0.00586        0.273

                uoutput_q=Choc moyen 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    2     17       1   0.9412  0.0571      0.83572        1.000
    7     16       1   0.8824  0.0781      0.74175        1.000
    9     15       2   0.7647  0.1029      0.58746        0.995
   13     13       1   0.7059  0.1105      0.51936        0.959
   14     12       1   0.6471  0.1159      0.45548        0.919
   25     11       1   0.5882  0.1194      0.39521        0.876
   26     10       1   0.5294  0.1211      0.33818        0.829
   29      9       1   0.4706  0.1211      0.28423        0.779
   37      8       1   0.4118  0.1194      0.23329        0.727
   41      7       1   0.3529  0.1159      0.18543        0.672
   49      6       1   0.2941  0.1105      0.14083        0.614
   52      5       2   0.1765  0.0925      0.06320        0.493
   85      3       1   0.1176  0.0781      0.03200        0.432
  119      2       1   0.0588  0.0571      0.00879        0.394
  130      1       1   0.0000     NaN           NA           NA

                uoutput_q=Choc fort 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1     20       1     0.95  0.0487       0.8591        1.000
    2     19       1     0.90  0.0671       0.7777        1.000
    3     18       4     0.70  0.1025       0.5254        0.933
    4     14       1     0.65  0.1067       0.4712        0.897
    8     13       1     0.60  0.1095       0.4195        0.858
   10     12       1     0.55  0.1112       0.3700        0.818
   11     11       1     0.50  0.1118       0.3226        0.775
   22     10       1     0.45  0.1112       0.2772        0.731
   23      9       1     0.40  0.1095       0.2339        0.684
   27      8       1     0.35  0.1067       0.1926        0.636
   32      7       1     0.30  0.1025       0.1536        0.586
   33      6       1     0.25  0.0968       0.1170        0.534
   35      5       1     0.20  0.0894       0.0832        0.481
   43      4       2     0.10  0.0671       0.0269        0.372
   44      2       1     0.05  0.0487       0.0074        0.338
  100      1       1     0.00     NaN           NA           NA
Code 
if (requireNamespace("survminer", quietly = TRUE)) {
  ggsurvplot(km_cycle, data = strikes,
             conf.int = FALSE, pval = TRUE,
             legend.title = "Choc d'activité",
             legend.labs  = levels(strikes$uoutput_q),
             xlab = "Durée de grève (jours)",
             ylab = "Probabilité de grève en cours")
}

2.4 Test du log-rank

  • Hypothèse nulle \(H_0\) : même fonction de survie dans tous les groupes
  • Statistique de test \(\chi^2\) :
    • compare, à chaque temps \(t_j\), les événements observés vs attendus dans chaque groupe.

En R :

Code 
survdiff(S_strikes ~ uoutput_q, data = strikes)
Call:
survdiff(formula = S_strikes ~ uoutput_q, data = strikes)

                       N Observed Expected (O-E)^2/E (O-E)^2/V
uoutput_q=Choc faible 25       24     32.4  2.198928  5.143950
uoutput_q=Choc moyen  17       17     16.9  0.000153  0.000219
uoutput_q=Choc fort   20       20     11.6  6.074740  8.183882

 Chisq= 9.2  on 2 degrees of freedom, p= 0.01

  • Le test du log-rank compare les courbes de survie des trois groupes de choc (Choc faible, Choc moyen, Choc fort) sous

    \(H_0\)​ : même fonction de survie pour les trois groupes.

  • La statistique de test \(\chi^2 = 9{,}2\) avec 2 ddl donne une p-valeur de 0,01. → Au seuil de 5 %, on rejette \(H_0\)​ : la durée des grèves dépend significativement du niveau de choc d’activité.

3 3. Modèle de Cox

3.1 Modèle de Cox – rappel

  • Modèle des risques proportionnels : \[h_i(t) = h_0(t)\,\exp(\beta_1 x_{i1} + \dots + \beta_p x_{ip})\]

  • \(h_0(t)\) : baseline hazard (fonction de risque de référence, non paramétrique)

  • Interprétation des coefficients :

    • \(\exp(\beta_k)\) = hazard ratio pour un accroissement d’une unité de \(x_k\)

    • \(> 1\) : grèves plus “rapides” à se terminer (durée plus courte)

    • \(< 1\) : grèves plus longues

3.2 Cox en R – un seul prédicteur

Code 
cox_uoutput <- coxph(S_strikes ~ uoutput, data = strikes)
summary(cox_uoutput)
Call:
coxph(formula = S_strikes ~ uoutput, data = strikes)

  n= 62, number of events= 61 

             coef exp(coef)  se(coef)     z Pr(>|z|)   
uoutput     9.216 10060.629     3.229 2.855  0.00431 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

        exp(coef) exp(-coef) lower .95 upper .95
uoutput     10061   9.94e-05     17.96   5634413

Concordance= 0.608  (se = 0.041 )
Likelihood ratio test= 8.43  on 1 df,   p=0.004
Wald test            = 8.15  on 1 df,   p=0.004
Score (logrank) test = 8.31  on 1 df,   p=0.004
  • coef : estimateur de \(\beta\)
  • exp(coef) : hazard ratio
  • coxph estime \(\beta\) par maximisation de la vraisemblance partielle.

3.3 Cox en R – plusieurs prédicteurs

On peut enrichir avec une version catégorielle du choc :

Code 
cox_multi <- coxph(S_strikes ~ uoutput + uoutput_q, data = strikes)
summary(cox_multi)
Call:
coxph(formula = S_strikes ~ uoutput + uoutput_q, data = strikes)

  n= 62, number of events= 61 

                        coef exp(coef) se(coef)     z Pr(>|z|)
uoutput               5.7304  308.0941  10.2171 0.561    0.575
uoutput_qChoc moyen   0.0587    1.0604   0.6119 0.096    0.924
uoutput_qChoc fort    0.3886    1.4750   1.0638 0.365    0.715

                    exp(coef) exp(-coef) lower .95 upper .95
uoutput               308.094   0.003246 6.192e-07 1.533e+11
uoutput_qChoc moyen     1.060   0.942994 3.196e-01 3.518e+00
uoutput_qChoc fort      1.475   0.677978 1.833e-01 1.187e+01

Concordance= 0.608  (se = 0.041 )
Likelihood ratio test= 8.79  on 3 df,   p=0.03
Wald test            = 9.01  on 3 df,   p=0.03
Score (logrank) test = 9.52  on 3 df,   p=0.02

On comparera les modèles (quanti vs quali) plus tard, dans les exercices, via AIC et tests de rapport de vraisemblance.

3.4 Vérifier l’hypothèse de risques proportionnels

Test “graphique” : courbes de log(-log(S(t))) par groupe.
Si les risques sont proportionnels, les courbes devraient être à peu près parallèles.

Code 
if (requireNamespace("survminer", quietly = TRUE)) {
  ggsurvplot(
    survfit(S_strikes ~ uoutput_q, data = strikes),
    fun = "cloglog",    # log(-log(S))
    xlab = "log(temps)",
    ylab = "log(-log(Survie))",
    legend.title = "Choc d'activité"
  )
}

3.5 Test de Schoenfeld

En pratique, on utilise souvent les résidus de Schoenfeld :

Code 
cox_zph <- cox.zph(cox_uoutput)
cox_zph
          chisq df    p
uoutput 0.00662  1 0.94
GLOBAL  0.00662  1 0.94
Code 
plot(cox_zph)
  • p-valeur petite → violation possible de l’hypothèse de PH pour le prédicteur concerné.

4 4. Modèles paramétriques

4.1 Weibull en risques proportionnels

  • Cas où : \[h(t) = \lambda p t^{p-1}\]
  • Si \(p=1\) → modèle exponentiel
  • Dans un cadre PH avec covariables : \[h_i(t) = \lambda p t^{p-1} \exp(\beta^\top x_i)\]

Intérêt? Plus d’efficacité si la forme est à peu près correcte:

  • Si cette forme est raisonnablement proche de la réalité :

    • estimation des \(\beta\) plus précise (variance plus faible),

    • meilleure performance en petits échantillons,

    • courbes de survie plus lisses, moins “en escalier” que Kaplan–Meier.

En R :

Code 
weib_ph <- survreg(S_strikes ~ uoutput, data = strikes, dist = "weibull")
summary(weib_ph)

Call:
survreg(formula = S_strikes ~ uoutput, data = strikes, dist = "weibull")
               Value Std. Error     z      p
(Intercept)  3.79012    0.13944 27.18 <2e-16
uoutput     -9.67700    3.00825 -3.22 0.0013
Log(scale)   0.00631    0.10180  0.06 0.9506

Scale= 1.01 

Weibull distribution
Loglik(model)= -285.4   Loglik(intercept only)= -290.2
    Chisq= 9.6 on 1 degrees of freedom, p= 0.002 
Number of Newton-Raphson Iterations: 6 
n= 62

Attention : survreg utilise la paramétrisation AFT. On peut revenir à l’écriture PH via des transformations (cf. aide R).

4.2 Modèles AFT (Accelerated Failure Time)

  • On suppose : \[\log T_i = \alpha + \beta^\top x_i + \varepsilon_i\]
  • Selon la loi de \(\varepsilon_i\), on obtient différents modèles :
    • normal → log-normal
    • logistique → log-logistique
    • Gumbel → Weibull (vu comme AFT)

En R, toujours via survreg() :

Code 
aft_lognormal   <- survreg(S_strikes ~ uoutput, data = strikes, dist = "lognormal")
aft_loglogistic <- survreg(S_strikes ~ uoutput, data = strikes, dist = "loglogistic")
summary(aft_lognormal)

Call:
survreg(formula = S_strikes ~ uoutput, data = strikes, dist = "lognormal")
              Value Std. Error     z      p
(Intercept)  3.2155     0.1606 20.02 <2e-16
uoutput     -9.3664     3.3945 -2.76 0.0058
Log(scale)   0.2047     0.0909  2.25 0.0243

Scale= 1.23 

Log Normal distribution
Loglik(model)= -287.4   Loglik(intercept only)= -291
    Chisq= 7.2 on 1 degrees of freedom, p= 0.0073 
Number of Newton-Raphson Iterations: 3 
n= 62
Code 
summary(aft_loglogistic)

Call:
survreg(formula = S_strikes ~ uoutput, data = strikes, dist = "loglogistic")
             Value Std. Error     z      p
(Intercept)  3.295      0.165 19.99 <2e-16
uoutput     -9.615      3.625 -2.65 0.0080
Log(scale)  -0.333      0.105 -3.16 0.0016

Scale= 0.717 

Log logistic distribution
Loglik(model)= -288.6   Loglik(intercept only)= -292
    Chisq= 6.85 on 1 degrees of freedom, p= 0.0088 
Number of Newton-Raphson Iterations: 3 
n= 62

5 5. Application : durée de grèves

5.1 Jeu de données StrikeDuration

  • Données issues de Kennan (1985), “The Duration of Contract Strikes in U.S. Manufacturing”, Journal of Econometrics.
  • 62 grèves dans l’industrie manufacturière US (1968–1976).
  • Variables :
    • duration : durée de la grève (jours)
    • uoutput : choc d’activité non anticipé (mesure macro, liée à la conjoncture)
Code 
summary(strikes_raw)
    duration         uoutput        
 Min.   :  1.00   Min.   :-0.10443  
 1st Qu.: 10.25   1st Qu.:-0.03143  
 Median : 27.00   Median : 0.01138  
 Mean   : 42.68   Mean   : 0.01102  
 3rd Qu.: 51.25   3rd Qu.: 0.06450  
 Max.   :216.00   Max.   : 0.07427

5.2 Construction des variables de TD

  • On introduit :

    • une censure artificielle à 200 jours ;
    Code 
    summary(strikes)
        duration         uoutput              time            status      
 Min.   :  1.00   Min.   :-0.10443   Min.   :  1.00   Min.   :0.0000  
 1st Qu.: 10.25   1st Qu.:-0.03143   1st Qu.: 10.25   1st Qu.:1.0000  
 Median : 27.00   Median : 0.01138   Median : 27.00   Median :1.0000  
 Mean   : 42.68   Mean   : 0.01102   Mean   : 42.42   Mean   :0.9839  
 3rd Qu.: 51.25   3rd Qu.: 0.06450   3rd Qu.: 51.25   3rd Qu.:1.0000  
 Max.   :216.00   Max.   : 0.07427   Max.   :200.00   Max.   :1.0000  
 uoutput_sign             uoutput_q 
 Length:62          Choc faible:25  
 Class :character   Choc moyen :17  
 Mode  :character   Choc fort  :20

  • deux codes pour le choc d’activité :

    • uoutput_sign : défavorable / favorable
    Code 
    table(strikes$uoutput_sign)
    
Cycle défavorable   Cycle favorable 
               25                37
    • uoutput_q : choc faible / moyen / fort (terciles)
Code 
table(strikes$uoutput_q)

Choc faible  Choc moyen   Choc fort 
         25          17          20

6 6. Exercices – Énoncé (survie de grève)

6.1 1. Présentation des données et analyse descriptive

On considère les 62 grèves de la base StrikeDuration.

  1. (Correspondance entre codages du cycle)
    On construit :

    • uoutput_sign : 3 modalités (cycle défavorable, favorable)
    • uoutput_q : 3 modalités (choc faible, moyen, fort, par terciles)
    1. Donner la table de contingence entre ces deux codages et la commenter (cohérence des deux “échelles” du cycle).
    Code 
    with(strikes, table(uoutput_sign, uoutput_q))
                       uoutput_q
uoutput_sign        Choc faible Choc moyen Choc fort
  Cycle défavorable          25          0         0
  Cycle favorable             0         17        20
    Code 
    prop.table(with(strikes, table(uoutput_sign, uoutput_q)), 1)
                       uoutput_q
uoutput_sign        Choc faible Choc moyen Choc fort
  Cycle défavorable   1.0000000  0.0000000 0.0000000
  Cycle favorable     0.0000000  0.4594595 0.5405405

6.2 1. Présentation des données et analyse descriptive (2)

  1. (Corrélations entre variables explicatives)
    On considère comme variables explicatives potentielles :

    • uoutput (continu)
    • une version numérique de uoutput_q (par ex. 1 = faible, 2 = moyen, 3 = fort)
    1. Construire une matrice de corrélations entre ces deux variables.
    2. Commenter : y a-t-il une information nouvelle dans la version catégorielle ?
    Code 
    strikes <- strikes |>
  mutate(uoutput_q_num = as.numeric(uoutput_q))

cor(dplyr::select(strikes, uoutput, uoutput_q_num), use = "complete.obs")
                    uoutput uoutput_q_num
uoutput       1.0000000     0.9245881
uoutput_q_num 0.9245881     1.0000000

6.3 1. Présentation des données et analyse descriptive (3)

  1. (Courbes de survie brutes)
    On trace les courbes de survie Kaplan–Meier de la durée de grève selon uoutput_q.

    1. Que constate-t-on sur les différences de durée de grève selon le choc d’activité ?
    Code 
    km_cycle <- survfit(S_strikes ~ uoutput_q, data = strikes)
plot(km_cycle)

  2. (Test du log-rank)
    On utilise un test du log-rank pour tester l’effet global de uoutput_q.

    1. Écrire les hypothèses \(H_0\) / \(H_1\).
    2. Interpréter la statistique de test et la p-valeur.
    Code 
    survdiff(S_strikes ~ uoutput_q, data = strikes)
    Call:
survdiff(formula = S_strikes ~ uoutput_q, data = strikes)

                       N Observed Expected (O-E)^2/E (O-E)^2/V
uoutput_q=Choc faible 25       24     32.4  2.198928  5.143950
uoutput_q=Choc moyen  17       17     16.9  0.000153  0.000219
uoutput_q=Choc fort   20       20     11.6  6.074740  8.183882

 Chisq= 9.2  on 2 degrees of freedom, p= 0.01

6.4 1. Présentation des données et analyse descriptive (4)

  1. (Vérification graphique des risques proportionnels)
    On trace les courbes de log(-log(S(t))) par niveau de uoutput_q.

    1. Les courbes semblent-elles approximativement parallèles ?
    2. Conclure sur la plausibilité du modèle de Cox avec cette variable.
    Code 
    if (requireNamespace("survminer", quietly = TRUE)) {
  ggsurvplot(
    survfit(S_strikes ~ uoutput_q, data = strikes),
    fun = "cloglog"
  )
}

7 Exercices – Correction (survie de grève)

7.1 Présentation des données et analyse descriptive

7.1.1 Q1 – Correspondance entre codages du cycle

On considère les 62 grèves de la base StrikeDuration.

On a construit :

  • uoutput_sign : 2 modalités (cycle défavorable, favorable)
  • uoutput_q : 3 modalités (choc faible, moyen, fort, par terciles)
Code 
with(strikes, table(uoutput_sign, uoutput_q))
                   uoutput_q
uoutput_sign        Choc faible Choc moyen Choc fort
  Cycle défavorable          25          0         0
  Cycle favorable             0         17        20
Code 
prop.table(with(strikes, table(uoutput_sign, uoutput_q)), 1)
                   uoutput_q
uoutput_sign        Choc faible Choc moyen Choc fort
  Cycle défavorable   1.0000000  0.0000000 0.0000000
  Cycle favorable     0.0000000  0.4594595 0.5405405

Commentaires :

  • Tous les épisodes en cycle défavorable sont classés en Choc faible.
  • Les épisodes en cycle favorable se répartissent entre Choc moyen et Choc fort.
  • La matrice est quasi diagonale : les deux codages racontent la même histoire (choc défavorable vs favorable), l’un en version “signe”, l’autre en version “force du choc”.
  • On peut donc considérer que les deux échelles sont cohérentes.

7.1.2 Q2 – Corrélations entre variables explicatives

Variables explicatives potentielles :

  • uoutput (continu)
  • uoutput_q_num : version numérique de uoutput_q (1 = faible, 2 = moyen, 3 = fort)
Code 
strikes <- strikes |>
  mutate(uoutput_q_num = as.numeric(uoutput_q))

cor(dplyr::select(strikes, uoutput, uoutput_q_num), use = "complete.obs")
                uoutput uoutput_q_num
uoutput       1.0000000     0.9245881
uoutput_q_num 0.9245881     1.0000000

Commentaires :

  • On obtient une corrélation positive forte entre uoutput et uoutput_q_num (la variable en terciles est une discrétisation monotone de uoutput).
  • La version catégorielle n’apporte donc pas d’information entièrement nouvelle, mais elle permet :
    • de capturer plus facilement des effets non linéaires du choc,
    • de produire des comparaisons simples (Choc faible vs Choc fort) en termes de hazard ratios.

7.1.3 Q3 – Courbes de survie brutes (Kaplan–Meier)

Code 
km_cycle <- survfit(S_strikes ~ uoutput_q, data = strikes)

plot(
  km_cycle,
  xlab = "Durée de grève (jours)",
  ylab = "Probabilité de grève en cours",
  col = c("black", "darkgray", "gray70"),   # 3 couleurs
  lty = c(1, 2, 3),                          # 3 types de ligne
  lwd = 2                                    # lignes plus épaisses
)

legend(
  "topright",
  legend = levels(strikes$uoutput_q),
  col    = c("black", "darkgray", "gray70"),
  lty    = c(1, 2, 3),
  lwd    = 2,
  bty    = "n"
)

Commentaires :

  • Pour les chocs forts (conjoncture très favorable), la courbe de survie décroît plus vite : les grèves ont tendance à être plus courtes.
  • Pour les chocs faibles (conjoncture défavorable), la courbe décroît plus lentement : les grèves restent plus longtemps en cours.
  • Les courbes sont bien séparées, ce qui suggère un effet du contexte macroéconomique sur la durée des grèves.

7.1.4 Q4 – Test du log-rank

Code 
survdiff(S_strikes ~ uoutput_q, data = strikes)
Call:
survdiff(formula = S_strikes ~ uoutput_q, data = strikes)

                       N Observed Expected (O-E)^2/E (O-E)^2/V
uoutput_q=Choc faible 25       24     32.4  2.198928  5.143950
uoutput_q=Choc moyen  17       17     16.9  0.000153  0.000219
uoutput_q=Choc fort   20       20     11.6  6.074740  8.183882

 Chisq= 9.2  on 2 degrees of freedom, p= 0.01

Commentaires :

  • Hypothèses du test du log-rank :
    • \(H_0\) : les fonctions de survie sont identiques dans les 3 groupes de choc (Choc faible, Choc moyen, Choc fort).
    • \(H_1\) : au moins une fonction de survie diffère.
  • On obtient une statistique \(\chi^2\) d’environ 9,2 avec 2 ddl et une p-valeur ≈ 0,01.
  • Au seuil de 5 %, on rejette \(H_0\) : la durée des grèves dépend significativement du niveau de choc d’activité.

7.1.5 Q5 – Vérification graphique des risques proportionnels

Code 
if (requireNamespace("survminer", quietly = TRUE)) {
  ggsurvplot(
    survfit(S_strikes ~ uoutput_q, data = strikes),
    fun = "cloglog",
    xlab = "log(temps)",
    ylab = "log(-log(Survie))"
  )
}

Commentaires :

  • Les courbes log(-log(S(t))) par niveau de uoutput_q sont grossièrement quasi parallèles (sans croisements extrêmes).
  • L’hypothèse de risques proportionnels pour uoutput_q est donc raisonnablement plausible dans ce jeu de données.
  • On peut donc utiliser un modèle de Cox avec cette variable, tout en gardant en tête que l’échantillon est petit.

7.2 2. Modèle de Cox – Correction

7.2.1 Q6 – Modèle de Cox continu

Modèle ajusté :

Code 
cox_full <- coxph(S_strikes ~ uoutput + uoutput_q_num, data = strikes)
summary(cox_full)
Call:
coxph(formula = S_strikes ~ uoutput + uoutput_q_num, data = strikes)

  n= 62, number of events= 61 

                  coef exp(coef) se(coef)     z Pr(>|z|)
uoutput         5.5249  250.8571  10.3206 0.535    0.592
uoutput_q_num   0.2007    1.2223   0.5385 0.373    0.709

              exp(coef) exp(-coef) lower .95 upper .95
uoutput         250.857   0.003986 4.117e-07 1.529e+11
uoutput_q_num     1.222   0.818155 4.254e-01 3.512e+00

Concordance= 0.608  (se = 0.041 )
Likelihood ratio test= 8.57  on 2 df,   p=0.01
Wald test            = 8.5  on 2 df,   p=0.01
Score (logrank) test = 8.83  on 2 df,   p=0.01
Code 
AIC(cox_full)
[1] 389.1666

Équation du modèle :

\(h_i(t) = h_0(t)\,\exp\big( \eta_1\,uoutput_i + \eta_2\,uoutput_{q,i} \big)\)

Interprétation du coefficient de uoutput :

  • Le terme \(\exp(\eta_1)\) est le hazard ratio associé à une augmentation d’une unité de uoutput.
  • Comme uoutput est un choc macro petit en amplitude, on interprète plutôt de petites variations (par ex. aller d’un quantile faible à un quantile élevé).
  • Si \(\exp(\eta_1) > 1\), un choc plus favorable est associé à une fin plus rapide de la grève (durées plus courtes).

AIC :

  • On retient la valeur d’AIC affichée et on la comparera aux autres spécifications (M1, M3).
  • Un AIC plus faible signale un meilleur compromis ajustement / complexité.

7.2.2 Q7 – Version purement catégorielle du cycle

Code 
cox_cat <- coxph(S_strikes ~ uoutput_q, data = strikes)
summary(cox_cat)
Call:
coxph(formula = S_strikes ~ uoutput_q, data = strikes)

  n= 62, number of events= 61 

                      coef exp(coef) se(coef)     z Pr(>|z|)   
uoutput_qChoc moyen 0.3556    1.4271   0.3253 1.093  0.27435   
uoutput_qChoc fort  0.9643    2.6229   0.3265 2.954  0.00314 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                    exp(coef) exp(-coef) lower .95 upper .95
uoutput_qChoc moyen     1.427     0.7007    0.7543     2.700
uoutput_qChoc fort      2.623     0.3813    1.3832     4.974

Concordance= 0.604  (se = 0.039 )
Likelihood ratio test= 8.47  on 2 df,   p=0.01
Wald test            = 8.81  on 2 df,   p=0.01
Score (logrank) test = 9.3  on 2 df,   p=0.01

Équation du modèle :

\(h_i(t) = h_0(t)\,\exp\big( \eta_2\,\mathbb{1}\{\,uoutput_{q,i} = \text{Choc moyen}\} + \eta_3\,\mathbb{1}\{\,uoutput_{q,i} = \text{Choc fort}\} \big).\)

avec Choc faible comme catégorie de référence.

Interprétation :

  • \(\exp(\eta_2)\) : hazard ratio “Choc moyen” vs “Choc faible”.
  • \(\exp(\eta_3)\) : hazard ratio “Choc fort” vs “Choc faible”.
  • Si au moins un coefficient est significatif, on conclut à une différence de durée entre au moins deux niveaux de choc.
  • En pratique, avec ce petit échantillon, les intervalles de confiance sont larges → significativité parfois fragile, mais la direction va vers des grèves plus rapides en cas de choc plus fort.

Comparaison des déviances :

  • On lit les log-vraisemblances de cox_full et cox_cat dans les sorties et l’on compare :
    • une déviance plus faible \((-2\log L)\) indique un meilleur ajustement (à nombre de paramètres donné).

7.2.3 Q8 – Nombre de paramètres

Modèles :

  • M1 : S ~ uoutput
    • 1 paramètre de pente ( \(\eta_1\) ).
  • M2 : S ~ uoutput + uoutput_q_num
    • 2 paramètres de pente ( \(\eta_1, \eta_2\) ).
  • M3 : S ~ uoutput_q (facteur 3 modalités)
    • 2 paramètres de pente (deux dummies vs référence : Choc moyen, Choc fort).

Emboîtement :

  • M1 est inclus dans M2 si l’on impose (eta_2 = 0).
  • M3 n’est pas emboîté dans M1/M2 (spécification complètement catégorielle vs continue).

7.2.4 Q9 – Choix de la meilleure mesure du cycle

Code 
cox_M1 <- coxph(S_strikes ~ uoutput, data = strikes)
cox_M2 <- coxph(S_strikes ~ uoutput + uoutput_q, data = strikes)
cox_M3 <- coxph(S_strikes ~ uoutput_q, data = strikes)

anova(cox_M1, cox_M2, test = "LRT")
Analysis of Deviance Table
 Cox model: response is  S_strikes
 Model 1: ~ uoutput
 Model 2: ~ uoutput + uoutput_q
   loglik  Chisq Df Pr(>|Chi|)
1 -192.65                     
2 -192.47 0.3621  2     0.8344
Code 
anova(cox_M1, cox_M3, test = "LRT")
Analysis of Deviance Table
 Cox model: response is  S_strikes
 Model 1: ~ uoutput
 Model 2: ~ uoutput_q
   loglik  Chisq Df Pr(>|Chi|)
1 -192.65                     
2 -192.63 0.0382  1      0.845

Commentaires :

  • LR(M1 vs M2) : si la p-valeur est élevée → l’ajout de uoutput_q_num ne permet pas une amélioration significative → on peut préférer le modèle plus simple M1.
  • LR(M1 vs M3) : si la p-valeur est faible → la version catégorielle du choc (M3) capture un effet non linéaire utile.
  • On combine ces tests avec les AIC pour choisir :
    • soit une spécification continue simple (M1),
    • soit une spécification catégorielle (M3) si les effets par classe sont plus lisibles et mieux ajustés.

7.3 3. Modèles paramétriques – Correction

7.3.1 Q10 – Modèle de Weibull

Code 
weib_best <- survreg(S_strikes ~ uoutput, data = strikes, dist = "weibull")
summary(weib_best)

Call:
survreg(formula = S_strikes ~ uoutput, data = strikes, dist = "weibull")
               Value Std. Error     z      p
(Intercept)  3.79012    0.13944 27.18 <2e-16
uoutput     -9.67700    3.00825 -3.22 0.0013
Log(scale)   0.00631    0.10180  0.06 0.9506

Scale= 1.01 

Weibull distribution
Loglik(model)= -285.4   Loglik(intercept only)= -290.2
    Chisq= 9.6 on 1 degrees of freedom, p= 0.002 
Number of Newton-Raphson Iterations: 6 
n= 62

Nombre de paramètres estimés :

  • Intercept ((\(\alpha\) ))
  • Coefficient de uoutput (( \(\eta_1\) ))
  • Paramètre de scale (lié au paramètre de forme du Weibull)

→ Au total, 3 paramètres dans ce modèle.

Comparaison AIC :

  • On calcule l’AIC du Weibull et on le compare à l’AIC du Cox retenu :
    • AIC(Weibull) plus faible → modèle paramétrique préféré.
    • AIC(Cox) similaire ou plus faible → on peut rester en Cox, plus souple (baseline non paramétrique).

7.3.2 Q11 – Courbes de survie pour différents scénarios

On utilise ici le modèle de Cox simple cox_uoutput pour illustrer les différences de survie selon le choc :

Code 
newdata <- data.frame(
  uoutput = quantile(strikes$uoutput, probs = c(.2, .5, .8))
)

surv_cox <- survfit(cox_uoutput, newdata = newdata)

plot(
  surv_cox,
  xlab = "Durée de grève (jours)",
  ylab = "Probabilité de grève en cours"
)

legend(
  "topright",
  legend = paste0("uoutput = ", round(newdata$uoutput, 3)),
  lty = 1,
  bty = "n"
)

Commentaires :

  • Pour les valeurs de uoutput élevées (choc favorable), la courbe de survie est en dessous des autres : les grèves ont une probabilité plus faible de “survivre” longtemps → elles se terminent plus vite.
  • Pour les valeurs faibles / négatives (choc défavorable), la courbe de survie est au-dessus : les grèves ont une probabilité plus forte de durer.
  • La forme du risque implicite (hazard) augmente avec le temps au début, puis peut se stabiliser, ce qui est cohérent avec un modèle de type Weibull.

7.3.3 Q12 – Interprétation économique

En combinant :

  • les hazard ratios des modèles de Cox,
  • les courbes de survie (Kaplan–Meier et prédictions de modèles),

on peut conclure :

  • En période de conjoncture favorable (choc d’activité positif / fort), les grèves ont plus de chances de se terminer rapidement → hazard plus élevé, survie plus faible.
  • En période de conjoncture défavorable, les grèves ont tendance à durer plus longtemps → hazard plus faible, survie plus élevée.

En termes de négociation :

  • Quand l’activité est bonne, les employeurs ont davantage de marges (les pertes liées à la grève sont plus coûteuses à court terme), ce qui peut favoriser des accords plus rapides.
  • Quand la conjoncture est mauvaise, les coûts d’opportunité sont plus faibles et les rapports de force différents, ce qui peut conduire à des grèves plus longues.