9月17日（水） - 2014年度マイコースプログラム

昨日は、「ロスマンの疫学」の第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)")

f:id:kouri_don:20140917151930p:plain

交絡なしのほうは

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

f:id:kouri_don:20140917170407p:plain

なお、層化した後全体のリスク差を計算するときには、以下に示すマンテル・ヘンツェルの公式というのを用いた。

${RD_{MH}=\frac{\displaystyle \sum_{i}^{}\frac{\displaystyle a_{i}N_{0i}-b_{i}N_{1i}}{\displaystyle T_i}}{\displaystyle \sum_{i}^{}\frac{\displaystyle N_{1i}N_{0i}}{\displaystyle T_i}}}$

厳密な導出方法は論文を読まなくては分からないが、なんとなく重み付平均っぽいのが出るのは分かるのでそれでよしとした。

結果としては、ばらつき具合は層化してもほとんど変わらなかった。
とにかく層化さえすれば交絡因子によるばらつきがキャンセルされる…と考えていたが、層化すること自体が結果がばらつく原因になるので（分母が小さくなるから）、必ずしもそんなにうまくはいかないということか（？）。

層化が最も威力を発揮するのは曝露・非曝露の間で交絡因子（ここでは年齢と考えているもの）の差が大きいときなので、明日は曝露率と交絡因子との間に関係を与えて（むしろそのほうが実際に近い）、層化によってどれだけ正確な値に近づくかやってみることにする。

以上