2010-03-24

メタアナリシスの等高線

メタアナリシスは複数のスタディの結果を統合する手法です。統合して、帰無仮説が棄却されるかどうかの判断と、採択するとして、「効果の大きさ」の推定値とその信頼区間を算出します。
２つの２ｘ２表でその基本的な雰囲気を感じてみます。
まず、２個のスタディの個々の検定結果が登場し
その２個の表を、「単純に足し合わせ」たときに、統計量の等高線がどうなるかを見て
ついで、固定効果モデルの代表としてマンテル-ヘンツェル法を、変量効果モデルの代表としてDerSimonian-Laired法の等高線を見て
それぞれの特徴を了解することにします。
掲載図は、横軸がスタディ１のセルの期待値からのずれ(スタディ１のカイ自乗値はこの値の二乗に比例)、縦軸がスタディ２のセルの期待値からのずれ(スタディ２のカイ自乗値はこの値の二乗に比例)、黒・青・赤の等高線がそれぞれ、「表の単純和(自由度２)」「マンテル-ヘンツェル(自由度１)」「DerSimonian-Laired(自由度１)」

t1<-matrix(c(10,10,10,10),nrow=2,byrow=TRUE)
#t2<-matrix(c(20,30,50,40),nrow=2,byrow=TRUE)
t2<-matrix(c(50,100,150,200),nrow=2,byrow=TRUE)
t2<-t1
m11<-apply(t1,1,sum)
m12<-apply(t1,2,sum)
M1<-sum(t1)

m21<-apply(t2,1,sum)
m22<-apply(t2,2,sum)
M2<-sum(t2)

e1<-m11%*%t(m12)/M1
e2<-m21%*%t(m22)/M2

x11<-seq(from=-M1,to=M1,by=1)

y11<-seq(from=-M2,to=M2,by=1)

x12<--x11+m11[1]
x21<--x11+m12[1]
x22<--x12+m12[2]

xbind<-cbind(x11,x12,x21,x22)
okx<-which(apply(xbind,1,min)>0)

x11<-x11[okx]
x12<-x12[okx]
x21<-x21[okx]
x22<-x22[okx]

y12<--y11+m21[1]
y21<--y11+m22[1]
y22<--y12+m22[2]

ybind<-cbind(y11,y12,y21,y22)
oky<-which(apply(ybind,1,min)>0)

y11<-y11[oky]
y12<-y12[oky]
y21<-y21[oky]
y22<-y22[oky]


chi1<-(x11-e1[1,1])^2/e1[1,1]+(x12-e1[1,2])^2/e1[1,2]+(x21-e1[2,1])^2/e1[2,1]+(x22-e1[2,2])^2/e1[2,2]

chi2<-(y11-e2[1,1])^2/e2[1,1]+(y12-e2[1,2])^2/e2[1,2]+(y21-e2[2,1])^2/e2[2,1]+(y22-e2[2,2])^2/e2[2,2]



z<-outer(chi1,chi2,FUN="+")

#filled.contour(x11-e1[1,1],y11-e2[1,1],z)
xlim<-ylim<-c(min(x11-e1[1,1],y11-e2[1,1]),max(x11-e1[1,1],y11-e2[1,1]))
#contour(x11-e1[1,1],y11-e2[1,1],sumz,xlim=xlim,ylim=ylim)
contour(x11-e1[1,1],y11-e2[1,1],z,xlim=xlim,ylim=ylim)
abline(h=0)
abline(v=0)

sum11<-outer(x11,y11,FUN="+")
sum12<-outer(x12,y12,FUN="+")
sum21<-outer(x21,y21,FUN="+")
sum22<-outer(x22,y22,FUN="+")

sume11<-e1[1,1]+e2[1,1]
sume12<-e1[1,2]+e2[1,2]
sume21<-e1[2,1]+e2[2,1]
sume22<-e1[2,2]+e2[2,2]

sumz<-(sum11-sume11)^2/sume11+(sum12-sume12)^2/sume12+(sum21-sume21)^2/sume21+(sum22-sume22)^2/sume22

xlim<-ylim<-c(min(x11-e1[1,1],y11-e2[1,1]),max(x11-e1[1,1],y11-e2[1,1]))
contour(x11-e1[1,1],y11-e2[1,1],sumz,xlim=xlim,ylim=ylim)
par(new=T)
contour(x11-e1[1,1],y11-e2[1,1],z,col="red",xlim=xlim,ylim=ylim)
abline(a=0,b=1)

library(rmeta)
n.case<-c(m11[1],m21[1])
n.ctrl<-c(m11[2],m21[2])
zMH<-z
zDSL<-z
for(i in 1:length(z[,1])){
 for(j in 1:length(z[1,])){
  mhout<-meta.MH(n.case,n.ctrl,c(x11[i],y11[j]),c(m12[1]-x11[i],m22[1]-y11[j]))
  zMH[i,j]<-mhout$MHtest[1]
  dslout<-meta.DSL(n.case,n.ctrl,c(x11[i],y11[j]),c(m12[1]-x11[i],m22[1]-y11[j]))
  zDSL[i,j]<-dslout$test[1]^2
 }
}
zlim<-c(0,max(sumz))
contour(x11-e1[1,1],y11-e2[1,1],sumz,xlim=xlim,ylim=ylim,zlim=zlim,nlevels=10)
par(new=T)
contour(x11-e1[1,1],y11-e2[1,1],zMH,col="red",xlim=xlim,ylim=ylim,zlim=zlim,nlevels=10)

contour(x11-e1[1,1],y11-e2[1,1],z,xlim=xlim,ylim=ylim,zlim=zlim,,nlevels=10)
par(new=T)
contour(x11-e1[1,1],y11-e2[1,1],zDSL,col="red",xlim=xlim,ylim=ylim,zlim=zlim,nlevels=10)