9月8日（月）午後

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

それぞれの結果に対する理由の説明を試みたが、思ったよりも難しくて、ちょっといい説明が分からなった。
ミーティングでの話題にするとして、一旦保留にすることにした。

結果を見て思ったのは、思ったよりもバラつきが大きいということだった。
ケース・コントロール研究は、特に発生率が非常に小さくコホート研究を行うと無駄が多い、という場合によく適応される。
しかし、この結果をみると、発生率が小さくて症例数を十分に集められないような疾患だと、研究をするたびにかなりの違いが生じそうである。
例えば、Ａという因子とＢという因子について、どちらがより影響が大きいか、というような検討をするにしても、これだけバラつきがあっては信頼度はかなり低いと思われる。
ましてやこれは、原集団の同定が完全にうまくいっていると仮定した結果である。
僕が予想していたよりも、ケース・コントロール研究による結果はあてにならないようである。
やはり手間のかかるコホート研究のほうが、結果の信頼性は高いらしい。

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

この調子で、ケース・コントロール研究の他の種類のシミュレーションもやってみる。

次はリスク期間が終わった後に行ってリスクの比を調べる、累積症例対照研究について扱う。

N人の未発症の集団のうち、曝露されている人数をeとし、ある期間における曝露群の発症のリスクをp,非曝露群の発症のリスクをqとおく。
このとき、期間後に発症している人からn1人を抽出して症例群、発症していない人からn2人を抽出して対照群とし、オッズを求める。

とりあえず作ったのが以下。

ex<-function(N,e,p){
r<-runif(e,0,1)
r[r>=1-p]<-1
r[r<1-p]<-0
return(rbind(r,1))
}

n.ex<-function(N,e,q){
r<-runif(N-e,0,1)
r[r>=1-q]<-1
r[r<1-q]<-0
return(rbind(r,0))
}

all<-function(N,e,p,q){
cbind(ex(N,e,p),n.ex(N,e,q))
}

case.sample<-function(N,e,p,q,n1){
a<-all(N,e,p,q)
b<-a[,which(a[1,]==1)]
sample.num<-sample(ncol(b),n1)
return(b[,sample.num])
}

control.sample<-function(N,e,p,q,n2){
a<-all(N,e,p,q)
b<-a[,which(a[1,]==0)]
sample.num<-sample(ncol(b),n2)
return(b[,sample.num])
}

myf<-function(c){
length(which(c))
}

odds<-function(N,e,p,q,n1,n2){
a<-case.sample(N,e,p,q,n1)
b<-control.sample(N,e,p,q,n2)
x<-(myf(a[2,]==1)*myf(b[2,]==0))/(myf(b[2,]==1)*myf(a[2,]==0))
return(x)
}

box<-function(N,e,p,q,n1,n2){
x<-c()
for(i in 1:10000){
x[i]<-odds(N,e,p,q,n1,n2)
}
boxplot(x)
}

しかし、案の定計算が遅すぎて使い物にならないので、速度の改善を試みた。

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){
v<-m[1,c(((i-1)*e+1):(i*e),(10000*e+(N-e)*(i-1)+1):(10000*e+(N-e)*i))]
a<-m[,sample(which(v==2),n1)]
b<-m[,sample(which(v!=2),n2)]
x[i]<-(length(which(a[2,]==2))*length(which(b[2,]!=2)))/(length(which(b[2,]==2))*length(which(a[2,]!=2)))
}
boxplot(x)
}

一応ここまでつくってみたが、まだちゃんとした箱ひげ図の作成には成功していない。
どこかにミスがあるようなので、明日もう少し探してみることにする。

効率化を考えるのはなかなかハマる。
後から見たら訳が分からないものになっていそうだが。

以上