14 Annexes

14.1 Tests Grambsch-Therneau OLS sur les résidus de Schoenfeld

Important

Attention il ne s’agit pas du test actuellement implémenté dans la nouvelle version de survival (v3) qui, malheureusement, lui a substitué la version dite exacte (moindres carrés généralisés). Le programme de la fonction du test OLS est néanmoins facilement récupérable et exécutable. lien.

Je continue de préconiser l’utilisation de cette version OLS du test, reproductible avec les autres applications statistiques (Stata,Sas,Python).

Le test dit “simplifié”, qui n’apparait pas dans le texte original de P.Gramsch et T.Thernau lien, répond à un soucis d’instabilité des variances des résidus de Schoenfeld en fin de durée d’observation lorsque peu d’observation restent soumises au risque. Cet argument est soulevé dans leur ouvrage de 2022 lien avant d’en présenter sa version.
Il est simplifié car on applique à tous les résidus bruts la variance du paramètre ($b$) estimés par le modèle de Cox.
Le test devient alors un simple test de corrélation entre les résidus et une fonction de la durée (centrée). Dans l’esprit, il peut être également approché par une regression linéaire par les moindre carrés ordinaires entre les résidus et une fonction de la durée (voir page 134 de l’ouvrage de Grambsch et Therneau).

Soit les données suivantes, avec t la variable de durées, Y la variable de censure et X la seule et unique covariable.

Pas d’évènement simultané (donc pas de correction de la vraisemblance)
Covariable de type indicatrice

$t_i$	$Y_i$	$X_i$
1	1	1
2	0	0
3	0	0
4	1	1
5	1	1
6	1	0
7	0	1

test = data.frame(time=  c(1,2,3,4,5,6,7),
                    Y=c(1,0,0,1,1,1,0),
                    X=     c(1,0,0,1,1,0,1))

Estimation du modèle de Cox:

library(survival)
fit = coxph(formula = Surv(time, Y) ~ X, data=test)
fit

Call:
coxph(formula = Surv(time, Y) ~ X, data = test)

    coef exp(coef) se(coef)    z     p
X 0.6217    1.8622   1.1723 0.53 0.596

Likelihood ratio test=0.31  on 1 df, p=0.5797
n= 7, number of events= 4

Calcul des résidus brut (si et seulement si $Y=1$) dans le cas d’une seule covariable avec $b$ égal à 0.62:

\[rs_{i}=X_{i}- \sum_{j\in R_i}X_{i}\frac{e^{0.62\times X}}{\sum_{j\in R_i}e^{0.62\times X}}= X_{i} - E(X_{j\in R_i})\] Il y a ici 4 résidus à calculer, pour $t=(1,4,5,6)$

Résidus pour $t=1$

$a_1= \sum_{j\in R_i}e^{0.62\times X} = e^{0.62} + 1 + 1 + e^{0.62} + 1 + e^{0.62}= 10.43$
$b_1= \sum_{j\in R_i}X_{i}\frac{e^{0.62\times X}}{\sum_{j\in R_i}e^{0.62\times X}} = 4\times\frac{e^{0.62}}{10.43} = 0.71$
$r_1 = 1 - 0.71 = 0.29$

Résidus pour $t=4$

$a_4 = e^{0.62} + e^{0.62} + 1 + e^{0.62} = 6.58$
$b_4 = 4\times\frac{e^{0.62}}{6.58} = 0.84$
$r_4 = 1 - 0.84 = 0.15$

Résidus pour $t=5$

$a_5 = e^{0.62} + e^{0.62} + 1 = 4.71$
$b_5 = 2\times\frac{e^{0.62}}{4.71} = 0.78$
$r_5 = 1 - 0.78 = 0.21$

Résidus pour $t=6$

$a_6 = e^{0.62} + 1 = 2.86$
$b_6 = \frac{e^{0.62}}{2.86} = 0.65$
$r_6 = 0 - 0.65 = -0.65$

Les résidus “standardisés”, ou plutôt scaled residuals (je cale sur une traduction correcte en français) sont égaux à:

\[sr_i = b + nd \times Var(b) \times r_i\] Avec $nd= \sum Y_i$

