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つに分けるのとで非常に差が大きい。
これは個人的には少し意外な結果だった。（モデルの作り方にもよるのかもしれないが…。）
層別化されたデータを統一化する方法のうち、上記のようにマンテル・ヘンツェルの公式によって計算するのが「プール化」と呼ばれる方法で、実はもう一つ、「標準化」と言われるものがあるらしい。
来週はそれについて考えてみることにして、そこまでを中間発表の内容としたいと思う。
残りの時間と今週末は、図書館で借りてきた教科書（「概説 確率統計」）で、統計学について体系的に勉強しようと思っている。
マイコース後半でさらに高度な内容を扱っていく上で、いまの知識量ではとても太刀打ち出来ないからである。
すでにある程度読み進めてはいるが、やはりちゃんと基礎からやったほうがよく理解できる気がする。
あと、余裕があれば線形代数の行列の知識も多少頭に入れておきたい。
以上
