餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

ベイズ推定の事前分布で多重共線性に対処しよう

線形回帰などでベイズ統計学を用いる利点は何なのかと言う御題が流れていました。ざっと思いつくところで、

  • 主観的事前分布で、データ以外の情報を加味できる
  • (同じことだが)逐次ベイズ推定にできる
  • 信頼区間よりも信用区間の方が解釈がしやすい
  • 有意水準のギリギリ前後で解釈を変えなくてよい
  • モデル選択でベイズファクターが使える

と言ったことがあげられますが、「データ以外の情報」のありそうな例を用意していなかったことを思い出しました。先行するデータ分析結果や、理論から演繹される情報を入れるのが一般的なようですが、もっと素朴に分析者が係数の符号の向きを知っている状態を考えてみましょう。

理屈から効果の正負の方向に予想がついている一方で、サンプルサイズが小さいために多重共線性で、推定結果の符号が予想と逆になることは多々あります。例えば、年収と(学部までの)学歴で既婚率を説明したときに、年収はプラスの効果があるのに、学歴はマイナスの効果が計測されたら、年収と学歴の多重共線性を疑うべきでしょう。

y = ꞵ₀ + ꞵ₁x₁ + ꞵ₂x₂ + ϵが真のモデルで、x₁とx₂で多重共線していて、かつ理論的にꞵ₁>0でꞵ₂<0と言う制約がある場合の推定をしてみましょう。

1. データセット

真のモデルを設定して、多重共線性のあるデータセットを作ります。

set.seed(105)
n <- 30
beta <- c(1, 2, -3)
names(beta) <- c("(Intercept)", "x1", "x2")
x1 <- runif(n, min=0, max=1)
x2 <- x1 + rnorm(n, mean=0, sd=0.2)
y <- beta[1] + beta[2]*x1 + beta[3]*x2 + rnorm(n, mean=0, sd=1)
cat(sprintf("cor(x1, x2)=%.3f\n", cor(x1, x2)))

もちろん、OLSが上手くいかないシード値を探した結果です。

2. OLSで推定

一般最小二乗法(OLS)で推定すると、x₁の係数が真のモデルと逆になっています。

r_lm <- lm(y ~ x1 + x2)
summary(r_lm)

なお、VIFは4.23と大きくは無いのですが、LOOCVをかけるとx2だけで回帰した方がよくなります*1

3. 事前分布に正規分布を置いてベイズ推定

誤差項を正規分布とし、x1とx2の係数の事前分布を平均aと-a, 標準偏差a/2の正規分布としてベイズ推定します。つまり、予想の逆の符号になるのは平均よりも2σ以上外側と言う主観を置きました。

library(MCMCpack)

a <- 5
post.nd <- MCMCregress(y ~ x1 + x2, b0=c(0, a, -a), B0=1/(a/2)^2)
summary(post.nd)

係数の事前分布に正規分布を置く線形回帰は標準的な関数なので簡単に推定できます。しかし、符号条件を満たす推定結果が得られますが、aを幾つにするか恣意的に思えるかも知れません。実際、aを変えると推定結果は大きく変化します。

4. 事前分布に一様分布を置いてベイズ推定

恐らく[0 a]と[-a 0]の一様分布を事前分布に置く方が説明は簡単になります。aは一定以上あれば、推定結果を左右しません。

library(MCMCpack)

llf_nd <- function(theta){
# concentrated loglikelihood functionにして分散の推定は省略
  e <- y - theta[1] - theta[2]*x1 - theta[3]*x2
  sd <- sqrt(sum((e - mean(e))^2)/length(e))
  sum(dnorm(e, mean=0, sd=sd, log=TRUE))
}

prior.u <- function(theta){
  a <- 100
  sum(
    dunif(theta[1], min=-a, max=a, log=TRUE)
    ,dunif(theta[2], min=0, max=a, log=TRUE)
    ,dunif(theta[3], min=-a, max=0, log=TRUE))
}

post.u <- MCMCmetrop1R(function(theta){
  v <- llf_nd(theta) + prior.u(theta)
  # -Infを返すとエラーが出るので絶対値が大きな負の数を戻す(大きすぎても特異になってエラー)
  if(!is.finite(v)){
    return(-1e+8)
  }
  return(v)
}, theta.init=c(0, 1, -1), optim.method="L-BFGS-B", optim.lower=c(-Inf, 1e-12, -Inf), optim.upper=c(Inf, Inf, -1e-12))
summary(post.u)


5. 比較

推定された係数を比較してみましょう。

cmp_coef <- rbind(beta, coef(r_lm), summary(post.nd)$statistics[1:3, 1], summary(post.u)$statistics[, 1])
rownames(cmp_coef) <- c("True Model", "OLS", "Bayes(gauss)", "Bayes(unif)")
cmp_coef

分析者の知識を事前確率として織り込むことで、OLSよりもベイズ推定の方が上手くいっていることが分かります*2。ただし主観的事前分布をどう置くかはやはり問題になります。正規分布を使う方がコードの量は減らせますが、裁量の余地は増えます。一様分布を置くと、裁量の余地は減るわけですが、不連続な点が出てくるためにコードの量が増えてしまいます。悩ましいですね。

*1:VIFの計算とLOOCVの具体的な方法は、「多重共線性を出してみよう」を参照してください

*2:Ridge回帰を使えば良い気がしますが、パラメーターの信用区間が出ることが利点になります。なお、主観的事前分布が受け入れ難い人には、主成分回帰でごちゃごちゃやる手があります。