simulations des courbes du taux de hasard conditionnel

**SikoKossi** · 16/01/2014, 12h44

Bonjour à tous, je suis un nouveau sur le forum et un non compétent en utilisation de R. Je voudrais une aide pour résoudre les problèmes que je rencontre en voulant effectuer une simulation avec R

En effet, je voulais comparer par simulations les courbes du taux de hasard conditionnel lambda(t|z) et de son estimateur lambdan(t|z) et obtenir la position du maximum s'il existe.
Nous simulons le modèle Log(T)=Z+U avec U suit N(0,1), Z suit N(1, 0.5).

Données simulées :
T=durée de vie
Z=covariable
T est soumise à une censure à droite par C. On observe donc (X,delta) avec
X=min(T,C)
delta=1 si T<=C et delta=0, sinon

Je joins en annexe en pdf la formule de l’estimateur.

Problèmes :
1) Nous n'obtenons pas de courbes estimées stables pour N simulations avec un échantillon de taille n
2) Il y a un décalage total des deux courbes

Code:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
size<-100           #taille de chaque échantillon
nsimul<-500		  # nombre de simulations
mu=2.51			# sera varié pour obtenir un taux de censures variable
 
 
##### Début de simulation
 
simul0<-function(n=10,z=0,mu=0){
 
# Simulation des variables explicatives
Z=rnorm(n,1,0.5)
 
# Simulation des temps Ti conditionnels à Z du model log(Ti)=Zi+Ui
U=rnorm(n)
V=Z+U
T=exp(V)
 
# Simulation des temps de censure
#lc=0.08
#C=rexp(n,lc)
C=rlnorm(n,mu,0.6)
 
# Données finales
X=rep(0,n)
delta=rep(0,n)
for (i in 1:n)
{
  X[i]=min(c(T[i],C[i]))
  delta[i]=1*(T[i]<C[i])   
}
 
#taux de censure
delt=delta[delta==0]
taux.censure=length(delt)/length(delta)
 
# tri en ordre croissant
o=order(X)
X=X[o]
Z=Z[o]
delta=delta[o]
 
#### Le programme de calcul de l'estimateur du taux de hasard conditionnel
 
## Calcul des Poids de Nadaraya-Watson W(n,hn,z)
 
# Definition du noyau d'Epanechnikov Kh
Kh = function(x,h=1)
{
  return(3*(1-(x/h)^2)/4*(abs(x/h)<=1))
}
 
# Définition du noyau de lissage Na
Na = function(x,a=1)
{
  return(dnorm(x,0,a))
}
 
# Calcul de la fenêtre optimal an=hn. On peut utiliser la fonction "density" de Z
# hn=density(Z)$bw
hn=sd(Z)*(3/(4*n))^(1/5)
an=hn
 
# fonction W(n,h,z)
NW=function(h=hn,n=10,z=z){
p=rep(0,n)
for (i in 1:n){
   p[i]=Kh(z-Z[i],h)
}
w=rep(0,n)
   for (i in 1:n){
       w[i]=p[i]/sum(p) 
       }
w=w[o]
 return(w) 
}
 
# Poids des individus à risques
PIR=function(z,n=10){
R=NW(h=hn,n,z)
R=R[o]
PIR=rep(0,n)
for(i in 1:n){
PIR[i]=sum(R[i:n])}
return(PIR)
}
 
 
# Lissage en temps
lissage=function(t,n=10,a=hn){
liss=rep(0,n)
for (i in 1:n){
    liss[i]=Na(t-X[i],a=hn)}
liss=liss[o]
return(liss)
}
 
# Fonction lambdan(t|z)
lambdan=function(t=1,z=z,n=10){
nw=NW(h=hn,n,z)
Liss=lissage(t,n)
p.risk=PIR(z,n)
Num=nw*delta*Liss
Terme=Num/p.risk
Terme<-Terme[!(Terme==Inf)]
Resultat=sum(Terme,na.rm=T)
return(Resultat)
}
 
# Fonction lambda(t|z) simulée de la loi de T avec log(T)=Z+U, U suit N(0,1)
lambda=function(t,z=0){
s<-1.5
f=(1/(s*t))*dnorm((log(t)-z)/s)
Survie=1-pnorm((log(t)-z)/s)
#S=pnorm((z-log(t))/s)
 
return(f/Survie)
}
 
#### Courbes obtenues avec indication de la valeur maximum
 
# Courbe de lamban(t|z=z)
t=seq(0,3,length=100)
y=rep(0,length(t))
v=rep(0,length(t))
 
for (i in 1:length(t)){
y[i]=lambdan(t[i],z=z,n)
v[i]=lambda(t[i],z=z)
}
plot(t,y,type="l",col="red",xlab="temps",ylab=paste("lambdan(t|z=",z,")"),ylim=c(0,max(y)+0.5),
main=paste("Conditional Hazard Rate estimate given z=",z))
lines(t,v, col="blue",type="l",lty=2)
t=seq(0,3,length=100)
F<-function(t){
lambdan(t,z=z,n=n)
}
max=optimize(F,interval=c(0,1),maximum=TRUE)
M=max$maximum
ordmax=max$objective
abline(v=M,lty=3)
#abline(c(M,0),c(M,ordmax),lty=3)
 
c(tauxcensure=taux.censure, HR.moy=mean(y),abscis.max=M, maxi=ordmax)
}
 
simul0(100,0.7,2.51)
 
SIMUL1<-function(nsimul,size,z){
 
x<-replicate(nsimul,simul0(size,z,mu))
list(ValeursMoyennes=rowMeans(x), variances=apply(x,1,var))
 
}
 
SIMUL1(100,100,0.7)

Merci d'avance.