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

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

結局土日は研究室には来なかった。
けど、土日はある程度休んでもいいんじゃないかと思い始めてきたので、進行度がやばい！ってなったらまた考えることにする。

一応土日の間ちょくちょくケース・コントロール研究のモデル作成について考えていたのだが、なんかしっくり来るのが見つからない。
とりあえず最低限度思いついたものをつくってみることにする。

まずは、本の中で「密度対症症例研究」とされているものから始める。

ある疾患のリスク因子に対する原集団の曝露率をe1、曝露群の疾患の発生率をp、非曝露群の発生率をqとおく。
このとき、疾患にかかっている人の曝露率e2=e1p/{e1p+(1-e1)q}と表せる。
（※本当はここに至るまでもシミュレーションできたほうが面白いと思ったが、密度対症症例研究の場合は原集団の規模は無限大と考える必要があるので、いいモデリングが思いつかなかった。）

対照群として原集団からN人、疾患群として疾患にかかっている人からn人を抽出してオッズを求めたとき、どの程度の正確さが出るのかをシミュレーションで検証してみる。

ef<-function(e){
r<-runif(1,0,1)
if(r<e)x<-1
else x<-0
return(x)
}

control.v<-function(e1,N){
x<-c()
for(i in 1:N)
x[i]<-ef(e1)
return(x)
}

case.v<-function(e1,p,q,n){
e2<-e1*p/(e1*p+(1-e1)*q)
x<-c()
for(i in 1:n)
x[i]<-ef(e2)
return(x)
}

odds<-function(e1,p,q,n,N){
a<-case.v(e1,p,q,n)
b<-control.v(e1,N)
x<-(length(a[a==1])*length(b[b==0]))/(length(b[b==1])*length(a[a==0]))
y<-p/q
return(c(x,y))
}

> odds(0.5,0.02,0.01,100,100)
[1] 1.909091 2.000000

一応それらしいものはできた。
今回の検証には、箱ひげ図を利用してみたいと思う。

最初は

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

こんなかんじでやっていたのだが、処理が遅かったので、勉強もかねて改善を試みた。

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<-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]))
boxplot(x)
}

こんなんで速くなるのか…？と半信半疑だったが、20秒くらいかかっていた処理が1秒くらいになってびっくりした。
これからは工夫してコードを書くことも意識していこうと思った。

いろいろ条件を変えて描いてみた結果は以下の通り。

◎原集団の曝露率e1について(p=0.02,q=0.01,n=1000,N=1000とした)

e1=0.1

f:id:kouri_don:20140908125357p:plain