9月10日（水）午前

累積症例対照研究のシミュレーション 3

先生からの助言を受けて、先生のおっしゃる通り1～4の4つの数字で集団を表すというのはもっとも効率のよいアイデアだと思ったので、それで箱ひげ図を描く関数を作り直してみた。

下は、同じ処理を今までのものにさせた場合との比較である。

mat<-function(N,e,p,q){
a<-N*10000
b<-e*10000
r<-runif(a,0,1)
r[1:b][r[1:b]>=1-p]<-2
r[(b+1):a][r[(b+1):a]>=1-q]<-2
s<-rep(c(1,0),c(b,a-b))
return(rbind(r,s))
}

box<-function(N,e,p,q,n1,n2){
m<-mat(N,e,p,q)
x<-c()
for(i in 1:10000){
n<-m[,c(((i-1)*e+1):(i*e),(10000*e+(N-e)*(i-1)+1):(10000*e+(N-e)*i))]
a<-n[,sample(which(n[1,]==2),n1)]
b<-n[,sample(which(n[1,]!=2),n2)]
x[i]<-(length(which(a[2,]==1))*length(which(b[2,]==0)))/(length(which(b[2,]==1))*length(which(a[2,]==0)))
}
boxplot(x)
}

system.time(box(10000,5000,0.1,0.05,50,50))

> system.time(box(10000,5000,0.1,0.05,50,50))
   ユーザ   システム       経過  
     19.77       3.57      49.89 

new.box<-function(N,e,p,q,n1,n2){
probs <- c(e*p,e*(1-p),(1-e)*q,(1-e)*(1-q))
s <- matrix(sample(1:4,N*10^4,prob=probs,replace=TRUE),
ncol=10^4,nrow=N)
odds.f<-function(x){
y <- sample(c(x[x==1],x[x==3]),n1)
z <- sample(c(x[x==2],x[x==4]),n2)
return((length(y[y==1])*length(z[z==4]))
/(length(z[z==2])*length(y[y==3])))
}
boxplot(apply(s,MARGIN=2,odds.f))
}

system.time(new.box(10000,0.5,0.1,0.05,50,50))

> system.time(new.box(10000,0.5,0.1,0.05,50,50))
   ユーザ   システム       経過  
     21.74       0.62      22.73 

劇的に速くなっていることがわかる。
一度使ってみたいと思っていたapply関数が利用可能になったことも大きいのかもしれない。

同じ処理でもこれだけ違うとは、プログラミングは奥が深い。

箱ひげ図を使ってだいたいの傾向は分かったので、あとは平均・分散の値をパラメータに応じてプロットしてみる。

成果は午後にまとめる。

以上