$\sum Y_i = 4$
$Var(b) = (1.1723)^2=1.37$
$sr_1 = 0.62 + 4\times 1.37 \times 0.29 = 2.20$
$sr_4 = 0.62 + 4\times 1.37 \times 0.15 = 1.47$
$sr_5 = 0.62 + 4\times 1.37 \times 0.21 = 1.78$
$sr_6 = 0.62 + 4\times 1.37 \times (-0.65) = -2.95$

Avec $g(t_i)$ une fonction de la durée ($g(t_i)=t_i$, $g(t_i)=1-KM(t_i)$…) et $\overline{g(t)}$ sa valeur moyenne, la statistique du test score simplifié pour une covariable est égale à :

\[\frac{[\sum_i(g(t_i) - \overline{g(t_i)}\times sr_i)]^2}{nd \times Var(b) \times (\sum_i(g(t_i) - \overline{g(t_i)})^2}\] Et suis un $\chi^2$ à 1 degré de liberté.

Avec $\overline{g(t_i)}=t_i$, le calcul de la statistique de test est:

$\overline{g(t_i)}= \frac{28}{7}=4$
$\frac{[(1-4)\times 2.20] + [(4-4)\times 1.47 + (5-4)\times 1.78 + (6-4)\times (-2.95)]^2 }{4\times 1.37 \times [(1-4)^2 + (4-4)^2 + (5-4)^2 + (6-4)^2] } = \frac{114.9}{76.72} = 1.49$

#source("D:/D/Marc/SMS/FORMATIONS/analyse_duree/cox.zphold/cox.zphold.R")

source("https://raw.githubusercontent.com/mthevenin/analyse_duree/master/cox.zphold/cox.zphold.R")

cox.zphold(fit, transform="identity")

Warning in is.R(): 'is.R' est obsolète.
Voir help("Deprecated") et help("base-deprecated").

     rho chisq     p
X -0.688  1.49 0.222

14.2 Fragilité et immunité

Seulement quelques remarques, le traitement de ces problématiques dépassant largement le contenu de la formation.

14.2.1 Fragilité (Frailty)

Pour la fragilité, je conseille fortement de lire la dernière section du document de travail de Simon Quantin (cf bibliographie), il n’y a pas meilleure présentation du problème que la sienne ^[petite maj par rapport à la version précédente: il ne traite que la fragilité individuelle stricto sensu et non la fragilité plus connu sous le terme de shared frailty* (proche modèle multiniveau). Problèmatique importante, car une des origines de la non proportionnalité des risques réside dans l’omission de variables. Ici on va être confronté une omission sur des traits non observables ou latents, qui accélèrent dès le début de la période d’exposition la survenue de l’évènement. L’introduction d’un facteur de fragilité se fait par l’introduction d’un effet aléatoire dans le modèle, de nature plus complexe, et rendant l’interprétation des modèles plus compliquée.

On peut distinguer deux types de modèles:

les modèles à fragilité partagée, c’est la situation la plus simple car la logique se rapproche des modèles multiniveaux, des groupes d’individus, identifiables, partagent une même fragilité, par exemple géographique.
les modèles à fragilité non partagée, avec des caractéristiques latentes non observable comme les préférences, ou en médecine certains traits génétiques non identifiés.

14.2.2 Immunité (Cure fraction)

Le phénomène d’immunité est un cas particulier du précédent, et a été étudié dès le début des années 1950, en questionnant l’exposition au risque d’une partie des observations. On s’interrogeait par exemple sur les risques de rechute et de décès après le traitement d’un premier cancer. Visuellement on peut commencer à se proser des questions sur la présence d’une fraction immunisée ou non susceptible de connaître l’évènement lorsque la fonction de séjour ne tend pas vers 0 mais présente une longue asymptote (plateau) sur une valeur supérieure à 0: $\lim_{t \to \infty}S(t)=a$.

Les modèles avec une fraction immunisée peuvent être de type mixte en associant une probabilité d’être immunisé aux observations censurées à droite à un modèle de durée ¹. Plus dans le vent je crois, on a également des modèles de type non mixte, avec il me semble une connotation bayesienne qui semble s’accroître. Il n’y a donc pas de méthode unifiée à ce jour [Si vous voulez vous en convaincre].

On peut également noter, c’est important, que cette problématique affecte les analyses avec des évènements dits récurrents. Ici, la stratégie classique qui consiste à introduire dans un modèle un simple effet aléatoire de type fragilité partagée (shared frailty) pour contrôler risque d’être insuffisante. Ici le groupe est constitué de chaque séquence de remise dans le risque set. Exemple pour la fécondité: une personne ayant eu un enfant est exposée au risque d’en avoir un autre, l’horloge temporelle étant alors simplement réinitialisée. Et donc, quid des préférences individuelles en terme de fécondité ².

Enfin, les modèles à fragilité ou à fraction immunisée repose tous sur une hypothèse très forte. La fragilité ou le degré d’immunité est toujours défini (estimé) en début d’exposition, et il ne varie pas. Cela peut ne pas toujours faire sens, en particulier pour les préférences, pas forcément stables ou fixes dans le temps.

Le plus classique utilise un algorithme Expectation Maximisation utilisé en imputation: on estime une probabilité d’être susceptible de connaitre l’évènement aux observations censurées à droite, qui intervient comme facteur de pondération dans le modèle de durée. Cette probabilité et le modèle de durée qui lui est associé est réévalué à chaque boucle de l’algorithme jusqu’à convergence. Le principale problème de cette méthode résite dans l’estimation de la variance, souvent effectué par bootstrap. Cette méthode à l’avantage d’être implémentable en durée discrète, bien qu’à ma connaissance aucun logiciel ne la propose (j’ai une commande Stata encore perfectible sous le coude). On trouve en revanche ce type d’estimation sous R, pour les modèles de Cox ou les modèles paramétriques dans le package smcure↩︎
en situation de récurrence, toujours penser à remettre à jour les conditions initiales, par exemple pour la fécondite l’âge de la mère à la naissance de l’enfant pour les rang supérieur à 1↩︎

--- # title: "Annexe" --- # Annexes ## Tests Grambsch-Therneau OLS sur les résidus de Schoenfeld ::: callout-important Attention il ne s'agit pas du test actuellement implémenté dans la nouvelle version de `survival` (v3) qui, malheureusement, lui a substitué la version dite *exacte* (moindres carrés généralisés). Le programme de la fonction du test OLS est néanmoins facilement récupérable et exécutable. [lien](https://github.com/mthevenin/analyse_duree/blob/main/cox.zphold/cox.zphold.R). Je continue de préconiser l'utilisation de cette version OLS du test, reproductible avec les autres applications statistiques (Stata,Sas,Python). ::: * Le test dit "simplifié", qui n'apparait pas dans le texte original de P.Gramsch et T.Thernau [lien](https://www.jstor.org/stable/2337123#metadata_info_tab_contents), répond à un soucis d'instabilité des variances des résidus de Schoenfeld en fin de durée d'observation lorsque peu d'observation restent soumises au risque. Cet argument est soulevé dans leur ouvrage de 2022 [lien](https://link.springer.com/book/10.1007/978-1-4757-3294-8) avant d'en présenter sa version. * Il est simplifié car on applique à tous les résidus bruts la variance du paramètre ($b$) estimés par le modèle de Cox. * Le test devient alors un simple test de corrélation entre les résidus et une fonction de la durée (centrée). Dans l'esprit, il peut être également approché par une regression linéaire par les moindre carrés ordinaires entre les résidus et une fonction de la durée (voir page 134 de l'ouvrage de Grambsch et Therneau). Soit les données suivantes, avec *t* la variable de durées, *Y* la variable de censure et *X* la seule et unique covariable. * Pas d'évènement simultané (donc pas de correction de la vraisemblance) * Covariable de type indicatrice | $t_i$ | $Y_i$ | $X_i$ | |-----|-----|-----| | 1 | 1 | 1 | | 2 | 0 | 0 | | 3 | 0 | 0 | | 4 | 1 | 1 | | 5 | 1 | 1 | | 6 | 1 | 0 | | 7 | 0 | 1 | ```{r} test = data.frame(time= c(1,2,3,4,5,6,7), Y=c(1,0,0,1,1,1,0), X= c(1,0,0,1,1,0,1)) ``` Estimation du modèle de Cox: ```{r} library(survival) fit = coxph(formula = Surv(time, Y) ~ X, data=test) fit ``` Calcul des résidus brut (si et seulement si $Y=1$) dans le cas d'une seule covariable avec $b$ égal à **0.62**: $$rs_{i}=X_{i}- \sum_{j\in R_i}X_{i}\frac{e^{0.62\times X}}{\sum_{j\in R_i}e^{0.62\times X}}= X_{i} - E(X_{j\in R_i})$$ Il y a ici 4 résidus à calculer, pour $t=(1,4,5,6)$ **Résidus pour $t=1$** * $a_1= \sum_{j\in R_i}e^{0.62\times X} = e^{0.62} + 1 + 1 + e^{0.62} + 1 + e^{0.62}= 10.43$ * $b_1= \sum_{j\in R_i}X_{i}\frac{e^{0.62\times X}}{\sum_{j\in R_i}e^{0.62\times X}} = 4\times\frac{e^{0.62}}{10.43} = 0.71$ * $r_1 = 1 - 0.71 = 0.29$ **Résidus pour $t=4$** * $a_4 = e^{0.62} + e^{0.62} + 1 + e^{0.62} = 6.58$ * $b_4 = 4\times\frac{e^{0.62}}{6.58} = 0.84$ * $r_4 = 1 - 0.84 = 0.15$ **Résidus pour $t=5$** * $a_5 = e^{0.62} + e^{0.62} + 1 = 4.71$ * $b_5 = 2\times\frac{e^{0.62}}{4.71} = 0.78$ * $r_5 = 1 - 0.78 = 0.21$ **Résidus pour $t=6$** * $a_6 = e^{0.62} + 1 = 2.86$ * $b_6 = \frac{e^{0.62}}{2.86} = 0.65$ * $r_6 = 0 - 0.65 = -0.65$ Les résidus "standardisés", ou plutôt *scaled residuals* (je cale sur une traduction correcte en français) sont égaux à: $$sr_i = b + nd \times Var(b) \times r_i$$ Avec $nd= \sum Y_i$ * $\sum Y_i = 4$ * $Var(b) = (1.1723)^2=1.37$ * $sr_1 = 0.62 + 4\times 1.37 \times 0.29 = 2.20$ * $sr_4 = 0.62 + 4\times 1.37 \times 0.15 = 1.47$ * $sr_5 = 0.62 + 4\times 1.37 \times 0.21 = 1.78$ * $sr_6 = 0.62 + 4\times 1.37 \times (-0.65) = -2.95$ Avec $g(t_i)$ une fonction de la durée ($g(t_i)=t_i$, $g(t_i)=1-KM(t_i)$...) et $\overline{g(t)}$ sa valeur moyenne, la statistique du test score simplifié pour une covariable est égale à : $$\frac{[\sum_i(g(t_i) - \overline{g(t_i)}\times sr_i)]^2}{nd \times Var(b) \times (\sum_i(g(t_i) - \overline{g(t_i)})^2}$$ Et suis un $\chi^2$ à 1 degré de liberté. Avec $\overline{g(t_i)}=t_i$, le calcul de la statistique de test est: * $\overline{g(t_i)}= \frac{28}{7}=4$ * $\frac{[(1-4)\times 2.20] + [(4-4)\times 1.47 + (5-4)\times 1.78 + (6-4)\times (-2.95)]^2 }{4\times 1.37 \times [(1-4)^2 + (4-4)^2 + (5-4)^2 + (6-4)^2] } = \frac{114.9}{76.72} = 1.49$ ```{r} #source("D:/D/Marc/SMS/FORMATIONS/analyse_duree/cox.zphold/cox.zphold.R") source("https://raw.githubusercontent.com/mthevenin/analyse_duree/master/cox.zphold/cox.zphold.R") ``` ```{r} cox.zphold(fit, transform="identity") ``` ## Fragilité et immunité Seulement quelques remarques, le traitement de ces problématiques dépassant largement le contenu de la formation. ### Fragilité (Frailty) Pour la *fragilité*, je conseille fortement de lire la dernière section du document de travail de ***Simon Quantin** (cf bibliographie), il n'y a pas meilleure présentation du problème que la sienne ^[petite maj par rapport à la version précédente: il ne traite que la fragilité individuelle stricto sensu et non la fragilité plus connu sous le terme de *shared frailty* (proche modèle multiniveau). Problèmatique importante, car une des origines de la non proportionnalité des risques réside dans l’omission de variables. Ici on va être confronté une omission sur des traits non observables ou latents, qui **accélèrent** dès le début de la période d’exposition la survenue de l’évènement. L’introduction d’un facteur de fragilité se fait par l’introduction d’un effet aléatoire dans le modèle, de nature plus complexe, et rendant l’interprétation des modèles plus compliquée. On peut distinguer deux types de modèles: - les modèles à *fragilité partagée*, c'est la situation la plus simple car la logique se rapproche des modèles multiniveaux, des groupes d'individus, identifiables, partagent une même *fragilité*, par exemple géographique. - les modèles à *fragilité non partagée*, avec des caractéristiques latentes non observable comme les préférences, ou en médecine certains traits génétiques non identifiés. ### Immunité (Cure fraction) Le phénomène d’immunité est un cas particulier du précédent, et a été étudié dès le début des années 1950, en questionnant l’exposition au risque d’une partie des observations. On s'interrogeait par exemple sur les risques de rechute et de décès après le traitement d'un premier cancer. Visuellement on peut commencer à se proser des questions sur la présence d'une fraction *immunisée* ou non *susceptible* de connaître l'évènement lorsque la fonction de séjour ne tend pas vers 0 mais présente une longue asymptote (plateau) sur une valeur supérieure à 0: $\lim_{t \to \infty}S(t)=a$. Les modèles avec une fraction immunisée peuvent être de ***type mixte*** en associant une probabilité d'être immunisé aux observations censurées à droite à un modèle de durée ^[Le plus classique utilise un algorithme *Expectation Maximisation* utilisé en imputation: on estime une probabilité d'être susceptible de connaitre l'évènement aux observations censurées à droite, qui intervient comme facteur de pondération dans le modèle de durée. Cette probabilité et le modèle de durée qui lui est associé est réévalué à chaque boucle de l'algorithme jusqu'à convergence. Le principale problème de cette méthode résite dans l'estimation de la variance, souvent effectué par bootstrap. Cette méthode à l'avantage d'être implémentable en durée discrète, bien qu'à ma connaissance aucun logiciel ne la propose (j'ai une commande Stata encore perfectible sous le coude). On trouve en revanche ce type d'estimation sous R, pour les modèles de Cox ou les modèles paramétriques dans le package **`smcure`**]. Plus dans le vent je crois, on a également des modèles de **type non mixte**, avec il me semble une connotation bayesienne qui semble s'accroître. Il n’y a donc pas de méthode unifiée à ce jour [[Si vous voulez vous en convaincre](https://www.annualreviews.org/doi/10.1146/annurev-statistics-031017-100101)]. On peut également noter, c'est important, que cette problématique affecte les analyses avec des évènements dits récurrents. Ici, la stratégie classique qui consiste à introduire dans un modèle un simple effet aléatoire de type fragilité partagée (shared frailty) pour contrôler risque d'être insuffisante. Ici le *groupe* est constitué de chaque séquence de remise dans le risque set. Exemple pour la fécondité: une personne ayant eu un enfant est exposée au risque d'en avoir un autre, l'horloge temporelle étant alors simplement réinitialisée. Et donc, quid des préférences individuelles en terme de fécondité ^[en situation de récurrence, toujours penser à *remettre à jour* les conditions initiales, par exemple pour la fécondite l'âge de la mère à la naissance de l'enfant pour les rang supérieur à 1]. Enfin, les modèles à fragilité ou à fraction immunisée repose tous sur une hypothèse très forte. La fragilité ou le degré d'immunité est toujours défini (estimé) en début d'exposition, et il ne varie pas. Cela peut ne pas toujours faire sens, en particulier pour les préférences, pas forcément stables ou fixes dans le temps.

\(t_i\)	\(Y_i\)	\(X_i\)
1	1	1
2	0	0
3	0	0
4	1	1
5	1	1
6	1	0
7	0	1