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}. \]

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

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.

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.

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,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)))

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

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")

Statistique non-paramétrique

TP3: estimation en présence de censure

Angelina Roche

Génération des données