Objectif — Question 1

Cette partie répond à la Question 1 : quel serait l’effet d’une exposition initiale universelle (\(a_0 = 1\) pour tous) par rapport à une absence totale d’exposition (\(a_0 = 0\) pour tous) ?

On souhaite estimer les courbes de survie contrefactuelles :

\[S^{a_0=1}(t) = P(T^{a_0=1} > t) \quad \text{et} \quad S^{a_0=0}(t) = P(T^{a_0=0} > t)\]

La méthode utilisée ici est la G-computation (standardisation par régression) : on modélise le résultat \(Y\) en fonction de l’exposition et des covariables, puis on prédit les devenirs contrefactuels pour l’ensemble de la population.

Note

Simplification pédagogique : pour cette partie, le critère de jugement est traité comme binaire (décès en fin de suivi, oui/non). La G-computation avec un vrai modèle de survie est possible mais plus complexe. La Partie 2 traitera correctement le format de survie via l’IPTW.

Note sur la session R

Les variables suivantes sont disponibles dans vos chunks interactifs (construites dans le contexte de setup) :

  • df : base de données au format long, avec A0 (exposition initiale) et L0 (facteur de confusion initial)

Étape 0 — Créer la base baseline

La G-computation s’applique sur une base une ligne par individu, avec les valeurs initiales des covariables et le statut final.

Créez df_base contenant : id, A0, L0, X, time (temps total de suivi = dernière valeur de T.stop), et D (statut final = dernière valeur de D).

Montrez moi comment faire !
df_base <- df |>
  group_by(id) |>
  summarise(
    A0   = first(A),
    L0   = first(L),
    X    = first(X),
    time = last(T.stop),
    D    = last(D)
  )

head(df_base)
nrow(df_base)

La base comporte 2000 lignes — une par individu. La variable D indique si l’individu est décédé avant la fin du suivi (3 ans).

Étape 1 — Estimer les modèles de résultat

La G-computation ajuste le modèle de résultat (\(D\)) séparément selon le groupe d’exposition \(A_0\), en contrôlant pour les facteurs de confusion \(X\) et \(L_0\).

Estimez deux modèles de régression logistique sur D : mod1 sur les individus avec A0 = 1, mod0 sur les individus avec A0 = 0.

Montrez moi comment faire !
mod1 <- glm(D ~ X + L0,
            data   = df_base[df_base$A0 == 1, ],
            family = binomial)

mod0 <- glm(D ~ X + L0,
            data   = df_base[df_base$A0 == 0, ],
            family = binomial)

summary(mod1)
summary(mod0)

Étape 2 — Prédire les devenirs contrefactuels

Pour chaque individu (quelle que soit son exposition réelle), on prédit le risque de décès comme s’il avait été exposé (via mod1) et comme s’il n’avait pas été exposé (via mod0).

Calculez les vecteurs y1 et y0 : probabilités de décès prédites sous \(a_0=1\) et \(a_0=0\).

Montrez moi comment faire !
## Probabilité de décès si tout le monde avait A0 = 1
y1 <- predict(mod1, newdata = df_base, type = "response")

## Probabilité de décès si tout le monde avait A0 = 0
y0 <- predict(mod0, newdata = df_base, type = "response")

head(data.frame(id = df_base$id, A0_obs = df_base$A0,
                pred_A0_1 = round(y1, 3),
                pred_A0_0 = round(y0, 3)))

Étape 3 — Calculer l’effet causal moyen (ATE)

Calculez \(\widehat{E}(Y^1) = \overline{y_1}\), \(\widehat{E}(Y^0) = \overline{y_0}\), puis l’ATE et le RR.

Montrez moi comment faire !
E_Y1 <- mean(y1)
E_Y0 <- mean(y0)

cat("E(Y^1) =", round(E_Y1, 3), "\n")
cat("E(Y^0) =", round(E_Y0, 3), "\n")
cat("ATE    =", round(E_Y1 - E_Y0, 3), "\n")
cat("RR     =", round(E_Y1 / E_Y0, 3), "\n")

On obtient : \(\widehat{E}(Y^1) \approx\) 0.125, \(\widehat{E}(Y^0) \approx\) 0.26, \(\widehat{ATE} \approx\) -0.135.

Interprétation

Que signifie cet ATE ? Comment le comparez-vous à l’analyse brute de la page précédente ?

