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\)).

n = 500
alpha = 2
lambda = 0.5
lambdaexp1 = 1
Tvec = rweibull(n,alpha,lambda)
C1 = rexp(n,lambdaexp1)

Y1 = pmin(Tvec,C1)
delta1 = Tvec<=C1

# Taux de censure 
1-mean(delta1)
## [1] 0.348
lambdaexp2 = 0.5
C2 = rexp(n,lambdaexp2)

Y2 = pmin(Tvec,C2)
delta2 = Tvec<=C2

# Taux de censure 
1-mean(delta2)
## [1] 0.208
lambdaexp3 = 3
C3 = rexp(n,lambdaexp3)

Y3 = pmin(Tvec,C3)
delta3 = Tvec<=C3

# Taux de censure 
1-mean(delta3)
## [1] 0.652
  1. Tracer des estimateurs de la densité :
plot(density(Tvec,bw='ucv'),main='Densité des données observées',ylim=c(0,3))
points(density(Y1,bw='ucv'),col=2,type='l')
points(density(Y2,bw='ucv'),col=3,type='l')
points(density(Y3,bw='ucv'),col=4,type='l')
legend('topright',c("sans censure","lambda_e = 1","lambda_e=0.5","lambda_e=3"),col=1:4,lty=rep(1,4))

# On observe une forte différence entre les densités des données pour différentes valeurs du paramètre lambda_e.
  1. Estimer la 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,3,2,4),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)))

  1. 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.
# estimateur de Kaplan-Meier de la censure
KLC <- function(Y,delta){
  Yordtot = sort.int(Y,index.return=T)
  Yord = Yordtot$x
  deltaord = delta[Yordtot$ix]
  hatSCY = cumprod((1-1/(seq(n+1,2,-1)))^(1-deltaord))
  list(Yord = Yord, hatSCY = hatSCY, deltaord=deltaord)
} 
KLC1 = KLC(Y1,delta1)
poids1 = KLC1$deltaord/KLC1$hatSCY
plot(density(KLC1$Yord,bw='ucv',weights = poids1/sum(poids1)),main='Estimation à noyau de la densité de T',col=2,ylim=c(0,2))
KLC2 = KLC(Y2,delta2)
poids2 = KLC2$deltaord/KLC2$hatSCY
points(density(KLC1$Yord,bw='ucv',weights = poids2/sum(poids2)),col=3,type='l')
KLC3 = KLC(Y3,delta3)
poids3 = KLC3$deltaord/KLC3$hatSCY
points(density(KLC3$Yord,bw=0.1,weights = poids3/sum(poids3)),col=4,type='l')
points(density(Tvec,bw='ucv'),type='l')
points(sort(Tvec),dweibull(sort(Tvec),alpha,lambda),type='l',lwd=2,col="darkgreen",lty=2)
legend('topright',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)))

  1. Données leukemia. Estimer la survie des patients atteints de leucémie avec et sans traitement par chimiothérapie.
leukemia.surv <- survfit(Surv(time, status) ~ x, data = aml) 
plot(leukemia.surv, lty = 2:3) 
legend(100, .9, c("Maintenance", "No Maintenance"), lty = 2:3) 
title("Kaplan-Meier Curves\nfor AML Maintenance Study")