宮川(2004)「統計的因果推論」(朝倉書店)の事例を使い、バックドア基準の効果を確かめてみる。
いま、p.107の図6.2左をもとに真のモデルを作り、X→Yの因果効果の強さを推定することを試みる
##まずは、バックドア 基準を使わず、素朴にyに対するxの線形モデルを作って最小二乗推定する場合
##ここで、真のモデルにおけるX→Yの強さをalpha, 各ノードのサンプルサイズをnとしてデータ生成→最小二乗推定する作業を1万回反復し、1万個の推定値をreturnさせる
non_back_sim<-function(alpha,n){
tmp<-numeric(10000)
for (i in 1:10000){
z2<-rnorm(n,0,10)
z1<-1.2*z2+rnorm(n,0,4)
w<--0.5*z1+rnorm(n,0,2)
x<-0.5*z2+rnorm(n,0,5)
y<-z1+w+alpha*x+rnorm(n,5)
tmp[i]<-solve(t(x)%*%x)%*%t(x)%*%y
}
return(tmp)}
##次に、バックドア基準を満たすz2を回帰モデルに入れた場合
>back_sim2<-function(alpha,n){
tmp<-numeric(10000)
for (i in 1:10000){
x<-matrix(0,n,2)
x[,2]<-rnorm(n,0,10)
z1<-1.2*x[,2]+rnorm(n,0,10)
w<--0.5*x[,2]+rnorm(n,0,2)
x[,1]<-0.5*x[,2]+rnorm(n,0,5)
y<-z1+w+alpha*x[,1]+rnorm(n,5)
tmp[i]<-solve(t(x)%*%x)%*%t(x)%*%y
}
return(tmp)}
##最後に、z1を使う場合
back_sim<-function(alpha,n){
tmp<-numeric(10000)
for (i in 1:10000){
x<-matrix(0,n,2)
z2<-rnorm(n,0,10)
x[,2]<-1.2*z2+rnorm(n,0,10)
w<--0.5*x[,2]+rnorm(n,0,2)
x[,1]<-0.5*z2+rnorm(n,0,5)
y<-x[,2]+w+alpha*x[,1]+rnorm(n,5)
tmp[i]<-solve(t(x)%*%x)%*%t(x)%*%y
}
return(tmp)}
ヒストグラムやカーネル密度推定を使うと、最初のシミュレーションでは不偏性を満たしていないことがわかる。
一方バックドア基準を満たす二つを比較すると、ともに真の値の周りに推定値がばらつくものの、その分散はz1を使った方が小さくなることがわかる。