L’ATE estimé par G-computation est la différence de risque de décès entre les scénarios contrefactuels “tout le monde exposé” et “personne exposé”, après ajustement sur \(X\) et \(L_0\).

Contrairement à l’analyse brute, cette estimation tient compte du fait que les individus exposés et non exposés ne sont pas comparables à l’inclusion.

Limites ici :

  • On traite le décès comme binaire, en ignorant l’information temporelle.
  • On n’ajuste que sur les confondeurs initiaux — les parties suivantes traitent les confondeurs dépendants du temps.
  • Le résultat repose entièrement sur la bonne spécification du modèle de résultat \(D\).

La Partie 2 répondra à la même Question 1 en utilisant une approche différente (IPTW), ce qui permettra un contrôle croisé des deux estimations.

Pour continuer, cliquez sur IPTW dans le menu à gauche.

BONUS — Intervalle de confiance par bootstrap

Estimez un intervalle de confiance à 95% par bootstrap (\(B = 200\) ré-échantillons).

Montrez moi comment faire !
set.seed(123)
B <- 200
ate_boot <- numeric(B)

for (b in 1:B) {
  idx  <- sample(nrow(df_base), replace = TRUE)
  df_b <- df_base[idx, ]

  m1 <- glm(D ~ X + L0, data = df_b[df_b$A0 == 1, ], family = binomial)
  m0 <- glm(D ~ X + L0, data = df_b[df_b$A0 == 0, ], family = binomial)

  y1_b <- predict(m1, newdata = df_b, type = "response")
  y0_b <- predict(m0, newdata = df_b, type = "response")

  ate_boot[b] <- mean(y1_b) - mean(y0_b)
}

cat("ATE    :", round(mean(ate_boot), 3), "\n")
cat("IC 95% :", round(quantile(ate_boot, c(0.025, 0.975)), 3), "\n")

BONUS 2 — G-computation sur critère censuré

Note

Cette section est à faire chez vous, à votre rythme. Elle n’est pas traitée en TP. Elle montre comment étendre la G-computation vue précédemment pour traiter correctement le critère de jugement censuré.

Pourquoi changer d’approche ?

La G-computation de la partie principale modélise \(D\) comme un indicateur binaire (décédé avant \(t = 3\) ans, oui/non). Deux limites importantes :

  1. Information temporelle perdue : on sait si l’individu est décédé, mais on ignore quand.
  2. Censure traitée comme un non-événement : un individu perdu de vue à \(t = 1\) an est traité de la même façon qu’un individu suivi 3 ans sans événement.

La G-computation complète contourne ces deux problèmes en estimant directement des courbes de survie contrefactuelles :

\[\hat{S}^{a_0 = 1}(t) = P(T^{a_0=1} > t) \qquad \text{et} \qquad \hat{S}^{a_0 = 0}(t) = P(T^{a_0=0} > t)\]

Principe

L’idée est simple : on remplace les deux modèles logistiques par deux modèles de Cox, un chez les exposés et un chez les non-exposés.

Dans un modèle de Cox, la survie individuelle est :

\[\hat{S}_i^{a_0}(t) = \exp\!\left(- \hat{H}_0^{a_0}(t) \times e^{\,\hat{\beta}^{a_0} \cdot (X_i,\, L_{0i})}\right)\]

  • \(\hat{H}_0^{a_0}(t)\) : hazard cumulatif de base (estimateur de Breslow), commun à tous les individus du groupe \(a_0\)
  • \(e^{\,\hat{\beta} \cdot (X_i, L_{0i})}\) : multiplicateur individuel de risque (prédicteur linéaire exponentié), propre à chaque individu

La courbe de survie marginale (contrefactuelle) s’obtient en moyennant sur l’ensemble de la population :

\[\hat{S}^{a_0}(t) = \frac{1}{n} \sum_{i=1}^n \hat{S}_i^{a_0}(t)\]

Étape 1 — Deux modèles de Cox

Note

df_base doit avoir été créée lors de l’Étape 0 ci-dessus. Si nécessaire, relancez ce chunk d’abord.

Estimez deux modèles de Cox séparément chez les exposés (\(A_0 = 1\)) et chez les non-exposés (\(A_0 = 0\)), en contrôlant sur \(X\) et \(L_0\). Utilisez coxph(Surv(time, D) ~ X + L0, data = ...).

Montrez moi comment faire !
mod_cox1 <- coxph(Surv(time, D) ~ X + L0,
                  data = df_base[df_base$A0 == 1, ])
