The French version of this document is no longer maintained: be sure to check the more up-to-date English version.

Analyse de la survie

Introduction
Quelques distributions
Estimation de Kaplan-Meier
Modèles à temps discret
Modèle de régression paramétrique
Proportional Hazard (PH) models
Accelerated Failure Time (AFT) models
Autres modèles
Cox Proportional Hazard model
A FAIRE

A FAIRE : terminer d'écrire cette partie (il me reste encore à comprendre la régression de Cox...)

A FAIRE : plan de cette partie

Analyse de la survie sans variables prédictives
  Fonction de survie, de risque, de risque cumulé
  Quelques distributions
  Méthodes non paramétriques : estimation de Kaplan-Meier
  Méthodes paramétriques : A FAIRE    
Modèles discrets
Analyse de la survie avec variables prédictives
  A FAIRE

A FAIRE :

Il convient de distinguer les données de survie discrète des données
de survie continues.
Très souvent, on ne parle que des données continues et on traite les
données discrètes comme si elles étaient continues -- ça n'est pas
une très bonne idée.

A FAIRE :

Les données de survie son assez difficiles à manipuler : on les
remplacera par la fonction de survie, ou le taux de survie cumulé ou
le taux de survie. Ce sont des fonctions, c'est plus facile.

A FAIRE :

Insuffisance de la régression classique : elle ne permet pas de
traiter la censure des données, ni de prendre en compte des
variables prédictives qui dépendent du temps.

A FAIRE :

