9月9日（火）午前

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

関数のバグはとれた。

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

これで、5分ほど時間はかかるものの箱ひげ図は描けるようになった。

集団が有限のため、密度症例対照研究の場合よりもパラメータの取り方に制限があってややこしい。
その上、なかなかよく分からない結果も出ている。

ここで、累積症例対照研究の検証の前に、前回の密度症例対照研究の結果についてもう少しまとめることにする。

密度症例対照研究のシミュレーションまとめ

まず、比べる箱ひげ図を横に描く方法があったようなので、それを使ってみた。

control.m<-function(e1,N){
r<-runif(10000*N,0,1)
r[r>=1-e1]<-1
r[r<1-e1]<-0
x<-matrix(r,nrow=10000,ncol=N)
return(x)
}

case.m<-function(e1,p,q,n){
e2<-e1*p/(e1*p+(1-e1)*q)
r<-runif(10000*n,0,1)
r[r>=1-e2]<-1
r[r<1-e2]<-0
x<-matrix(r,nrow=10000,ncol=n)
return(x)
}

box.v<-function(e1,p,q,n,N){
x<-c()
c<-case.m(e1,p,q,n)
d<-control.m(e1,N)
for(i in 1:10000){
a<-c[i,]
b<-d[i,]
x[i]<-(length(a[a==1])*length(b[b==0]))/(length(b[b==1])*length(a[a==0]))
}
return(x)
}

a<-box.v(0.1,0.02,0.01,1000,1000)
b<-box.v(0.5,0.02,0.01,1000,1000)
c<-box.v(0.9,0.02,0.01,1000,1000)

boxplot(a,b,c,ylim=c(0.5,4),names=c("0.1","0.5","0.9"),
main="p=0.02,q=0.01,N=1000,n=1000",ylab="OR",xlab="e1")

このほうが圧倒的にわかりやすい。
今回すでになんとなくの傾向が分かっているので全部描き直すことはしないが、今度からは比べるときは並べて描くことにする。

各パラメータにおける分散について、また午後にまとめる。

以上