9月4日（木）午後 - 2014年度マイコースプログラム

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

${r_i}$ について、

ri<-function(d,t){
if (length(d[,1][d[,1]>=t])==0)return(1)
else return(length(d[,1][d[,1]>=t])
}

としておかないと、length(d[,1][d[,1]>=t])==0のときにS.hatがNaNとなってしまうことに気が付いた。（「1」にはどんな定数を入れてもよい。）

というわけで、一応これに変更することにしたが、よく考えるとそもそもカプラン・マイヤー法は

$\frac{d_i}{r_i}$

が死亡率を近似していると仮定できるくらいに ${r_i}$ が大きくないと成立しない。

午前の時点では気づいていなかったが、 ${t=l}$ のときも十分な ${r_i}$ の大きさが保てるような ${p}$ と ${l}$ を設定する必要があるようである。

まずは、 ${E}$ の平均値について調べた。（同じ条件で100回サンプルをとった。）

Em<-function(p,q,l,N){
x<-c()
for(i in 1:100) x<-c(x,E(p,q,l,N))
return(mean(x))
}

xq<-seq(0,0.5,0.05)

yq<-function(j){
yq<-c()
for(i in xq)yq<-c(yq,Em(0.01*j,i,20,50))
return(yq)
}

yq1<-yq(1)
yq2<-yq(2)
yq3<-yq(3)
yq4<-yq(4)
yq5<-yq(5)

plot(xq,yq1,xlim=c(0,0.5),ylim=c(0,0.1), 
ann=F,type="l", lty=1,col=1)

par(new=T)
plot(xq,yq2,xlim=c(0,0.5),ylim=c(0,0.1), 
ann=F,type="l", lty=1,col=2) 

par(new=T) 
plot(xq,yq3,xlim=c(0,0.5),ylim=c(0,0.1), 
ann=F,type="l", lty=1,col=3)

par(new=T) 
plot(xq,yq4,xlim=c(0,0.5),ylim=c(0,0.1), 
ann=F,type="l", lty=1,col=4)

par(new=T) 
plot(xq,yq5,xlim=c(0,0.5),ylim=c(0,0.1), 
,type="l", lty=1,col=5,xlab="q 脱落率" , ylab="E" ,main="p,qによるEの変化")

legend(0.01,0.095,c("p=0.01","p=0.02","p=0.03","p=0.04","p=0.05"),lty=1,col=1:5)

${N=50}$ のとき

f:id:kouri_don:20140904200140p:plain

${N=100}$ のとき

f:id:kouri_don:20140904200207p:plain

（前述の理由から、 ${l=20}$ は固定とした。）

本当は3次元ヒストグラムをつくってみたかったのだが、パッケージのダウンロードがうまくいかなかったので、苦肉の策としてこのようなグラフになった。

考察としては、

・ ${p}$ , ${q}$ が大きいほど ${E}$ は大きくなる。
・ ${p}$ , ${q}$ が ${E}$ に及ぼす影響は、どちらも同じくらいに思われる。
・ ${N}$ が大きいほど ${E}$ は小さくなる。

これらはすべて、推定の精度に重要なのは ${r_i}$ の大きさであることを示している。
脱落率が高いとしても、 ${r_i}$ が十分に大きければ推定が不正確であるとはいえないと分かる。（死亡と脱落が独立して起こると仮定できればだが）
逆に言えば、いかに脱落率が低いデータが得られたとしても、死亡率が高かったり、対象者の数が少なかったりすれば推定の信憑性が高いとはいえない。

続いて、 ${E}$ の分散。（同じく100回サンプルをとった。）

Ev<-function(p,q,l,N){
x<-c()
for(i in 1:100) x<-c(x,E(p,q,l,N))
return(var(x)*99/100)
}

xq<-seq(0,0.5,0.05)

yq<-function(j){
yq<-c()
for(i in xq)yq<-c(yq,Ev(0.01*j,i,20,50))
return(yq)
}

yq1<-yq(1)
yq2<-yq(2)
yq3<-yq(3)
yq4<-yq(4)
yq5<-yq(5)

plot(xq,yq1,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=1)

par(new=T)
plot(xq,yq2,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=2) 

par(new=T) 
plot(xq,yq3,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=3)

par(new=T) 
plot(xq,yq4,xlim=c(0,0.5),ylim=c(0,0.01), 
ann=F,type="l", lty=1,col=4)

par(new=T) 
plot(xq,yq5,xlim=c(0,0.5),ylim=c(0,0.01), 
,type="l", lty=1,col=5,xlab="q 脱落率" , ylab="Eの分散値" ,main="p,qによるEの分散")

legend(0.01,0.0095,c("p=0.01","p=0.02","p=0.03","p=0.04","p=0.05"),lty=1,col=1:5)

${N=50}$ のとき

f:id:kouri_don:20140904194922p:plain

${N=100}$ のとき

f:id:kouri_don:20140904194931p:plain

これはよくわからない結果になった。

${N}$ の値による有意差はほとんど見られない。
${N}$ が大きいほど分散は小さくなりそうなものだが…？

また、 ${p}$ の値が小さいほど ${q}$ が大きくなった時の分散の変動が大きい。
これも感覚からすると意外な結果である。

${p}$ と ${q}$ でとっている変域が違うので、 ${p}$ と ${q}$ の対称性についてはちょっとよくわからない。

どうやら分散についてはまだ追及不足なようなので、明日ももう少し分析してみることにする。

以上