Pas de variables prédictives :
- méthodes paramétriques (MLE, avec les modèles exposés ci-après)
- méthodes non paramétriques (Kaplan-Meier (cas continu) ou
  LifeTable (cas discret)

Avec des variables prédictives : 
- méthodes paramétriques
- méthodes semi-paramétriques (régression de Cox : paramétrique pour
  la partie qui dépend des variables prédictives, non paramétriques
  pour le reste).

A FAIRE : frailty = unobserved heterogeneity.

A FAIRE : quelle différence entre données de survie et données longitudinales ?

Introduction

Définitions et exemples

On s'intéresse ici à des variables de survie : une variable de survie, c'est à peu près comme une variable quantitative, mais elle peut prendre des valeurs numériques précises, comme "13 ans", mais aussi des valeurs moins précises, du genre "plus de 15 ans".

On qualifie de "valeurs censurées" ces valeurs moins précises. On les note généralement ainsi :

10+    Plus de 10 (donnée censurée à droite)
5-     Moins de 5 (censurée à gauche)
[5,10] Entre 5 et 10

On les rencontre par exemple quand on regarde la durée de survie après un cancer : si l'étude dure 10 ans, on pourra donner la durée de survie de tous les patients qui sont morts, mais pour les autres, on pourra juste dire que leur durée de survie est supérieure à 10 ans. La censure à gauche correspond à des évènements qui ont eu lieu avant le début de l'étude (mais tous les individus ne rentrent pas dans l'étude en même temps). A FAIRE : donner un exemple concret de censure à gauche. La censure par intervalle se rencontre par exemple quand on constate l'apparition d'un symptome lors d'une visite médicale : il est apparu depuis la visite précédente, mais on ne sait pas exactement quand.

Voici d'autres exemples : durée de survie d'un patient, durée de survie d'un mariage, durée de survie d'une entreprise, durée de fonctionnement sans panne d'un produit manufacturé, durée pendant laquelle un produit reste suffisemment populaire pour continuer à être produit, durée de survie d'un enseignant en lycée (jusqu'à sa démission, son internement en hopital psychiatrique, son suicide), durée jusqu'à l'arret de prise d'un médicament, etc.

On peut aussi voir les variables de survie comme des variables qualitatives binaires (celles qu'on cherche à prédire dans la régression logistique) indiquant la survenue ou non d'un évènement, auxquelles on a ajouté l'instant de survenue de cet évènement. Les valeurs censurées correspondent à << l'évènement ne s'est pas produit >> : dans les exemples précédents, cela revient à remplacer les valeurs << Vivant / Mort >> par << Vivant / Date de la mort >>.

Les variables de survie peuvent avoir plusieurs valeurs : par exemple sain, malade et mort (si un patient guérit entre deux maladies, ça devient encore plus compliqué). Dans ces situations plus complexes, on parle parfois de données transitionnelles (car on s'intéresse à la transition d'un état vers un autre) ou de "spell duration data".

On peut aussi voir une donnée de survie comme une série temporelle qualitative, dans laquelle on s'intéresse aux transitions d'un état à l'autre.

A FAIRE :

Competing risk models
Un état de départ et plusieurs états absorbants
(exemple : passage du chomage à un emploi ou du chomage à une sortie
du marché de l'emploi ; autre exemple : concurrence entre mort du
cancer et mort de maladie cardiovasculaire).

Survival function, Hazard function

Si T est une variable de survie, sa fonction de survie est

S(t) = P(T>t) = 1-F(t).

Par exemple, si T désigne le nombre d'années de survie après la diagnostic d'un cancer, S(t) est la probabilité de survivre au moins t années.

Sa "hazard function", ou "hazard rate" (je ne sais pas comment on dit en français : fonction de risque ?) est

                    P( t < T <= t + u | T > t )
lambda(t) =  lim   -----------------------------
            u -> 0              u

             dF/dt
          = -------
               S

               d ln S
          = - --------
                 dt

La formule ressemble beaucoup à la définition de la densité, mais on a ajouté la condition "T>t".

A la différenciation près, c'est très semblable à la différence entre l'espérance de vie à la naissance et l'espérance de vie à 60 ans, ou la différence entre la probabilité de mourir entre 12 et 13 ans et la probabilité de mourir entre 12 et 13 ans sachant que le patient a déjà survécu 12 ans.

On rencontre aussi parfois la "cumulative hazard function", primitive de la précédente, qui représente le risque cumulé.

Lambda(t) = - ln S(t)

Quelques distributions

Distribution exponentielle

Elle correspond à une "hazard function" constante,

lambda(t) = lambda
Lambda(t) = lambda t
S(t) = exp( - lambda t )

op <- par(mfrow=c(3,1))
n <- 20
lambda <- rep(2,n)
x <- seq(0,2,length=n)
plot(lambda ~ x, type='l', main=expression(lambda))
plot(lambda*x ~ x, type='l', main=expression(Lambda))
plot(exp(-lambda*x) ~ x, type='l', main="S")
par(op)

*

Distribution de Weibull

Elle correspond à une "hazard function" de la forme

Lambda(t) = a * t ^ gamma
S(t) = exp( - a * t ^ gamma )

op <- par(mfrow=c(3,1))
n <- 20
alpha <- 1
g <- rep(2,n)
x <- seq(0,2,length=n)
plot(g * alpha * x^(g-1) ~ x, type='l', main=expression(lambda (gamma==2)))
plot(alpha * x^g ~ x, type='l', main=expression(Lambda))
plot(exp(-alpha*x^g) ~ x, type='l', main="S")
par(op)

*

op <- par(mfrow=c(3,1))
n <- 20
alpha <- 1
g <- rep(.5,n)
x <- seq(0,2,length=n)
plot(g * alpha * x^(g-1) ~ x, type='l', main=expression(lambda (gamma==.5)))
plot(alpha * x^g ~ x, type='l', main=expression(Lambda))
plot(exp(-alpha*x^g) ~ x, type='l', main="S")
par(op)

*

C'est par exemple ce qu'on utilise pour modéliser la durée de fonctionnement sans panne d'une machine : plus la machine est vieille, plus le risque de panne est grand ; le risque de panne dépend de l'âge de la machine.

Distribution Gamma

Considérons une machine très fiable : elle possède des circuits de secours. Pour qu'elle tombre réellement en panne, il faut que plusieurs circuits tombent en panne. Le temps d'attente avant une panne de cette machine est le temps d'attente de trois pannes de ces circuits.

A FAIRE : distribution Gamma

Distribution log-normale

A FAIRE : expliquer

op <- par(mfrow=c(3,1))
n <- 200
x <- seq(0,5,length=n)
plot(1-plnorm(x) ~ x, type='l', main="S")
L <- function (x) { -log(1-plnorm(x)) }
plot(L(x) ~ x, type='l', main=expression(Lambda))
h <- .001
plot( (L(x+h)-L(x))/h ~ x, type='l', main=expression(lambda))
par(op)

*

Autres modèles

Weibul: lambda = alpha * t^(alpha-1) * exp(+X béta)
???
A FAIRE : il y a une variable prédictive, qui n'a pas encore sa
place ici...

log-logistic model:

                  1
S(t) = --------------------------  with lambda = exp(- X béta)
        1 + (lambda t)^(1/gamma)


Gompertz: lambda = exp( exp(+X béta) + gamma t )
(important)

Generalized Gamma
(généralise les modèles précédents: kappa=1:Weibul,
kappa=0:lognormal, kappa=sigma:gamma)

Il y a aussi des modèles à temps discret :

logistic
complementary log-log

Estimation de Kaplan-Meier

Estimation de la fonction de survie selon Kaplan--Meier (product-limit Kaplan--Meier estimator)

S'il n'y a pas de censure, on peut estimer la fonction de survie ainsi :

        Card( i : T_i>t )
S(t) = -------------------
              n

En présence de censure, on tient compte des données censurées jusqu'à ce qu'on arrive à leur valeur, et après on les oublie. Plus précisément, si le temps est discret :

S(t0) = P(vivant à t=t0 | vivant à t=t0-1) * P(vivant à t=t0-1)

Par exemple, la fonction de survie de la série

1 2 2 2+ 3+ 3+ 3+ 4 4 4 4 ... (100 individus)

se calcule ainsi :

temps sujets morts censures       p1     p2    S=p1*p2
------------------------------------------------------
1        100     1        0   99/100      1        .99 
2         99     2        1   97/100    .99      .9603
3         96     0        3   96/100  .9603      .9219
4         93     3        0   90/100  .9219      .8297

Les valeurs de S sont des produits de plus en plus gros :

S(0) = 1
S(1) = 1 * 99/100
S(2) = 1 * 99/100 * 97/100
S(3) = 1 * 99/100 * 97/100 * 96/100
S(4) = 1 * 99/100 * 97/100 * 96/100 * 90/100
etc.

On peut donner une formule générale :

                              di
S(t) =    Produit      ( 1 - ---- )
       i tel que ti<t         ni

On peut calculer des intervalles de confiance en supposant que Lambda (qui vaut -log(S)) est normale.

A FAIRE : comment a-t-on la variance de Lambda ???

Il existe d'autres estimateurs de la fonction de survie, par exemple celui d'Altschuler--Nelson--Aalen--Flemming--Harrington :

                             di
Lambda(t) =      Somme     ------.
            i tel que ti<t   ni

Ces estimateurs ont une variance élevée : si on a un modèle raisonnable, il vaut mieux préférer un estimateur paramétrique obtenu par le maximum de vraissemblance (on peut signaler que l'estimateur de Kaplan--Meier correspond à un maximum de vraissemblance (non paramétrique)).

Analyse de survie sous R

La plupart des fonctions sont dans le paquetage "survival".

La commande Surv permet de créer des variables de survie.

> Surv(c(1,2,2,2,3,3,3,4,4,4,4),
+      c(1,1,1,0,0,0,0,1,1,1,1))
[1] 1  2  2  2+ 3+ 3+ 3+ 4  4  4  4

La commande survfit permet de calculer la fonction de survie.

> x <- Surv(c(1,2,2,2,3,3,3,4,4,4,4),
+           c(1,1,1,0,0,0,0,1,1,1,1))

> survfit(x)
Call: survfit(formula = x)
        n    events     rmean se(rmean)    median   0.95LCL   0.95UCL
   11.000     7.000     3.364     0.322     4.000       Inf       Inf

> summary(survfit(x))
Call: survfit(formula = x)
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    1     11       1    0.909  0.0867        0.754            1
    2     10       2    0.727  0.1343        0.506            1
    4      4       4    0.000      NA           NA           NA

ou de la tracer

set.seed(87638)
library(survival)
n <- 200
x <- rweibull(n,.5)
y <- rexp(n,1/mean(x))
s <- Surv(ifelse(x<y,x,y), x<=y)
plot(s) # peu informatif

*

plot(survfit(s))
lines(survfit(s, type='fleming-harrington'), col='red')
r <- survfit(s)
lines( 1-pweibull( r$time, .5 ) ~ r$time, lty=3, lwd=3, col='blue' )
legend( par("usr")[2], par("usr")[4], yjust=1, xjust=1,
        c("Kaplan-Meier", "Fleming-Harrington", "Fonction de survie théorique"),
        lwd=c(1,1,3), lty=c(1,1,3),
        col=c(par("fg"), 'red', 'blue'))
title(main="Fonction de survie (estimation de Kaplan-Meier)")

*

op <- par(mfrow=c(3,1))
r <- survfit(s)
plot(r$surv ~ r$time, type='l', main="S")
curve( 1-pweibull(x,.5,1), col='red', lty=2, add=T )
plot(-log(r$surv) ~ r$time, type='l', main=expression(Lambda))
curve( -log(1-pweibull(x,.5,1)), col='red', lty=2, add=T )
# Avant de dériver, on commence par lisser Lambda
a <- -log(r$surv)
b <- r$time
library(modreg)
l <- loess(a~b)
bb <- seq(min(b),max(b),length=200)
aa <- predict(l, data.frame(b=bb))
plot( diff(aa) ~ bb[-1], type='l', main=expression(lambda) )
aa <- -log(1-pweibull(bb,.5,1))
lines( diff(aa) ~ bb[-1], col='red', lty=2 )
par(op)

*

Remarque : on peut obtenir certains des graphiques ci-dessus en modifiant l'option "fun".

plot(r, fun="log", main="log-survival curve")

*

plot(r, fun="event", main="cumulative events: f(y)=1-y")

*

plot(r, fun="cumhaz", main="cumulative hazard: f(y)=-log(y)=Lambda")

*

try(
  plot(r, fun="cloglog", main="complementary log-log plot")
)
# f(y)=log(-log(y)), log-scale on the x-axis

*

On peut aussi avoir un (ou plusieurs) facteurs -- nous en reparlerons plus loin.

n <- 200
x1 <- rweibull(n,.5)
x2 <- rweibull(n,1.2)
f <- factor( sample(1:2, n, replace=T), levels=1:2 )
x <- ifelse(f==1,x1,x2)
y <- rexp(n,1/mean(x))
s <- Surv(ifelse(x<y,x,y), x<=y)
plot(s, col=as.numeric(f))

*

plot( density(s[,1][ s[,2] == 1 & f == 1]), lwd=3,
      main="Analyse de la survie, un facteur" )
lines( density(s[,1][ s[,2] == 0 & f == 1]), lty=2, lwd=3 )
lines( density(s[,1][ s[,2] == 1 & f == 2]), col='red', lwd=3 )
lines( density(s[,1][ s[,2] == 0 & f == 2]), lty=2, col='red', lwd=3 )
legend( par("usr")[2], par("usr")[4], yjust=1, xjust=1,
        c("non censored, f = 1", "censored, f = 1",
          "non censored, f = 2", "censored, f = 2"),
        lty=c(1,2,1,2),
        lwd=1,
        col=c(par('fg'), par('fg'), 'red', 'red') )

*

plot(survfit(s ~ f), col=as.numeric(levels(f)))

*

A FAIRE : avec deux facteurs
survfit(s~f1+f2)

A FAIRE : et si f n'est pas un facteur ??? (mettre ça un peu plus loin, quand on parlera de régression d'une variable de survie).

A FAIRE : choisir deux exemples, l'un avec l'autre sans covariables, que l'on utilisera jusqu'à la fin de ce chapitre.

data(lung)
x <- Surv(lung$time, lung$status)
plot(x)

*

plot(survfit(x))

*

Mantel-Cox Log-rank test

On peut demander si une variable de survie dépend d'un facteur.

data(lung)
y <- Surv(lung$time, lung$status)
x <- lung$sex
summary(coxph(y~x))

Rsquare= 0.046   (max possible= 0.999 )
Likelihood ratio test= 10.6  on 1 df,   p=0.00111
Wald test            = 10.1  on 1 df,   p=0.00149
Score (logrank) test = 10.3  on 1 df,   p=0.00131

Modèles à temps discret

A FAIRE

Discrete hazard rate : 
h(aj) = P[  a(j-1) < T <= a(j) | T > a(j-1)  ]
A FAIRE : exprimer la fonction de survie en fonction de h.

A FAIRE
LifeTable : c'est l'équivalent discret de l'estimateur de Kaplan-Meier

Modèle de régression paramétrique

A FAIRE : Introduction

Nous allons maintenant ajouter des variables prédictives.

Certaines ne dépendent pas du sujet (mais dépendent du temps) :
inflation, faux de chomage, température extérieure, etc.

D'autres dépendent du sujet et peuvent dépendre du temps (CSP, age,
pression artérielle) ou pas (sexe).

Méthodes : 
- paramétrique
- semi-paramétrique (régression de Cox) : paramétrique pour la
  partie qui dépend des variables prédictives, non paramétrique
  (Kaplan-Meier) pour le reste.

A FAIRE : classification des différents modèles

Continu : exponentielle (PH et AFT), Weibull (PH et AFT), 
Gompertz (PH), lognormal (AFT), loglogistic (AFT), gamma généralisée (AFT)  

discret : logistic (proportionnal odds), cloglog (PH)

PH (Proportionnal Hazards) : lambda = lambda0(t) * exp(X béta)

AFT (Accelerated Failure Time) : T = exp(-X béta) * exp(z)

A FAIRE : que fait la fonction survreg ???

On peut tout d'abord effectuer ces calculs à la main, en écrivant explicitement la log-vraissemblance

Prod( f(b,x_i) où i non censuré ) * Prod( S(b,x_i) où i censuré )

(il manque des termes correspondant à la probabilité de censure, mais on suppose que la censure est indépendante du résultat de l'expérience, donc on peut l'oublier) à l'aide de la commande optim. Le seul point délicat est le choix des paramètres initiaux : s'il sont mal choisis (j'avais commencé par prendre c(1,1)), la vraissemblance est nulle...

Pour vérifier graphiquement que c'est bien un minimum, on peut faire un dessin (ici, il y a deux paramètres à estimer : s'il n'y en avait qu'un, on ferait un dessin en dimension un (on en a déjà fait), s'il y en avait plus, on aurait plus de problèmes (on pourrait s'en tirer avec une animation ou un treillis))

data(lung)
x <- Surv(lung$time, lung$status)
f <- function (p,t) { dweibull(t,p[1],p[2]) }
S <- function (p,t) { 1-pweibull(t,p[1],p[2]) }
ll <- function (p) {
  time <- x[,1]
  status <- x[,2]
  censored <- 0
  dead <- 1
  # cat(p); cat("\n"); str(time); cat("\n");
  -2*( sum(log(f(p,time[status==dead]))) + sum(log(S(p,time[status==censored]))) )
}

# Quelques estimations du second paramètre
m <- 1     # Ne marche pas
m <- 100
s <- survfit(x)
m <- mean(s$time)
m <- max(s$time[s$surv>.5])

r <- optim( c(1,m), ll )

# Dessinons la log-vraissemblance
myOuter <- function (x,y,f) {
  r <- matrix(nrow=length(x), ncol=length(y))
  for (i in 1:length(x)) {
    for (j in 1:length(y)) {
      r[i,j] <- f(x[i],y[j])
    }
  }
  r
}
lll <- function (u,v) { 
  r <- ll(c(u,v)) 
  if(r==Inf) r <- NA
  r
}
a <- seq(1,1.6,length=50)
b <- seq(100,700,length=50)
ab <- myOuter(a,b,lll)
persp(a,b,ab)

*

op <- par(mfrow=c(3,3))
for (i in seq(0,360,length=10)[-10]) {
  persp(a,b,ab,theta=i)
}
par(op)

*

On pourrait aussi utiliser xgobi :

xgobi(data.frame(
  x=rep(a,1,each=length(b)), 
  y=rep(b,length(a)), 
  z=as.vector(ab))
)

Mais on voit peut-être mieux en restant en dimension deux :

image(a,b,ab)
points(r$par[1],r$par[2],lwd=3,cex=3)

*

Avec les couleurs par défaut, on ne voit rien :

n <- 255
image(a,b,ab, col=topo.colors(n), breaks=quantile(ab,(0:n)/n, na.rm=T))
points(r$par[1],r$par[2],lwd=3,cex=3)

*

Voire même, pour souligner les différences au voisinage du minimum :

image(a,b,ab, col=topo.colors(n), breaks=quantile(ab,((0:n)/n)^2,na.rm=T))
points(r$par[1],r$par[2],lwd=3,cex=3)

*

On pourrait faire les mêmes dessins avec la bibliothèque "lattice".

A FAIRE
library(lattice)
?levelplot
?contourplot

On peut comparer la courbe de Kaplan-Meier avec la courbe théorique :

plot(survfit(x))
curve( 1-pweibull(x,r$par[1],r$par[2]), add=T, col='red', lwd=3, lty=2 )

*

Prenons un exemple et essayons de le modéliser à l'aide d'une distribution exponentielle, de Weibull et Gamma.

ph.mle.weibull <- function (x) {
  f <- function (p,t) { dweibull(t,p[1],p[2]) }
  S <- function (p,t) { 1-pweibull(t,p[1],p[2]) }
  m <- mean(survfit(x)$time)
  ph.mle.optim(x,f,S,c(1,m))
}
ph.mle.exp <- function (x) {
  f <- function (p,t) { dexp(t,p) }
  S <- function (p,t) { 1-pexp(t,p) }
  m <- mean(survfit(x)$time)
  ph.mle.optim(x,f,S,1/m)
}
ph.mle.gamma <- function (x) {
  f <- function (p,t) { dgamma(t,p[1],p[2]) }
  S <- function (p,t) { 1-pgamma(t,p[1],p[2]) }
  m <- mean(survfit(x)$time)
  ph.mle.optim(x,f,S,c(1,1/m))
}
ph.mle.optim <- function (x,f,S,m) {
  ll <- function (p) {
    time <- x[,1]
    status <- x[,2]
    censored <- 0
    dead <- 1
    -2*( sum(log(f(p,time[status==dead]))) + sum(log(S(p,time[status==censored]))) )
  }
  optim(m,ll)
}
eda.surv <- function (x) {
  r <- survfit(x)
  plot(r)
  a1 <- ph.mle.exp(x)$par
  lines( 1-pexp(r$time,a1) ~ r$time, col='red' )
  a2 <- ph.mle.weibull(x)$par
  lines( 1-pweibull(r$time,a2[1],a2[2]) ~ r$time, col='green' )
  a3 <- ph.mle.gamma(x)$par
  lines( 1-pgamma(r$time,a3[1],a3[2]) ~ r$time, col='blue' )
  legend( par("usr")[2], par("usr")[4], yjust=1, xjust=1,
          c(paste("Exponentielle(", signif(a1,2), ")", sep=''), 
            paste("Weibull(", signif(a2[1],2), ", ", signif(a2[2],2), ")", sep=''), 
            paste("Gamma(", signif(a3[1],2), ", ", signif(a3[2],2), ")", sep='')
          ),
          lwd=1, lty=1,
          col=c('red', 'green', 'blue'))
  title(main="Parametric estimation of PH models")
}

data(lung)
x <- Surv(lung$time, lung$status)
eda.surv(x)

*

Pour regarder la qualité de l'estimation, on peut aussi adapter le qqplot.

x <- Surv(lung$time, lung$status)
r <- survfit(x)
a1 <- ph.mle.exp(x)$par
t1 = 1-pexp(r$time,a1)
a2 <- ph.mle.weibull(x)$par
t2 <- 1-pweibull(r$time,a2[1],a2[2])
a3 <- ph.mle.gamma(x)$par
t3 <- 1-pgamma(r$time,a3[1],a3[2])
plot( t1 ~ r$surv, col='red', xlab='sample', ylab='model')
points( t2 ~ r$surv, col='green')
points( t3 ~ r$surv, col='blue' )
abline(0,1)
legend( par("usr")[1], par("usr")[4], yjust=1, xjust=0,
        c(paste("Exponentielle(", signif(a1,2), ")", sep=''), 
          paste("Weibull(", signif(a2[1],2), ", ", signif(a2[2],2), ")", sep=''), 
          paste("Gamma(", signif(a3[1],2), ", ", signif(a3[2],2), ")", sep='')
        ),
        lwd=1, lty=1,
        col=c('red', 'green', 'blue'))
title(main="Parametric estimation of PH models")

*

On peut aussi représenter les résidus (la différence entre l'estimation non paramétrique (Kaplan-Meier) et l'estimation paramétrique de la fonction de survie.

plot( t1 - r$surv ~ t1, col='red', xlab='predicted values', ylab='residuals')
points( t2 - r$surv ~ t2, col='green')
points( t3 - r$surv ~ t3, col='blue' )
abline(h=0, lty=3)
legend( par("usr")[2], par("usr")[4], yjust=1, xjust=1,
        c(paste("Exponentielle(", signif(a1,2), ")", sep=''), 
          paste("Weibull(", signif(a2[1],2), ", ", signif(a2[2],2), ")", sep=''), 
          paste("Gamma(", signif(a3[1],2), ", ", signif(a3[2],2), ")", sep='')
        ),
        lwd=1, lty=1,
        col=c('red', 'green', 'blue'))
title(main="Parametric estimation of PH models")

*

plot( abs(t1 - r$surv) ~ t1, col='red', 
      xlab='predicted values', ylab=expression( abs(residuals) ))
points( abs(t2 - r$surv) ~ t2, col='green')
points( abs(t3 - r$surv) ~ t3, col='blue' )
abline(h=0, lty=3)
legend( par("usr")[2], par("usr")[4], yjust=1, xjust=1,
        c(paste("Exponentielle(", signif(a1,2), ")", sep=''), 
          paste("Weibull(", signif(a2[1],2), ", ", signif(a2[2],2), ")", sep=''), 
          paste("Gamma(", signif(a3[1],2), ", ", signif(a3[2],2), ")", sep='')
        ),
        lwd=1, lty=1,
        col=c('red', 'green', 'blue'))
title(main="Parametric estimation of PH models")

*

A FAIRE : ajouter un test pour vérifier si, sous l'hypothèse d'un modèle
de Weibull, les données sont vraiment différentes d'un modèle exponentiel
(ie, tester si le premier paramètre est vraiment différent de un).
(mettre ce test dans eda.surv et afficher le résultat sur le dessin)

Si on reprend l'exemple simulé utilisé plus haut (Weibull de paramètres .5 et 1), on obtient des valeurs proches des valeurs attendues.

> n <- 200
> x <- rweibull(n,.5)
> y <- rexp(n,1/mean(x))
> s <- Surv(ifelse(x<y,x,y), x<=y)
> ph.mle.weibull(s)$par
[1] 0.4891810 0.9924432

A FAIRE : tests...

A FAIRE : diagnostic plots

Proportional Hazard (PH) models

Si on cherche à prédire une variable de survie à l'aide d'une autre variable, X, on peut utiliser le modèle suivant

lambda(t|X) = lambda(t) * exp(X*b)

où lambda est une Hazard function sous-jacente (correspondant par exemple à un modèle de Weibull) et b un paramètre à estimer. Pour trouver b et les paramètres du modèle sous-jacent, on utilise toujours le maximum de raissemblance.

A FAIRE : exemple

data(lung)
y <- Surv(lung$time, lung$status)
x1 <- lung$age
x2 <- lung$meal.cal

A FAIRE : ajouter b dans le code suivant...
Je n'arrive pas à obtenir d'expression simple et générale de la
log-vraissemblance (j'ai une intégrale de lambda qui apparait...)

f <- function (p,t) { dweibull(t,p[1],p[2]) }
S <- function (p,t) { 1-pweibull(t,p[1],p[2]) }
m <- mean(survfit(y)$time)
ll <- function (p) {
  time <- y[,1]
  status <- y[,2]
  censored <- 0
  dead <- 1
  -2*( sum(log(f(p,time[status==dead]))) + sum(log(S(p,time[status==censored]))) )
}
optim(m,ll)

En fait, la fonction survreg nous donne ce résultat.

A FAIRE : comprendre.
Expliquer en quoi c'est le même résultat.

> summary(survreg(y~x1+x2))
Call:
survreg(formula = y ~ x1 + x2)
                Value Std. Error      z        p
(Intercept)  6.82e+00    0.57891 11.785 4.64e-32
x1          -1.35e-02    0.00809 -1.667 9.55e-02
x2           2.37e-05    0.00018  0.131 8.95e-01
Log(scale)  -2.54e-01    0.06973 -3.641 2.71e-04
Scale= 0.776
Weibull distribution
Loglik(model)= -930.6   Loglik(intercept only)= -932.2
        Chisq= 3.17 on 2 degrees of freedom, p= 0.21
Number of Newton-Raphson Iterations: 4
n=181 (47 observations deleted due to missing)

A FAIRE : pour vérifier aue je comprends bien le sens de ce résultat,
il faudrait faire une simulation (Weibul de paramètres .5 et 3, b=c(1 2 -1))
et la donner à survreg.

A FAIRE : un exemple avec des variables prédictives qui ne dépendent pas du temps
???
A FAIRE : un exemple avec des variables prédictives qui dépendent du temps ???
???

Accelerated Failure Time (AFT) models

On suppose que la variable de survie T est donnée par :

ln(T) = béta X + erreur

où le terme d'erreur suit une certaine distribution (par exemple, si c'est une distribution normale, T est log-normale)

A FAIRE : extreme value distribution with 2 parameters (this yields
the Weibul distribution)

On peut aussi écrire ces modèles :

ln(psi * T) = erreur
psi = exp(-béta X)

On remarquera que psi ne dépend pas de T ; de plus, si psi>1, on meurt plus vite, alors que si psi<1, on meurt plus lentement -- d'où le nom "Accelerated Failure Time".

Ce sont simplement des modèles linéaires généralisés.

A FAIRE 
?glm

A FAIRE

Autres modèles

PCE model (Piecewise Constant Exponential)

C'est un modèle PH dans lequel lambda (la Hazard function sous-jacente) est constante par morceaux (c'est l'expérimentateur qui choisit les instants en lesquels la constante change).

C'est une première analyse, un peu simpliste, que l'on peut faire quand les prédicteurs dépendent du temps.

A FAIRE : traiter un exemple

Cox Proportional Hazard model

Dans un "Proportional Hazard Model",

lambda = lambda0 * exp(b * X)

on peut en fait estimer les coefficients de régression b sans tenir compte de la forme de lambda0, à l'aide d'une vraissemblance partielle.

A FAIRE (pas tout de suite) : je n'ai pas compris d'où vient la
vraissemblance partielle.

A FAIRE

A FAIRE

library(survival)
?survfit
?survreg
?coxph

Proportionnal hazard:  lambda(t|X) = lambda(t) exp(X beta)

Cox proportional hazard model
A popular semi-parametric model for survival analysis -- as efficient as parametric models,
even when those apply.
The model is still lambda(t|X) = lambda(t) exp(X beta), but we do not assign a prescribed form to lambda.

Variables prédictives : elles peuvent dépendre du temps (ça
complique beaucoup les choses).

A FAIRE : différence entre 
PH: Proportional Hazard models
AFT: Accelerated Failure Time models

Évènements multiples

On peut avoir plusieurs évènements dans la série, par exemple plusieurs types de mort, ou un évènement final (mort) éventuellement précédé d'autres évènements non léthaux (AVC, infarctus, diagnostic d'HTA, opération, etc.).

A FAIRE : exemple
S'il y a deux types d'évènements léthaux indépendants, on peut les traiter séparément.
  Mort accidentelle : 1 1 1 3
  Mort non accidentelle : 1 2 4 5 5 6 
  Vie : 7+ 7+ 7+
Pour étudier les morts accidentelles : 
  1  1  1  3  1+ 2+ 4+ 5+ 5+ 6+ 7+ 7+ 7+ 
Pour étudier les morts non accidentelles : 
  1+ 1+ 1+ 3+ 1  2  4  5  5  6  7+ 7+ 7+ 

Si les évènements ne sont pas indépendants (par exemple, mort par
cancer ou par infarctus, qui peuvent avoir ue cause commune
(tabac)), les résultats seront biaisés.

Frailty models (unobserved heterogeneity)

A FAIRE : est-ce le bon endroit ?

On peut aussi imaginer la situation suivante : tous les patients ont une fonction de risque qui a la même forme, à un facteur multiplicatif près. Si on disposait de suffisemment de variables prédictives, on essaierait de prédire ce facteur à l'aide de ces variables. Mais si on n'en a pas, ou pas suffisemment, on ne peut pas bien prédire ce facteur. On s'intéresse alors à des modèles de la forme

lambda(t) = lambda0(t) * u

où u est une variable aléatoire, positive, de moyenne 1. Comme la valeur de u change d'un sujet à l'autre, on ne pourra pas la prédire, mais on pourra au moins estimer sa variance.

A FAIRE : comment fait-on ?

A FAIRE : et sous R ?

frail.fit(eha)          Fits a frailty proportional hazards model
frailty(survival)       (Approximate) Frailty models
mlefrailty.fit(survrec)
                        Survival function estimator for correlated
                        recurrence time data under a Gamma Frailty
                        Model

A FAIRE

Recurrent failure times
library(survrec)

Penalized maximum likelihood
?pspline
?frailty
?ridge

The Kaplan-Meier curves should not cross: it indicates that 
the PH model does NOT hold.

On peut utiliser l'analyse de la survie même si les observations ne
sont pas censurées : les notions (hazard function) et les graphiques
que l'on manipule peuvent être pertinents.

Frailty models (c'est important)

Vincent Zoonekynd
<zoonek@math.jussieu.fr>
latest modification on Wed Oct 13 22:33:13 BST 2004