mod_cox0 <- coxph(Surv(time, D) ~ X + L0,
                  data = df_base[df_base$A0 == 0, ])

summary(mod_cox1)
summary(mod_cox0)

Étape 2 — Hazard cumulatif de base et prédicteurs linéaires individuels

L’estimation des courbes contrefactuelles nécessite deux ingrédients issus de chaque modèle de Cox :

  • Le hazard cumulatif de base \(\hat{H}_0(t)\) : c’est une fonction en escalier qui augmente à chaque temps d’événement. On l’obtient avec basehaz(mod, centered = FALSE), qui renvoie un data frame avec les colonnes time et hazard.
  • Le prédicteur linéaire individuel \(\hat{\text{lp}}_i = \hat{\beta} \cdot (X_i, L_{0i})\) : c’est un scalaire par individu, qui mesure à quel point cet individu a un risque plus ou moins élevé que la baseline. On l’obtient avec predict(mod, newdata = df_base, type = "lp").

Extrayez le hazard cumulatif de base de chaque modèle, et calculez le prédicteur linéaire pour chaque individu de df_base sous chaque scénario.

Montrez moi comment faire !
## Hazard cumulatif de base (estimateur de Breslow)
bh1 <- basehaz(mod_cox1, centered = FALSE)
bh0 <- basehaz(mod_cox0, centered = FALSE)

head(bh1)  # time = temps d'événement, hazard = H0(t) cumulé

## Prédicteur linéaire individuel lp_i = beta * (X_i, L0_i)
## Prédit pour chaque individu de df_base, quel que soit son A0 observé
lp1 <- predict(mod_cox1, newdata = df_base, type = "lp")
lp0 <- predict(mod_cox0, newdata = df_base, type = "lp")

head(data.frame(
  id         = df_base$id,
  A0_obs     = df_base$A0,
  lp_si_A1   = round(lp1, 3),
  lp_si_A0   = round(lp0, 3)
))

Étape 3 — Courbes de survie individuelles et marginalisation

On dispose maintenant de tout ce qu’il faut pour calculer \(\hat{S}_i^{a_0}(t)\) pour chaque individu \(i\) et chaque temps \(t\), puis d’en faire la moyenne.

Astuce

stepfun(knots, values) crée une fonction en escalier à partir des nœuds et valeurs associées. Ici, stepfun(bh1$time, c(0, bh1$hazard)) crée la fonction \(\hat{H}_0^{(a_0=1)}(t)\) : elle vaut 0 avant le premier événement, puis prend la valeur de Breslow correspondante à chaque nœud. On peut ensuite l’évaluer sur n’importe quelle grille de temps.

outer(u, v) calcule le produit extérieur de deux vecteurs : outer(u, v)[i, j] = u[i] * v[j]. Cela permet de calculer \(-e^{\hat{\text{lp}}_i} \times \hat{H}_0(t_j)\) pour tous les couples \((i, j)\) d’un coup, sans boucle, et d’obtenir directement la matrice des \(n \times T\) valeurs \(\log \hat{S}_i(t_j)\).

Calculez les courbes de survie marginales \(\hat{S}^{a_0=1}(t)\) et \(\hat{S}^{a_0=0}(t)\) sur une grille de temps.

Montrez moi comment faire !
## 1. Grille de temps : tous les temps d'événement de df_base
t_grid <- sort(unique(df_base$time))

## 2. H0(t) évalué sur t_grid via des fonctions en escalier
H1_fun  <- stepfun(bh1$time, c(0, bh1$hazard))
H0_fun  <- stepfun(bh0$time, c(0, bh0$hazard))
H1_grid <- H1_fun(t_grid)   # vecteur longueur T
H0_grid <- H0_fun(t_grid)

## 3. Survie individuelle : S_i(t_j) = exp(-H0(t_j) * exp(lp_i))
##    outer(-exp(lp), H_grid)[i, j] = -exp(lp[i]) * H_grid[j]
S1_mat <- exp(outer(-exp(lp1), H1_grid))  # matrice n × T
S0_mat <- exp(outer(-exp(lp0), H0_grid))

## 4. Courbes marginales : moyenne sur les individus (ligne = individu)
S1_marg <- colMeans(S1_mat)
S0_marg <- colMeans(S0_mat)

head(data.frame(t     = round(t_grid, 2),
                S_a1  = round(S1_marg, 3),
                S_a0  = round(S0_marg, 3)))

