9月22日（月）
今日は、「ロスマンの疫学」の第10章の残りと、第11章（交互作用についての章）を読んだ。
標準化についてＲでやってみようと思っていたが、基準を決めて重み付き平均を出すだけで、そんなにややこしい話ではないみたいだったので、スルーした。
交互作用は興味深いテーマで、実際の疫学研究でもかなり重要になりそうである。
ただ、詳しい解析は第12章の回帰分析でやるようだ。
交絡についても、回帰分析によるアプローチもあるらしい。
何やら回帰分析が一つの大きなテーマになりそうである。
実際に自分でＲでいろいろ試す題材も豊富そうに思える。
回帰分析については中間発表会後にじっくりやることにして、中間発表会までの残りの時間はいままでやってきたことのおさらいと、木曜の発表会の準備をすることにする。
今日は一応マンテル・ヘンツェルの公式についてもう一度考え直して、Ｒで計算してみたりもしたのだが、結局まとまった結論は出ずに終わった。
次また考え直す機会があれば（最後のまとめのときとか）、参考文献の論文を読んでみようと思う。
以上
9月19日（金）
データの層化 2
kを0から1までの実数として、曝露率をe=0.5+k*xとし、k（すなわち交絡と曝露率の相関の強さ）によって層化の影響がどのように変わるか調べてみた。
N <- 1000 e <- 0.5 b <- 0.05 c.max <- 0.5 alpha <- 0.1 n <- 7 myfunc <- function(k){ deta <- matrix(nrow=N,ncol=3) deta[,1] <- runif(N,0,c.max) f1 <- function(x){ sample(0:1,1,prob=c(0.5-k*x,0.5+k*x)) } deta[,2] <- apply(deta[,1,drop=F],1,f1) f2 <- function(x){ p <- b+x[1]+x[2]*alpha sample(0:1,1,prob=c(1-p,p)) } deta[,3] <- apply(deta[,1:2],1,f2) d <- deta RD.MH <- function(n){ deta.list <- list() for(i in 1:n){ deta. <- subset(d,d[,1]>=(i-1)*c.max/n&d[,1]<i*c.max/n) deta.list[[i]] <- deta. } RD <- function(d){ RD.u <- function(d){ (length(which(d[which(d[,2]==1),3]==1))*length(which(d[,2]==0)) -length(which(d[which(d[,2]==0),3]==1))*length(which(d[,2]==1)))/nrow(d)} RD.l <- function(d){ length(which(d[,2]==0))*length(which(d[,2]==1))/nrow(d) } return(c(RD.u(d),RD.l(d))) } deta2 <- sapply(deta.list,RD) x <- sum(deta2[1,])/sum(deta2[2,]) return(x) } x <- c() for(i in 2:n) x <- c(x,abs(alpha-RD.MH(1))-abs(alpha-RD.MH(i))) return(x) } myfunc2 <- function(k){ mat <- numeric(0) for(i in 1:10000){ mat <- rbind(mat,myfunc(k))} return(mat) } mat0 <- myfunc2(0) mat0.1 <- myfunc2(0.1) mat0.2 <- myfunc2(0.2) mat0.3 <- myfunc2(0.3) mat0.4 <- myfunc2(0.4) mat0.5 <- myfunc2(0.5) mat0.6 <- myfunc2(0.6) mat0.7 <- myfunc2(0.7) mat0.8 <- myfunc2(0.8) mat0.9 <- myfunc2(0.9) x <- 2:n z0 <- apply(mat0,2,mean) z0.1 <- apply(mat0.1,2,mean) z0.2 <- apply(mat0.2,2,mean) z0.3 <- apply(mat0.3,2,mean) z0.4 <- apply(mat0.4,2,mean) z0.5 <- apply(mat0.5,2,mean) z0.6 <- apply(mat0.6,2,mean) z0.7 <- apply(mat0.7,2,mean) z0.8 <- apply(mat0.8,2,mean) z0.9 <- apply(mat0.9,2,mean) plot(x,z0,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=1,xlab="層数",ylab="改善度") par(new=T) plot(x,z0.1,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=2,ann=F) par(new=T) plot(x,z0.2,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=3,ann=F) par(new=T) plot(x,z0.3,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=4,ann=F) par(new=T) plot(x,z0.4,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=5,ann=F) par(new=T) plot(x,z0.5,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=6,ann=F) par(new=T) plot(x,z0.6,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=7,ann=F) par(new=T) plot(x,z0.7,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=8,ann=F) par(new=T) plot(x,z0.8,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),lty=2,ann=F) par(new=T) plot(x,z0.9,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),lty=3,ann=F)
この図のグラフは下から順にk=0,0.1,…,0.9の場合で、kの値が大きくなるほどに、層化による改善度は大きくなることを示している。
しかし、kの値が大きくても、分けるべき層数はせいぜい3くらいで、多くの層に分けなければならないわけではないらしい。
続いて、交絡と曝露率の関係は昨日のシミュレーションのままで、交絡と発症率の関係を変えた場合にどうなるかやってみた。
kを0から1.5までの実数とおき、発症率p=b+k*c(+alpha)とおいた。
N <- 1000
e <- 0.5
b <- 0.05
c.max <- 0.5
alpha <- 0.1
myfunc <- function(k){
deta <- matrix(nrow=N,ncol=3)
deta[,1] <- runif(N,0,c.max)
f1 <- function(x){
sample(0:1,1,prob=c(0.75-x,x+0.25))
}
deta[,2] <- apply(deta[,1,drop=F],1,f1)
f2 <- function(x){
p <- b+k*x[1]+x[2]*alpha
sample(0:1,1,prob=c(1-p,p))
}
deta[,3] <- apply(deta[,1:2],1,f2)
d <- deta
RD.MH <- function(n){
deta.list <- list()
for(i in 1:n){
deta. <- subset(d,d[,1]>=(i-1)*c.max/n&d[,1]<i*c.max/n)
deta.list[[i]] <- deta.
}
RD <- function(d){
RD.u <- function(d){
(length(which(d[which(d[,2]==1),3]==1))*length(which(d[,2]==0))
-length(which(d[which(d[,2]==0),3]==1))*length(which(d[,2]==1)))/nrow(d)}
RD.l <- function(d){
length(which(d[,2]==0))*length(which(d[,2]==1))/nrow(d)
}
return(c(RD.u(d),RD.l(d)))
}
deta2 <- sapply(deta.list,RD)
x <- sum(deta2[1,])/sum(deta2[2,])
return(x)
}
x <- c()
for(i in 2:7)
x <- c(x,abs(alpha-RD.MH(1))-abs(alpha-RD.MH(i)))
return(x)
}
myfunc2 <- function(k){
mat <- numeric(0)
for(i in 1:10000){
mat <- rbind(mat,myfunc(k))}
return(mat)
}
mat0.1 <- myfunc2(0.1)
mat0.2 <- myfunc2(0.2)
mat0.4 <- myfunc2(0.4)
mat0.6 <- myfunc2(0.6)
mat0.8 <- myfunc2(0.8)
mat1 <- myfunc2(1)
mat1.2 <- myfunc2(1.2)
mat1.5 <- myfunc2(1.5)
x <- 2:7
z0.1 <- apply(mat0.1,2,mean)
z0.2 <- apply(mat0.2,2,mean)
z0.4 <- apply(mat0.4,2,mean)
z0.6 <- apply(mat0.6,2,mean)
z0.8 <- apply(mat0.8,2,mean)
z1 <- apply(mat1,2,mean)
z1.2 <- apply(mat1.2,2,mean)
z1.5 <- apply(mat1.5,2,mean)
plot(x,z0,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=1,xlab="層数",ylab="改善度",main="交絡と発症率")
par(new=T)
plot(x,z0.1,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=2,ann=F)
par(new=T)
plot(x,z0.2,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=3,ann=F)
par(new=T)
plot(x,z0.3,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=4,ann=F)
par(new=T)
plot(x,z0.4,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=5,ann=F)
par(new=T)
plot(x,z0.5,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=6,ann=F)
par(new=T)
plot(x,z0.6,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=7,ann=F)
par(new=T)
plot(x,z0.7,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),col=8,ann=F)
par(new=T)
plot(x,z0.8,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),lty=2,ann=F)
par(new=T)
plot(x,z0.9,type="l",xlim=c(2,7),ylim=c(0.00085,0.07),lty=3,ann=F)
下から順に、k=0.2,0.1,0.4,0.6,0.8,1,1.2,1.5の場合である。
やはりこちらも、関係が強いほうが改善度も大きいという結果になった。
また、こちらもkの値に関わらず、層数を上げてもたいして効果が増すわけではないようである。
交絡と発症率との関係が大きいほうが、層数を多くする必要がある、という予想は間違ってはいなかったと思うが、どのみち3つくらいに分ければそれで十分のようである。
特に、全く分けないのと2つに分けるのとで非常に差が大きい。
これは個人的には少し意外な結果だった。（モデルの作り方にもよるのかもしれないが…。）
層別化されたデータを統一化する方法のうち、上記のようにマンテル・ヘンツェルの公式によって計算するのが「プール化」と呼ばれる方法で、実はもう一つ、「標準化」と言われるものがあるらしい。
来週はそれについて考えてみることにして、そこまでを中間発表の内容としたいと思う。
残りの時間と今週末は、図書館で借りてきた教科書（「概説 確率統計」）で、統計学について体系的に勉強しようと思っている。
マイコース後半でさらに高度な内容を扱っていく上で、いまの知識量ではとても太刀打ち出来ないからである。
すでにある程度読み進めてはいるが、やはりちゃんと基礎からやったほうがよく理解できる気がする。
あと、余裕があれば線形代数の行列の知識も多少頭に入れておきたい。
以上
9月18日（木）
データの層化
交絡因子と曝露率に関係があったほうが層化の威力が分かりやすいと思ったので、e=c+0.25として計算してみた。
N <- 1000 e <- 0.5 b <- 0.05 c.max <- 0.5 alpha <- 0.1 n <- 5 myfunc <- function(){ deta <- matrix(nrow=N,ncol=3) deta[,1] <- runif(N,0,c.max) f1 <- function(x){ sample(0:1,1,prob=c(0.75-x,x+0.25)) } deta[,2] <- apply(deta[,1,drop=F],1,f1) f2 <- function(x){ p <- b+x[1]+x[2]*alpha sample(0:1,1,prob=c(1-p,p)) } deta[,3] <- apply(deta[,1:2],1,f2) d <- deta deta.list <- list() for(i in 1:n){ deta. <- subset(d,d[,1]>=(i-1)*c.max/n&d[,1]<i*c.max/n) deta.list[[i]] <- deta. } RD <- function(d){ RD.u <- function(d){ (length(which(d[which(d[,2]==1),3]==1))*length(which(d[,2]==0)) -length(which(d[which(d[,2]==0),3]==1))*length(which(d[,2]==1)))/nrow(d)} RD.l <- function(d){ length(which(d[,2]==0))*length(which(d[,2]==1))/nrow(d) } return(c(RD.u(d),RD.l(d))) } deta2 <- sapply(deta.list,RD) RD.MH <- sum(deta2[1,])/sum(deta2[2,]) return(RD.MH) } x <- c() for(i in 1:10000){ x <- c(x,myfunc())} y <- c() for(i in 1:10000){ y <- c(y,myfunc())} boxplot(x,y,names=c("n=5","n=1"))
このように、層化によって計算結果が真の値に近づくことが確認できた。
ここで、nの値によってどれだけ値が真の値に近づくかを知りたいと思ったので、
同じデータに対して、
（alphaと層数１のときの計算結果との差の絶対値）-（alphaと層数nのときの計算結果との差の絶対値）
を真の値に近づいた量の指標として、nを変えたときの変化を調べた。
N <- 1000 e <- 0.5 b <- 0.05 c.max <- 0.5 alpha <- 0.1 n <- 30 myfunc <- function(){ deta <- matrix(nrow=N,ncol=3) deta[,1] <- runif(N,0,c.max) f1 <- function(x){ sample(0:1,1,prob=c(0.75-x,x+0.25)) } deta[,2] <- apply(deta[,1,drop=F],1,f1) f2 <- function(x){ p <- b+x[1]+x[2]*alpha sample(0:1,1,prob=c(1-p,p)) } deta[,3] <- apply(deta[,1:2],1,f2) d <- deta RD.MH <- function(n){ deta.list <- list() for(i in 1:n){ deta. <- subset(d,d[,1]>=(i-1)*c.max/n&d[,1]<i*c.max/n) deta.list[[i]] <- deta. } RD <- function(d){ RD.u <- function(d){ (length(which(d[which(d[,2]==1),3]==1))*length(which(d[,2]==0)) -length(which(d[which(d[,2]==0),3]==1))*length(which(d[,2]==1)))/nrow(d)} RD.l <- function(d){ length(which(d[,2]==0))*length(which(d[,2]==1))/nrow(d) } return(c(RD.u(d),RD.l(d))) } deta2 <- sapply(deta.list,RD) x <- sum(deta2[1,])/sum(deta2[2,]) return(x) } v <- c(5,10,20,30) x <- c() for(i in v) x <- c(x,abs(alpha-RD.MH(1))-abs(alpha-RD.MH(i))) return(x) } x <- numeric(0) for(i in 1:10000){ x <- rbind(x,myfunc())} boxplot(x[,1],x[,2],x[,3],x[,4],names=c("n=5","n=10","n=20","n=30"),ylab="真の値への近づき")
この結果から、nをかなり大きくしても正確さの程度はほとんど変わらないことが分かる。
さすがにnを30くらいにしたら大きく正確さが損なわれるのかと思っていたので、意外だった。
ただ、今回はnが5のときが若干ではあるがもっともバラつきが小さいようである。
とにかく細かく層化すればよい、というものではないということは確認できた。
明日は、交絡因子と曝露、交絡因子と発症に対しての関係にそれぞれ強弱をつけてみて、層化した場合にどうなるかを検証してみることにする。
あとの時間は、先生に勧められた「統計的因果推論」の第6章の最初のほうを読んでいたが、なかなか難しかった。
ほとんど哲学的とも思えた。
とりあえず、Simpsonのパラドックスは交絡の問題を考えるうえで非常によい題材だと思った。
Simpsonのパラドックス自体は、「因果関係」というものが定義できればあっさり解消されるみたいだが、交絡の影響によってこんなことも起こる、ということは覚えておいて損はないと思う。
第6章の残りもまた少しずつ読むことにする。
以上
9月17日（水）
昨日は、「ロスマンの疫学」の第7章から第10章の途中まで読み進めた。
第9章で、この前のケース・コントロール研究のシミュレーションの結果にも関連しそうな話題も出てきたが、統計学的な理解が追い付かなかったので、一旦保留にした。
とりあえず今日は、第10章の内容である層化データによる交絡の制御についてＲのシミュレーションを行ってみたいと思う。
客観的に観測できる交絡因子（例えば年齢）をcとし、非曝露の場合の発症リスクをb+c、曝露の場合はそこにalphaを加えたb+c+alphaとするモデルを考える。（総人数はN、曝露率はeとする。）
はじめはリスク比を扱おうかと思ったが、コホート研究のリスク比は以前やったので、今回はリスク差で計算してみる。
まずは、このモデルに対し、粗データを用いてリスク差を求めてみた。
N <- 100 e <- 0.5 b <- 0.05 alpha <- 0.1 myfunc <- function(){ c <- runif(N,0,0.5) expo. <- sample(0:1,N,prob=c(1-e,e),replace=TRUE) out. <- c() for(i in 1:N){ p <- b+c[i]+expo.[i]*alpha out. <- c(out.,sample(0:1,1,prob=c(1-p,p)))} deta <- matrix(c(c,expo.,out.),nrow=N,ncol=3) RD.f <- function(d){ length(which(d[which(d[,2]==1),3]==1))/length(which(d[,2]==1))-length(which(d[which(d[,2]==0),3]==1))/length(which(d[,2]==0)) } return(RD.f(deta)) } x <- c() for(i in 1:10000){ x <- c(x,myfunc())} boxplot(x,main="粗データのリスク差（alpha=0.1)")
交絡なしのほうは
c <- rep(0,N)
としている。
このように、平均はだいたい合っているが（曝露・非曝露への交絡因子の交絡が平均化されるため）、交絡があるほうがバラつきが大きくなってしまう。
ここで、層化して計算してみる。
N <- 1000 e <- 0.5 b <- 0.05 c.max <- 0.5 alpha <- 0.1 n <- 5 myfunc <- function(){ c <- runif(N,0,c.max) expo. <- sample(0:1,N,prob=c(1-e,e),replace=TRUE) out. <- c() for(i in 1:N){ p <- b+c[i]+expo.[i]*alpha out. <- c(out.,sample(0:1,1,prob=c(1-p,p)))} deta <- matrix(c(c,expo.,out.),nrow=N,ncol=3) d <- deta deta.list <- list() for(i in 1:n){ deta. <- subset(d,d[,1]>=(i-1)*c.max/n&d[,1]<i*c.max/n) deta.list[[i]] <- deta. } RD.u <- function(d){ (length(which(d[which(d[,2]==1),3]==1))*length(which(d[,2]==0)) -length(which(d[which(d[,2]==0),3]==1))*length(which(d[,2]==1)))/nrow(d)} RD.l <- function(d){ length(which(d[,2]==0))*length(which(d[,2]==1))/nrow(d) } x <- 0 for(i in 1:n){ x <- x+RD.u(deta.list[[i]]) } y <- 0 for(i in 1:n){ y <- y+RD.l(deta.list[[i]]) } RD.MH <- x/y return(RD.MH) } x <- c() for(i in 1:10000){ x <- c(x,myfunc())}
なお、層化した後全体のリスク差を計算するときには、以下に示すマンテル・ヘンツェルの公式というのを用いた。
厳密な導出方法は論文を読まなくては分からないが、なんとなく重み付平均っぽいのが出るのは分かるのでそれでよしとした。
結果としては、ばらつき具合は層化してもほとんど変わらなかった。
とにかく層化さえすれば交絡因子によるばらつきがキャンセルされる…と考えていたが、層化すること自体が結果がばらつく原因になるので（分母が小さくなるから）、必ずしもそんなにうまくはいかないということか（？）。
層化が最も威力を発揮するのは曝露・非曝露の間で交絡因子（ここでは年齢と考えているもの）の差が大きいときなので、明日は曝露率と交絡因子との間に関係を与えて（むしろそのほうが実際に近い）、層化によってどれだけ正確な値に近づくかやってみることにする。
以上
9月12日（金）午後
リード・フロストモデル 拡張版
いまちょうどデング熱が流行っているということで、蚊によって媒介される病気のモデルをつくってシミュレーションしてみた。
（完全に「やってみた」だけなので、だからどうということはないのだが…。）
仮定として、
（1）集団内の人々と蚊は感染・未感染を問わずランダムにまじりあう。
（2）感染者と未感染蚊、感染蚊と未感染者が接触して伝播する一定の確率がそれぞれに存在する。
（3）一度感染した人はもう感染しない。
（4）一度感染した蚊は治癒しない。ただし、一定の割合で死んで新しい未感染蚊と入れ替わる。
（5）その集団は他の集団から隔離されている。
式は
（Ct:時刻tにおける感染者数、St：時刻tにおける未感染者、M.Ct:時刻tにおける感染蚊数、M.St：時刻tにおける未感染蚊数、p：ある1人の感染者が1世代時間にある1匹の未感染蚊に感染を伝播する確率、q：ある1匹の感染蚊が1世代時間にある1人の未感染者に感染を伝播する確率、r:蚊が1世代時間に入れ替わる割合）
とした。
N <- 100 n <- 10 p <- 0.05 q <- 0.05 r <- 0.2 a <- c(1,N-1,n-10,10) l <- 30 mat <-function(t){ x <- matrix(nrow=t,ncol=4) x[1,] <- a for (i in 1:(t-1)){ x[i+1,3] <- (1-r)*x[i,3]+x[i,4]*(1-(1-p)^x[i,1]) x[i+1,1] <- x[i,2]*(1-(1-q)^x[i+1,3]) x[i+1,2] <- N-sum(x[1:i+1,1]) x[i+1,4] <- n-x[i+1,3] } return(x) } x <- 0:l plot(x,c(0,mat(l)[,1]),xlim=c(0,l),ylim=c(0,N),type="l", col=1,xlab="時間",ylab="人数") par(new=T) plot(x,c(N,mat(l)[,2]),xlim=c(0,l),ylim=c(0,N),ann=F,type="l", col=2) par(new=T) plot(x,c(0,N-mat(l)[,1]-mat(l)[,2]),xlim=c(0,l),ylim=c(0,N),ann=F,type="l", col=3) par(new=T) plot(x,c(0,mat(l)[,3]),xlim=c(0,l),ylim=c(0,N),ann=F,type="l", col=4)
全人数を100人、蚊を10匹、p=q=0.05,r=0.2としたときのグラフが下の通り。
他にもいろいろ描いてみて、だいたい予想した通りの結果になったので、満足した。
リード・フロストモデルは非現実的と言われているとはいえ、拡張もしやすいしよくできていると感じた。
来週からはいよいよバイアスとかの話に入る。
一応ここから先が一番興味があった部分だし、少しでも何か成果があればよいと思う。
以上
9月12日（金）午前
リード・フロストモデルによるシミュレーション
一応ネットでも調べておくか…と思って「リード・フロストモデル」と検索したら、自分の昨日の記事が4番目くらいに出てきて唖然とした。
ネット検索において、価値のある情報が上位にくるわけではないというのはたまに言われることだが、それを実感した。
リード・フロストモデルとは、各期間内において、
（1）集団内の人々は感染・未感染を問わずランダムにまじりあう。
（2）感染者と未感染者が接触して伝播する一定の確率が存在する。
（3）一度感染した人はもう感染しない。
（4）その集団は他の集団から隔離されている。
という仮定のもと、
（:時刻tにおける感染者数、
：時刻t+1における感染者数、
：時刻tにおける未感染者、
：ある1人の感染者が1世代時間にある1人の未感染者に感染を伝播する確率）
という式を用いるモデルである。
時間の単位は「世代時間」で、感染してから人に伝播するまでの時間である。
Ct <- function(t,N,n,p){
if(t==0){
return(0)
}else if(t==1){
return(1)
}else{
x <- 0
y <- 1
for(i in 1:(t-1)){
x <- x+Ct(i,N,n,p)
y <- (N-x)*(1-(1-p)^y)
}
return(y)
}}
Rt <- function(t,N,n,p){
if(t==0){
return(0)
}else{
x <- c()
for(i in 0:(t-1)) x <- c(x,Ct(i,N,n,p))
sum(x)
}}
Nt <- function(t,N,n,p){
N-Rt(t,N,n,p)-Ct(t,N,n,p)
}
N <- 100
n <- 1
p <- 0.04
x <- 0:15
y <- c()
for(i in x)
y[i+1] <- Ct(i,N,n,p)
z <- c()
for(i in x)
z[i+1] <- Rt(i,N,n,p)
w <- c()
for(i in x)
w[i+1] <- Nt(i,N,n,p)
plot(x,y,xlim=c(0,15),ylim=c(0,100),xlab="時間",ylab="人数",main="p=0.04のとき",type="l",lty=1)
par(new=T)
plot(x,z,xlim=c(0,15),ylim=c(0,100),ann=F,type="l", lty=2)
par(new=T)
plot(x,w,xlim=c(0,15),ylim=c(0,100),ann=F,type="l", lty=3)
text(8,90,"免疫獲得者")
text(1.2,80,"感受性ある者")
text(7,20,"感染者")
一応、「ロスマンの疫学」に書いてあった通りのシミュレーションはできた。
せっかくなので午後は自分なりのアレンジを加えてちょっとやってみる。
以上
9月11日（木）
累積症例対照研究のシミュレーション まとめ 2
結論から言うと、どうして分散がこのように変化するのかの説明は出来なかった。
ある多項分布にしたがう4つの確率変数のかけ算・割り算をしたときに、計算結果はどういった分布を示すか、という問題に帰着できて、理論上はちゃんとそれぞれのパラメータの変化による分散の説明もできるような気がするのだが、ちょっと僕の手には負えなかった。
それに、おとといの記事で
「いま、オッズ比の計算は、
（症例群の曝露者数）/（症例群の非曝露者数）を（対照群の曝露者数）/（対照群の非曝露者数）で割って出しているが、
真の値が（症例群の曝露者数）＝（症例群の非曝露者数）かつ（対照群の曝露者数）＝（対照群の非曝露者数）に近い値のほうが、より収束するのが早いと考えられる。」
と書いたが、それは確かにそうでも、（確率変数）/（確率変数）の真の平均値が小さければ小さいほど（すなわち、分子の平均値が小さくて分母の平均値が大きいほど）、分散の値も小さくなりやすいはずである。
そのため、いまのままではe1=0.1のときのほうが分散の値が大きいことの説明が十分ではない。
ちょっと解決の目途がいまのところたたない上、かなり時間を消費してしまっているので、この問題は一旦保留にしようと思う。
保留にしている問題が増えてきたので、マイコース終盤のまとめの段階で保留にしている問題について再考する日をつくることにする。
解決には到らなかったが、それに関連する話題で、いろいろ寄り道したりしつつ統計に関する理解を深めることが出来たとは思うので、一応考えて意味はあったかなと思う。
あと、一般教養で取り扱った数学が思ったよりも多くのところで応用されていることを知り、1回生のときにもっと真面目に勉強しておけばよかったと後悔した。
勉強に割いた時間が長かったので、今日はほとんど記事にすべき作業はなかった。
明日は第6章の「感染症疫学」に進み、遊び半分でリード・フロストモデルおよび自分のオリジナルのシミュレーションをＲでやってみることにする。
その次の章からバイアスや交絡といった範囲に入っていくので、そのへんの勉強がひととおり済んだら論文を読んでみようと思う。
以上
