9月10日（水）午後
累積症例対照研究のシミュレーション まとめ
まず、「ロスマンの疫学」に累積症例対照研究は、発症する人が多くなるとオッズ比を過大評価としてしまうと書いてあったので、中央値を見てその確認をしてみる。
p/q=2のとき（グラフのx軸はpとなっているが、qの間違いである。）
確かに発症者が多い方が中央値が大きくなっている。
これはp>qのため、発症する人が多いほど、対照群の曝露率が過小評価されてしまうためである。
ちなみにp/q=1のとき
ぴったり一定になった。
密度症例対照研究のときは、p/qが一定であれば分散の値は変わらなかったが、今回は（これもグラフのx軸はpとなっているが、qの間違いである。）
このように、qが大きくなると分散は大きくなった。
気になったのでp/q=1のときも調べてみると、
こちらはバラバラである。
とここで、昨日の分と合わせて、これらの現象について詳しい説明をするには、二項分布についてちゃんと理解する必要があるな…と思い勉強を始めて、ついでに今までなんとなくスルーしてきた正規分布、ポアソン分布、指数分布等についても勉強した。
なぜ今までこんな重要な概念をスルーしてたんだろうと愕然とした。
このまま一気に疑問に思っていたことを理解してしまおう…と思っていろいろ読んでいるうちに今日は用事がある時間になってしまったので、また明日まとめることにする。
以上
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関数が利用可能になったことも大きいのかもしれない。
同じ処理でもこれだけ違うとは、プログラミングは奥が深い。
箱ひげ図を使ってだいたいの傾向は分かったので、あとは平均・分散の値をパラメータに応じてプロットしてみる。
成果は午後にまとめる。
以上
9月9日（火）午後
密度症例対照研究のシミュレーションまとめ 2
各パラメータにおける分散についてグラフにしてみた。
◎e1について
e1が0.5に近い値をとっているときのほうが分散の値が小さいことは確かなようである。
◎p/qについて
p/qの値が大きいほうが分散の値は大きくなった。
ちなみに、q=2*pとして、pの値を変えてプロットした場合
規則性は見られなかった。大事なのはp/qの値のようである。
◎n,Nについて
このグラフの特徴といえるのは、nの値に対する傾きが、Nの値によらずほぼ同じことだと思う。
すなわち、Nとnの値は、どちらも大きくなれば分散を小さくするようにはたらくものの、そのはたらきは独立している…？ということか。
N,nの両方とも一定の大きさをとらないと、信頼の得られる結果は得られないようである。
どうしてそのような結果になるのか、という背景に関する考察だが、
1:1の割合でＡかＢの事象が起こる、というのと、1:9の割合で起こる、というのを何度も繰り返した場合、前者のほうが真の割合に収束するのが早く、分散が小さい。
いま、オッズ比の計算は、
（症例群の曝露者数）/（症例群の非曝露者数）を（対照群の曝露者数）/（対照群の非曝露者数）で割って出しているが、
真の値が（症例群の曝露者数）＝（症例群の非曝露者数）かつ（対照群の曝露者数）＝（対照群の非曝露者数）に近い値のほうが、より収束するのが早いと考えられる。
分散の値が最少になるのがe1=0.5のときよりもやや左よりなのも（p>qなので、（症例群の曝露者数）＝（症例群の非曝露者数）により近くなるから）、e1=0.1のときよりもe1=0.9のときのほうが分散の値が大きいのも（p>qなので、（症例群の曝露者数）＞（症例群の非曝露者数）がより顕著になるから）、これで説明がつく。
同じように、p,qの結果については、p+qを大きくしても（症例群の曝露者数）/（症例群の非曝露者数）の比率は変わらないが、p/qの値が大きくなればより（症例群の曝露者数）＞（症例群の非曝露者数）が顕著になることと一致している。
N,nについては、Nが影響しているのはオッズの分母のみ、nが影響しているのはオッズの分子のみで、それぞれの影響に相関がないから、このような結果になったと考えられる。
箱ひげ図を描いたときに、N=1000、n=10のときのほうがN=10,n=1000のときよりもオッズ比が高く出て「？？」となったが、これも気合でなんとなく説明ができて、
N=1000.n=10のときのオッズ比のばらつきは、ほとんど分子の値の大きな変動によるものである。
ここで、いま分子（症例群の曝露者数）/（症例群の非曝露者数）の真の値は、分子のほうが大きい状態になっている。
ここから同じ程度でばらつきが生じるとき、（症例群の曝露者数）/（症例群の非曝露者数）の値は小さくなりにくいが、大きくなりやすい。
一方で、N=10,n=1000のときのオッズのばらつきは、ほとんどが分母（対照群の曝露者数）/（対照群の非曝露者数）の値のばらつきによるものだが、
これは真の値を1に設定してあるので、小さい方向へも大きい方向へも同じ程度にばらつく。
このような理由で、N=1000、n=10のときのほうがN=10,n=1000のときよりもオッズ比が高く出たと考えられる。
はじめはちんぷんかんぷんだったが、一つ一つ考えれば案外分かってくることもあるみたいである。
ただ、シミュレーションすることよりも、その結果を整理して考察を加えることのほうがはるかに難しいと感じた。
明日は累積症例対照研究に戻って、今回のような調子で成果をまとめようと思う。
以上
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")
このほうが圧倒的にわかりやすい。
今回すでになんとなくの傾向が分かっているので全部描き直すことはしないが、今度からは比べるときは並べて描くことにする。
各パラメータにおける分散について、また午後にまとめる。
以上
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)
}一応ここまでつくってみたが、まだちゃんとした箱ひげ図の作成には成功していない。
どこかにミスがあるようなので、明日もう少し探してみることにする。
効率化を考えるのはなかなかハマる。
後から見たら訳が分からないものになっていそうだが。
以上
9月8日（月）午前
密度対症症例研究のシミュレーション
結局土日は研究室には来なかった。
けど、土日はある程度休んでもいいんじゃないかと思い始めてきたので、進行度がやばい！ってなったらまた考えることにする。
一応土日の間ちょくちょくケース・コントロール研究のモデル作成について考えていたのだが、なんかしっくり来るのが見つからない。
とりあえず最低限度思いついたものをつくってみることにする。
まずは、本の中で「密度対症症例研究」とされているものから始める。
ある疾患のリスク因子に対する原集団の曝露率を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
e1=0.5
e1=0.9
e1は0.5に近い方が正確な値が出やすいといえる。
◎p/qについて（e1=0.5,n=1000,N=1000とした)
p/q=1.1のとき
p/q=2のとき
p/q=20のとき
p/qの値が小さい方が正確な値が出やすいといえる。
なお、p=0.2,q=0.01のときとp=0.02,q=0.001のときとで比べた場合には、有意な差は見られなかった。
◎n,Nについて（e1=0.5,p=0.02,q=0.01とした）
n,Nのそれぞれが大きいほど正確な値が出るであろうことは自明なので、Nとnの対称性について検証した。
n=10,N=1000のとき
n=1000,N=10のとき
何度やってもだいたいこのような結果になった。
どうやら、nが小さくNが大きいときのほうがオッズが大きく出るようである。
ちょっと理由がすぐには思いつかないので、全体の考察とまとめて午後に考えようと思う。
以上
9月5日（金）午後
ガリガリ計算して確率の対数をとって、もっともらしいpとqの値をそれぞれ0.01～0.05まで0.0050刻みで選べる関数をつくった。
p<-runif(1,0.01,0.05)
q<-runif(1,0.01,0.05)
Surv.time<-function(p,q,l){
cnt<-1
repeat{
if (runif(1)< p) {x<- c(cnt,1);break}
else if (runif(1)< q) {x<- c(cnt,0);break}
else if (cnt==l){x<- c(cnt,0);break}
else cnt<-cnt+1
}
return(x)
}# 生存時間のシミュレーション
A.St<-function(p,q,l){
Surv.time(p,q,l)
}
data.sample<-function(p,q,l,N){
mat<-numeric(0)
for(i in 1:N)
mat<-rbind(mat,A.St(p,q,l))
return(mat)
}# N人を対象にしたデータを作成
d<-data.sample(p,q,20,100)
pt<-function(p,q,t){
log(prod((((1-p)*(1-q))^((1:t)-1))*p))
}# tに死亡が起こる確率
qt<-function(p,q,t){
if(t==20)x<-log((((1-p)*(1-q))^((1:t)-1))*(1-p))
else x<-log(prod((((1-p)*(1-q))^((1:t)-1))*(1-p)*q))
return(x)
}# tに脱落が起こる確率
dp<-function(d){
subset(d,d[,2]==1)
}
dq<-function(d){
subset(d,d[,2]==0)
}
pa<-function(p,q,d){
x<-c()
for(i in 1:length(dp(d)[,1]))x<-c(x,pt(p,q,dp(d)[i,1]))
return(sum(x))
}
qa<-function(p,q,d){
x<-c()
for(i in 1:length(dq(d)[,1]))x<-c(x,qt(p,q,dq(d)[i,1]))
return(sum(x))
}
A<-function(p,q,d){
pa(p,q,dp(d))+qa(p,q,dq(d))
}
mat<-matrix(nrow=9,ncol=9)
ps<-(seq(0.01,0.05,0.005))
qs<-(seq(0.01,0.05,0.005))
for(i in 1:9)
for(j in 1:9)
mat[i,j]<-A(ps[i],qs[j],d)
ex<-function(m){
{if(which(mat==max(mat))%%9==0)x<-ps[9]
else x<-ps[which(mat==max(mat))%%9]}
y<-qs[(which(mat==max(mat))%/%9)+1]
return(c(x,y))
}
p
q
mat
ex(mat)
> p
[1] 0.02417261
>
> q
[1] 0.0497098
>
> mat
[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] -3299.853 -3150.888 -3056.414 -2991.888 -2946.378 -2914.053 -2891.436
[2,] -3252.458 -3103.493 -3009.019 -2944.493 -2898.983 -2866.658 -2844.041
[3,] -3231.308 -3082.343 -2987.869 -2923.343 -2877.832 -2845.508 -2822.890
[4,] -3224.638 -3075.673 -2981.199 -2916.673 -2871.162 -2838.838 -2816.220
[5,] -3227.208 -3078.244 -2983.770 -2919.243 -2873.733 -2841.409 -2818.791
[6,] -3236.227 -3087.262 -2992.788 -2928.261 -2882.751 -2850.427 -2827.809
[7,] -3250.025 -3101.060 -3006.587 -2942.060 -2896.550 -2864.226 -2841.608
[8,] -3267.531 -3118.566 -3024.092 -2959.566 -2914.055 -2881.731 -2859.113
[9,] -3288.010 -3139.045 -3044.571 -2980.045 -2934.534 -2902.210 -2879.592
[,8] [,9]
[1,] -2876.281 -2867.059
[2,] -2828.886 -2819.664
[3,] -2807.736 -2798.514
[4,] -2801.066 -2791.844
[5,] -2803.637 -2794.415
[6,] -2812.655 -2803.433
[7,] -2826.454 -2817.232
[8,] -2843.959 -2834.737
[9,] -2864.438 -2855.216
>
> ex(mat)
[1] 0.025 0.050できるだけfor-loopを使わないようにしていたら、なんだかすごくつくるのに時間がかかった気がする。
処理を早くできたのかはわからないが…。
（前回のコードも書き直そうかと思ったが、for-loopをなくせそうなところがすぐには見つからなかったので、追究するのはやめにした。）
ちなみにこれ、
> exp(sum(mat)) [1] 0
となってしまうので、尤度そのものを出すことは出来なかった。
しかし、尤度の計算の要領はつかめたと思うので、今回はこれでよしとしようと思う。
次は、ケース・コントロール研究のシミュレーションをしていく予定である。
以上