Étape 4 — Visualisation des courbes de survie contrefactuelles

Tracez les deux courbes de survie marginales contrefactuelles sur le même graphique.

Montrez moi comment faire !
plot(t_grid, S1_marg, type = "s", col = "blue", lwd = 2,
     ylim = c(0, 1),
     xlab = "Temps (années)",
     ylab = "Probabilité de survie",
     main = "G-computation (Cox) — Courbes contrefactuelles")
lines(t_grid, S0_marg, type = "s", col = "red", lwd = 2)
legend("bottomleft",
       c("Scénario a₀=1 (tous exposés)",
         "Scénario a₀=0 (aucun exposé)"),
       col = c("blue", "red"), lty = 1, lwd = 2, bty = "n")

Étape 5 — Différence de survie à \(t = 3\) ans

Calculez \(\hat{S}^{a_0=1}(3)\), \(\hat{S}^{a_0=0}(3)\), et la différence de survie à 3 ans. Comparez avec l’ATE obtenu par G-computation logistique.

Montrez moi comment faire !
## Dernier temps <= 3 dans la grille
idx3 <- max(which(t_grid <= 3))

S1_3 <- S1_marg[idx3]
S0_3 <- S0_marg[idx3]

cat("S^(a0=1)(3) =", round(S1_3, 3), "\n")
cat("S^(a0=0)(3) =", round(S0_3, 3), "\n")
cat("Différence de survie à 3 ans :", round(S1_3 - S0_3, 3), "\n")
cat("Rapport de survie à 3 ans    :", round(S1_3 / S0_3, 3), "\n")

Les deux G-computations (logistique et Cox) ciblent en principe le même estimand : \(P(T^{a_0} \leq 3)\), ce qui correspond à \(1 - S^{a_0}(3)\). En pratique, la différence de survie à \(t = 3\) obtenue par les deux méthodes est souvent proche, mais peut légèrement différer car :

  • la version Cox utilise toute l’information temporelle (pas seulement le statut binaire à 3 ans) ;
  • la version Cox traite correctement la censure : un individu censuré à \(t = 1\) an ne contribue qu’à l’estimation jusqu’à son temps de censure, sans être compté comme “survivant” ou “décédé” à 3 ans ;
  • les deux modèles imposent des hypothèses fonctionnelles différentes (effets proportionnels des hazards pour Cox, effets log-linéaires sur la cote pour la logistique).

La version Cox est donc plus rigoureuse et constitue l’approche standard pour la G-computation avec un critère de survie.

Intervalle de confiance par bootstrap (modèle de Cox)

Estimez un intervalle de confiance à 95% par bootstrap (\(B = 200\) ré-échantillons) pour la différence de survie à \(t = 3\) ans.

Montrez moi comment faire !
set.seed(42)
B <- 200
S1_boot <- numeric(B)
S0_boot <- numeric(B)

for (b in 1:B) {
  ## 1. Ré-échantillonnage avec remise
  idx <- sample(nrow(df_base), replace = TRUE)
  db  <- df_base[idx, ]

  ## 2. Modèles de Cox
  m1 <- coxph(Surv(time, D) ~ X + L0, data = db[db$A0 == 1, ])
  m0 <- coxph(Surv(time, D) ~ X + L0, data = db[db$A0 == 0, ])

  ## 3. Basehaz et prédicteurs linéaires
  bh1b <- basehaz(m1, centered = FALSE)
  bh0b <- basehaz(m0, centered = FALSE)
  lp1b <- predict(m1, newdata = db, type = "lp")
  lp0b <- predict(m0, newdata = db, type = "lp")

  ## 4. H0 au temps 3 (dernière valeur <= 3 dans la baseline hazard)
  H1_3b <- if (any(bh1b$time <= 3)) tail(bh1b$hazard[bh1b$time <= 3], 1) else 0
  H0_3b <- if (any(bh0b$time <= 3)) tail(bh0b$hazard[bh0b$time <= 3], 1) else 0

  S1_boot[b] <- mean(exp(-H1_3b * exp(lp1b)))
  S0_boot[b] <- mean(exp(-H0_3b * exp(lp0b)))
}

diff_boot <- S1_boot - S0_boot
cat("Différence de survie à 3 ans (G-computation Cox, bootstrap) :\n")
cat("  Estimée :", round(mean(diff_boot), 3), "\n")
cat("  IC 95%  :", round(quantile(diff_boot, c(0.025, 0.975)), 3), "\n")