Statistique non-paramétrique

TP3: estimation en présence de censure

Génération des données

Nous allons étudier des données générées selon la loi de Weibull qui est une loi classique pour modéliser des durées. Sa densité dépend de deux paramètres \(\alpha\) (paramètre de forme) et \(\lambda\) (paramètre d’échelle) qui sont des réels positifs, \[ f_T(t ) = \frac\alpha\lambda \left(\frac t\lambda\right)^{\alpha-1}e^{-\left(\frac t\lambda\right)^{\alpha}}\mathbf 1_{t>0}. \]

1. (à faire chez vous) Représenter \(f_X\) pour différentes valeurs de \(\alpha\) et \(\lambda\). On prétera une attention particulière aux cas \(\alpha<1\), \(\alpha=1\) (loi exponentielle) et \(\alpha>1\).

2. Générer indépendamment, pour \(n=500\), \(T_1,...,T_n\) i.i.d. selon la loi de Weibull de paramètres \(\alpha = 2\) et \(\lambda=1/2\) et \(C_1,...,C_n\) i.i.d. selon une loi exponentielle de paramètre \(\lambda_e=1\). Puis générer les observations associées \[ Y_i = \min(T_i, C_i), i=1,...,n \] et, l’indicateur de non censure, \[ \delta_i = \mathbf 1_{T_i\leq C_i}, i=1,...,n. \] Quel est le taux de censure ? Faire varier le taux de censure en prennant d’autres valeurs du paramètre \(\lambda_e\) de la loi exponentielle (en gardant le même vecteur \(T_1,...,T_n\)).

3. Tracer des estimateurs de la densité :

• des données originales (sans censure)
• des données censurées pour différentes valeurs du paramètre \(\lambda_e\) de la loi exponentielle.

4. Estimer la fonction de répartition :

• des données originales (sans censure)
• des données censurées pour différentes valeurs du paramètre \(\lambda_e\) de la loi exponentielle. Pour calculer l’estimateur de Kaplan-Meier, on pourra s’aider de la fonction survfit() du package survival où coder l’estimateur à la main à partir de la formule du cours (si vous faites les deux vous pouvez comparer les deux courbes). Comparer à la vraie fonction de répartition.

# estimateur de Kaplan-Meier à la main
  KL <- function(Y,delta){
  
  Yordtot = sort.int(Y,index.return=T)
  Yord = Yordtot$x
  deltaord = delta[Yordtot$ix]
  hatSTY = cumprod((1-1/(seq(n,1,-1)))^deltaord)
  list(Yord = Yord, hatSTY = hatSTY)
  } 
# KL calcule l'estimateur en les temps d'évènements
KL1 = KL(Y1,delta1)
KL2 = KL(Y2,delta2)
KL3 = KL(Y3,delta3)
plot(KL1$Yord,1-KL1$hatSTY,type='s',col=2,main="Estimateur de Kaplan-Meier (à la main)",xlab='t',ylab=expression(hat(F)[T](t)))
points(KL2$Yord,1-KL2$hatSTY,col=3,type='s')
points(KL3$Yord,1-KL3$hatSTY,col=4,type='s')
plot(ecdf(Tvec),add=T,do.points=F,verticals=T)
points(sort(Tvec),pweibull(sort(Tvec),alpha,lambda),lwd=2,type='l',col="darkgreen",lty=2)
legend('bottomright',c(expression(F[T]), "sans censure","lambda_e = 0.5","lambda_e=1","lambda_e=3"),col=c("darkgreen",1,2,4,3),lty=c(2,rep(1,4)),lwd=c(2,rep(1,4)))

# Package survival
library(survival)
sample = c(rep(1,n),rep(2,n),rep(3,n))
fit = survfit(Surv(c(Y1,Y2,Y3),c(delta1,delta2,delta3))~sample)
plot(fit,fun='F',col=2:4,main="Estimateur de Kaplan-Meier (avec survfit)",xlab='t',ylab=expression(hat(F)[T](t)))
plot(ecdf(Tvec),add=T,do.points=F,verticals=T)
points(sort(Tvec),pweibull(sort(Tvec),alpha,lambda),lwd=2,type='l',col="darkgreen",lty=2)
legend('bottomright',c(expression(F[T]), "sans censure","lambda_e = 0.5","lambda_e=1","lambda_e=3"),col=c("darkgreen",1,3,2,4),lty=c(2,rep(1,4)),lwd=c(2,rep(1,4)))

5. Estimer la densité de \(T\) à partir des différents échantillons. On pourra utiliser la fonction density en attribuant un poids à chaque individu (option weights) pour mettre en place la correction de censure.

6. Données leukemia. Estimer la survie des patients atteints de leucémie avec et sans traitement par chimiothérapie.