Stanで疑似相関を見抜く - 統計コンサルの議事メモ

ベイズ統計をちゃんと勉強しようと思い、最近Statistical Rethinkingという本を読んでいます。

Statistical Rethinking: A Bayesian Course with Examples in R and Stan (Chapman & Hall/CRC Texts in Statistical Science)

作者: Richard McElreath
出版社/メーカー: Chapman and Hall/CRC
発売日: 2016/02/19
メディア: ハードカバー
この商品を含むブログを見る

その中でタイトルの通り面白い話があったので紹介がてら実行してみます。以下は５章１節からの引用ですが、x_spurとx_realはそれぞれyと相関の「ない」変数と相関の「ある」変数を表しています：

But when you include both x variables in a linear regression predicting y, the posterior mean for the association between y and x_spur will be close to zero, while the comparable mean for x_real will be closer to 1.

この記事の内容を端的に言えば「yと相関のあるx1と、x1と相関のあるx2（つまりyとは直接相関がない）が説明変数として与えられている時、stanを用いることで不要な変数（この場合x2）を除外できる」というものになります。
またこちらの本では全体を通して{rethinking}パッケージによってパラメータの推定が行われていますが、ここではせっかくなので{rstan}でやってみます。

前準備

まずは{rstan}を読み込みます。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

続いてy、x1、x2をそれぞれ以下のように定義します。

## 乱数のシードの設定
set.seed(15)
N  <- 100 # サンプルサイズ
x1 <- rnorm(N) # yと相関のある（真の）説明変数
x2 <- rnorm(N, x_real) # x1と相関することでyと相関する（偽の）説明変数
y  <- rnorm(N, x_real) # x1と相関する目的変数

## データフレームにする
dat <- data.frame(y, x_real, x_spur)

定義からx2はyと相関しないため、本来であれば統計モデルから除外されるべきものです。しかし通常はデータを取得した段階で上記のような情報は得られないため、モデリングした結果を見て判断することになります。ちょっとやってみましょう。

## 回帰分析の実行（一部省略）
> summary(res.lm <- lm(y ~ ., data = dat))

Call:
lm(formula = y ~ ., data = dat)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.1262 -0.5336  0.0269  0.6975  1.8557 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.10912    0.10485  -1.041    0.301    
x_real       1.27236    0.14407   8.832 4.45e-14 ***
x_spur      -0.19652    0.09488  -2.071    0.041 *

このデータではx1、x2ともにyに対して有意に（P < 0.05）影響していることとなり*1、本来であれば不要であるx2をモデルから落とすことができません。当然かもしれませんが、両者ともに有意であるためstep関数などを用いた変数選択も意味がありません。

## stepによる変数選択の実行
> step(res.lm)

Start:  AIC=10.25
y ~ x_real + x_spur

         Df Sum of Sq    RSS    AIC
<none>                104.34 10.248
- x_spur  1     4.615 108.95 12.576
- x_real  1    83.901 188.24 67.255

Call:
lm(formula = y ~ x_real + x_spur, data = dat)

Coefficients:
(Intercept)       x_real       x_spur  
    -0.1091       1.2724      -0.1965

Stanによる推定

それでは上記のデータセットを用いてstanで推定を実行した場合はどのようになるか試してみましょう。まずはstanに渡すデータを定義します。

## stanのdataブロックに渡すデータの定義
dat_stan <- list(
   N = nrow(dat),
   D = 2,
   X = dat[, 2:3],
   Y = dat[, 1]
)

stanのスクリプトは、以下のような内容になります：

data {
   int<lower=0> N;
   int<lower=0> D;
   matrix[N, D] X;
   vector[N] Y;
}

parameters {
   real a;
   vector[D] beta;
   real<lower=0> sig;
}

model {
   for (d in 1:D) {
      beta[d] ~ normal(0, 1);
   }
   Y ~ normal(a + X * beta, sig);
}

特に工夫した点などもなく、単に重回帰を推定するためのスクリプトとなっていますが、このスクリプトを用いてstan()を実行してみましょう*2。

## フィッティング
fit_01 <- stan(file = './StanModel/Detect_Confound_Factor.stan', 
               data = dat_stan, 
               iter = 3000,
               chains = 4,
               seed = 1234)

結果はこちらとなります*3

> summary(fit_01)$summary[, c("mean", "2.5%", "50%", "97.5%", "Rhat")]

               mean        2.5%         50%         97.5%      Rhat
a        -0.1055297  -0.3129441  -0.1052853   0.102064821 0.9998115
beta[1]   1.2470319   0.9587653   1.2454155   1.530978156 1.0001602
beta[2]  -0.1827657  -0.3682030  -0.1828429   0.003176939 1.0005802
sig       1.0497539   0.9158756   1.0444427   1.210491329 0.9995408
lp__    -54.9435610 -58.5631720 -54.6322411 -53.160356164 1.0010219