9月4日（木）午前

生存曲線のシミュレーション

時刻を0,1,2…とし、時間1の間の死亡率を、脱落率をとするモデルを考えると、

Surv.time<-function(p,q,l){
cnt<-1
repeat{ 
if (cnt==l){x<- c(cnt,0);break}
else if (runif(1)< p) {x<- c(cnt,1);break}
else if (runif(1)< q) {x<- c(cnt,0);break}
else cnt<-cnt+1
}
return(x)
}# 生存時間のシミュレーション

A.St<-function(p,q,l){
Surv.time(p,q,l)
}

data.sample<-function(p,q,l,N){
mat<-numeric(0)
for(i in 1:N)
mat<-rbind(mat,A.St(p,q,l))
return(mat)
}# N人を対象にしたデータを作成

data.sample(0.05,0.01,100,100)

> data.sample(0.05,0.01,100,100)
       [,1] [,2]
  [1,]    5    1
  [2,]    6    1
  [3,]   36    0
  [4,]    1    1
  [5,]   29    1
  [6,]    1    1
  [7,]   21    1
  [8,]   27    1
  [9,]   26    1
 [10,]   11    1
（以下略）

このように、コホート研究のシミュレーションができる。

このシミュレーションの場合、理論上は

となるはずである。

カプラン・マイヤー法によって求めた推定値との誤差の指標は、

とできる。

ri<-function(d,t){
length(d[,1][d[,1]>=t])
}#rの値を出す関数

di<-function(d,t){
cnt<-0
for(i in 1:nrow(d)){
if(d[i,1]==t&&d[i,2]==1)cnt<-cnt+1
else cnt<-cnt
}
return(cnt)
}#diを出す関数

S.hat<-function(d,t){
x<-1
for(i in 0:t){
x<-x*(1-di(d,i)/rt(d,i))
}
return(x)
}#カプラン・マイヤー法に適用

S.true<-function(p,t){
(1-p)^t
}

E<-function(p,q,l,N){
d<-data.sample(p,q,l,N)
x<-0
for(i in 1:l)x<-x+(S.hat(d,i)-S.true(p,i))^2
return(x/l)
}

ここで試しに、の場合について、S.hatとS.trueのそれぞれを実際にグラフにしてみた。

この規模と脱落率だと、だいたい同じくらいになっているはず。

x<-0:100

A.d<-data.sample(0.05,0.01,100,100)

y.Sh<-c()
for(i in 0:100) y.Sh<-c(y.Sh,S.hat(A.d,i))

y.St<-c()
for(i in 0:100) y.St<-c(y.St,S.true(0.05,i))
[f:id:kouri_don:20140904123135p:plain]
 
plot(x,y.Sh,xlim=c(0,100),ylim=c(0,1), 
ann=F,type="s", lty=1) 

par(new=T)                                          
plot(x,y.St,xlim=c(0,100), ylim=c(0,1),
xlab="時間" , ylab="生存率" ,main="生存曲線",type="s",lty=2 ) 

legend(70,0.95,c("計算値","理論値"),lty=1:2)

結果は下の通り。

これはあくまで一例ではあるが、カプラン・マイヤー法の正当性を見せつけられたと言ってよいのではないか。

あとは、パラメーターをいろいろ変えてみつつ、の平均値・分散がどの程度変化するか記述していこうと思う。

成果は午後にまとめる。

以上