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章の残りもまた少しずつ読むことにする。

以上