餡子付゛録゛

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

Rで線形混合モデル×多重代入法

反復測定分散分析をlme4パッケージとmiceパッケージを使った線形混合モデル*1×多重代入法で行なう作業を確認したので公開します*2。以前のエントリー「反復測定分散分析を主観ベイズ推定に置き換えよう」で、こっちの方が楽と書いたのですが、欠損値補定がやや煩雑だったので。

1. データセット

前回とほとんど同じモノにしますが、比較のために欠損値なしのデータセットも作っておきます。細かい説明は以前のエントリーを参照してください。

set.seed(604)
rm(list=ls())

n <- 120 # 各時点のサンプルサイズ
t <- 1:3 # 時点
g <- c("C", "T") # 群

dimnames <- list(g, sprintf("t=%d", t))

# 以下のmu, sigma, betaは真の値になるので、本当は分からない
mu <- matrix(c(0+0*t, 0+1*(t-1)), ncol=length(t), byrow=TRUE, dimnames=dimnames)
sigma <- matrix(c(1+0*t, 1+0.1*(t-1)), ncol=length(t), byrow=TRUE, dimnames=dimnames)
beta <- 1

x <- rep(runif(n), length(t))

df01 <- data.frame(
    id=rep(1:n, length(t)),
    time=rep(t, each=n),
    group=rep(rep(g, each=n/length(g)), length(t)),
    x=x,
    z=x/2 + rnorm(n, sd=1),
    score=numeric(n*length(t))
)

rm(x)

for(j in t){
    for(i in 1:length(g)){
        s <- n*(j-1)+n/length(g)*(i-1) + 1
        e <- n*(j-1)+n/length(g)*i
        df01$score[s:e] <- rnorm(n/length(g), mean=mu[i, j], sd=sigma[i, j]) + beta*df01$x[s:e]
    }
}

# xに欠損値をつくる
df00 <- df01
df01$x[df01$id %in% sample(1:n, n/2)] <- NA

2. 補定と推定

miceパッケージのお手軽関数は、lme4パッケージが出力するオブジェクトを現時点で処理できません。miceが補定*3したデータセットに逐次推定をまわし、その結果となる係数と分散共分散行列の推定量をリストにまとめて、miceaddsパッケージ、mitoolsパッケージに渡して推定結果を統合します。こう書くと煩雑で大変そうですが、偉い人が情報をまとめておいてくれました*4

# 利用パッケージを呼び出す
library(lme4)
library(mice)
library(miceadds)
library(mitools)

# miceで補定したm個のデータセットをつくる
# 欠損値があるのは変数xのみ
mice.out <- mice(df01, blocks=c("x"), print=FALSE, m=100)

# 補定したm個のデータセットをそれぞれ推定する
mods <- lapply(1:mice.out$m, function(i){
    # m個の補定データセットの中からi番目を取り出す
    dfm <- complete(mice.out, action=i)
    # i番目で推定をかけて結果を戻すとリストに格納される
    lmer(score ~ x + factor(time)*factor(group) + (1|group), data = dfm)
})

# 推定結果の集合から係数を抜き出す
cmod <- MIextract(mods, fun=function(r_lmer){
    # lme4の癖でcoef(r_lmer)が複数行返ってくるため
    summary(r_lmer)$coefficients[,1]
})

# 推定結果の集合から分散共分散行列を抜き出す
vmod <- MIextract(mods, fun=function(r_lmer){
    # lme4の癖でvcov(r_lmer)が行列にならないため
    as.matrix(vcov(r_lmer))
})

# 自由度を計算
df <- with(summary(mods[[1]]), length(residuals) - length(coefficients[,1]))

# 推定結果を統合
r_pool <- pool_mi(qhat=cmod, u=vmod, dfcom=df)

# coef(r_pool)とvcov(r_pool)ができるので、後はまとめる

欠損値補定後にクロス項をつくっていますが、今回はクロス項に使う変数は補定されていないので問題ないです。unable to evaluate scaled gradientや Model failed to converge: degenerate Hessian with 1 negative eigenvaluesと言う警告が出るのですが、どうも内部で生成している乱数が原因で、推定結果を大きく左右するものでは無さそうなので無視します。

3. 結果の比較

欠損値補定、劇的に推定結果が代わるというものではないのですが、欠損値なしの場合と、listwise法の場合と、miceの場合で比較してみましょう。

# 欠損値なしデータでの推定(比較用)
mixed.lmer.c <- lmer(score ~ x + factor(time)*factor(group) + (1|group), data = df00)

# listwise法(比較用)
mixed.lmer.m <- lmer(score ~ x + factor(time)*factor(group) + (1|group), data = df01)

cmp_beta <- cbind(summary(mixed.lmer.c)$coefficients[,1], summary(mixed.lmer.m)$coefficients[,1], coef(r_pool))
cmp_se <- cbind(summary(mixed.lmer.c)$coefficients[,2], summary(mixed.lmer.m)$coefficients[,2], sqrt(diag(vcov(r_pool))))
colnames(cmp_beta) <- colnames(cmp_se) <- c("No Missing", "listwise", "MI/mice")

listwise法よりも欠損値なしに近づけば改善ですが、係数の比較は改善しているか、悪化しているのか良く分かりません。

標準誤差は改善していると言えるでしょう。

A. Ameliaを使う

lme4パッケージとmiceパッケージを併用するにはどうしたらいいかと言う質問に、Ameliaパッケージを使えと回答がついていたのを見かけて試してみたのですが、微妙でした。Amelia向きのデータ生成のはずですが。

library(Amelia)
library(merTools)
library(plyr)

# Ameliaで補定したm個のデータセットをつくる
amelia.out <- amelia(df01, idvars="id", ts="time", cs="group", m=100)

# 補定したm個のデータセットをそれぞれ推定する
mods <- lapply(amelia.out$imputations, function(dfa){
    lmer(score ~ x + factor(time)*factor(group) + (1|group), data=dfa)
})

# 推定結果の集合から係数を抜き出す
beta <- sapply(mods, function(r_lmer){
    summary(r_lmer)$coefficients[, 1]
})

# 推定結果の集合から標準誤差を抜き出す
se <- sapply(mods, function(r_lmer){
    summary(r_lmer)$coefficients[, 2]
})

# 合成する(分散共分散行列でなくて良いのか感があるが、良いらしい)
r_mi <- mi.meld(q=beta, se=se, byrow=FALSE)

# 比較表に追加
cmp_beta <- cbind(cmp_beta, c(r_mi$q.mi))
cmp_se <- cbind(cmp_se, c(r_mi$se.mi))
colnames(cmp_se)[4] <- colnames(cmp_beta)[4] <- "Amelia"

*1:線形混合モデルとその応用例 - Qiita

*2:データセットの都合で、一般化線形混合モデルではなくて、線形混合モデルになっています。lme4パッケージの場合、lmerをglmerにすれば一般化線形混合モデルになるので、だいたい同じノリで処理できますが。

*3:階層構造を指定したりと利用がややこしいので、今回は使いませんでしたが、miceaddsパッケージのmice.impute.ml.lmerを使うと、線形混合モデルを使った補定が可能になるようです。

*4:[R][初心者の質問] mice で多重代入法の結果を統合できない時の対処法 - ill-identified diary

OpenBLASをリンクしたWindows版R 4.2 PR

標準のNetlib BLASよりも高速な線形代数ライブラリOpenBLASをリンクしたR 4.1.xを使ってきたのですが、今年から4.2系が標準になるようです。古いRを使い続けていると、更新されたパッケージを使うときに警告で出て目障りなので、バージョンアップをしてみました。ちょっとだけチートペーパーと手順を変えたので、記録しておきます。

1. 補助ツールのインストール

1.1. Inno Setupのインストール

公式サイトからインストーラーinnosetup-6.2.1.exeをダウンロードしてきて、Inno Setupをインストール*1。インストール先はデフォルトで*2

1.2. MikTeXのインストール

公式サイトからbasic-miktex-22.3-x64.exeをダウンロードしてきて、MikTeXをインストール。付属ユーティリティMiKTeX Consoleで、administration modeを選択し、Check for updates、Update now、メニューのTasksのRefresh file name database、Refresh font map files、Update package databaseを順番に実行。

1.3. QPDFのインストール

SourceForgeのQPDFのリポジトリからバイナリqpdf-10.6.3-bin-mingw32.zipをダウンロードしてきて、C:で展開。

2. Rtools42のインストールから、R 4.2PRのコンパイルまで

Rの公式サイトのRtools42のフォルダーからRtools42 installer(rtools42-5253-5107.exe)をダウンロードしてきてインストールします。
まず、Rtools42 Bashを起動して、wgetコマンドを追加して、Rtoolsのアップデートをします。

pacman -Sy wget
pacman -Syuu

自動終了するので、Rtools42 Bashを再び起動して、(一時的に使う)PATHを設定します

export PATH=/c/rtools42/x86_64-w64-mingw32.static.posix/bin:/c/rtools42/usr/bin:$PATH
export PATH=/c/Program\ Files/MiKTeX/miktex/bin/x64:$PATH
export TAR="/usr/bin/tar"
export TAR_OPTIONS="--force-local"

Rtools42 Bashを終了したら、再設定になるので注意してください。
リリース候補版の最新ソースコードをダウンロードしてきて、C:に置きます*3

cd /c
wget https://cran.r-project.org/src/base-prerelease/R-latest.tar.gz
tar zxf R-latest.tar.gz --force-local
# 「シンボリックリンクが作れません」とエラーが出た場合は、展開したファイルを残したまま、同じコマンドで展開すると誤魔化せる*4

C:\R-rcが出来ます*5

Rの公式サイトのRtools42のフォルダーにあるTcl/Tkのソースコードをダウンロードしてきて、C:\R-rcに展開します。

export wdir=/c/R-rc # 展開先がC:\R-rcでない場合は、ここを修正
cd $wdir
wget https://cran.r-project.org/bin/windows/Rtools/rtools42/files/tcltk-5253-5175.zip
unzip tcltk-5253-5175.zip

なお、ファイル名tcltk-5253-5175.zipは、今後、変わっていく可能性があるので、適時変更してください。
そして、src/gnuwin32に移動して、MkRules.localをつくります。qpdfのフォルダ名に注意してください。

cd $wdir/src/gnuwin32
cat <<EOF >MkRules.local
USE_ATLAS = YES
EOPTS = -march=native -pipe -mno-rtm -mno-fma
LTO = -flto -ffat-lto-objects -fuse-linker-plugin
LTO_OPT = -flto -ffat-lto-objects -fuse-linker-plugin
LTO_FC_OPT = -flto -ffat-lto-objects -fuse-linker-plugin
QPDF = C:/qpdf-10.6.3
OPENMP = -fopenmp
EOF

QPDFのパスはインストール先にあわせてください。
src/extra/blas/Makefile.winをsedで置換します。notepadで編集したいのですが、UNIX改行コードを認識しないので。

cd $wdir/src/extra/blas
mv Makefile.win Makefile.win.old
sed 's/-L"$(ATLAS_PATH)" -lf77blas -latlas/-fopenmp -lopenblas/' < Makefile.win.old > Makefile.win

Makefile.winが変わっているとパターンマッチしないので、cat Makefile.winをして、しっかり書き換わっているかチェックしてください。

PATHが通っているか確認します。

which make gcc pdflatex tar

エラーが出なければ問題なしです。
コンパイルを開始します。

cd $wdir/src/gnuwin32
make distribution

無事、終われば C:\R-rc\src\gnuwin32\installerにR-4.2.1rc-win.exe*6が出来上がっています。私はMikTeXのアップデート後の処理が抜けて!pdfTeX error: pdflatex.exe (file ts1-zi4r): Font ts1-zi4r at 540 not foundと言われたり、QPDFのパス指定を誤ったりして、修正後、make distributionをやり直しました。

3. 動作確認

Windows版ではsessionInfo()をしても使っているBLASの種類を教えてくれないのですが、マルチコアの計算機で以下のように行列演算をさせて、ユーザー時間が経過時間よりも大きければ並行処理ができているので、OpenBLASとリンクしているのが分かります。

set.seed(1013)
n <- 5000
M1 <- matrix(runif(n^2, min=1, max=10), n, n)
M2 <- matrix(runif(n^2, min=1, max=10), n, n)
system.time({ M3 <- M1 %*% M2 })

*1: 本当は古いバージョンのInno Setupが残っていたので、それで済ましています。

*2: インストール先をカスタマイズした場合は、Mkrules.localに、ISDIR = /path/to/Inno を付け加える必要があります。

*3: https://cran.r-project.org/src/base/R-latest.tar.gzを落としてくれば、リリース最新版がコンパイルできます。

*4: 無視すると「make[2]: *** 'stamp-recommended' に必要なターゲット 'MASS.ts' を make するルールがありません. 中止.」のようなことを言われて止まる.

*5: リリース版だと-rc無しで-4.2.1のようなバージョン番号がつきます。4.2リリース後は、R-patchedに展開されるようになりました。

*6:4.2リリース後はR-4.2.1patched-win.exe

反復測定分散分析を主観ベイズ推定に置き換えよう

TJO氏にふられた質問で、反復測定分散分析(repeated measures anova)をベイズ推定に置き換えられないかと言う話があったので、質問文から想定される状況を整理し、欠損値補定をFully Bayesian Imputation Method *1で、モデル比較をBayes Factorで行なうコードを描いて試してみました。

結論はRのlme4パッケージあたりで時点ダミー×介入群ダミーを入れた線形混合モデルを推定した方がコード量も計算時間も1/100で済むといった感じなのですが *2、各時点各群の平均と分散をそれぞれ推定できること、欠損値の補定方法のパラメーターも同時推定できる利点もあるので、ファンがうぃーんと言うのが好きな人には良いかもです。うぃーん。

1. 分析する状況の想定

反復測定分散分析が使われているところなので、医療系の実験で、薬剤を与える介入群Tと、偽薬を与える操作群Cの2群を、第1時点、第2時点、第3時点それぞれ調べ、介入群と操作群の違いが出る時点を探すことを考えます。
数式で書くとこんなモデルになります。

\mbox{score}_{gti} = \mu_{gt} + \beta x_i + \epsilon_{gt}

 \mbox{score}は検査指標、 \muは群の各時点の状態、 xは個体ごとの欠損値もある変数、 gは介入群と操作群を現し、 tは時点、 iは個体、 \epsilonは誤差項を示します。
こう書くと簡素なモデルが1つしかないと思うかも知れませんが*3 \mu_{gt} \mbox{VAR}(\epsilon_{gt}) で評価して、

  •  \mbox{T}_{t=1}=\mbox{C}_{t=1}, \mbox{T}_{t=2}=\mbox{C}_{t=2}, \mbox{T}_{t=3}=\mbox{C}_{t=3}(T群とG群にずっと差異なし)
  •  \mbox{T}_{t=1}=\mbox{C}_{t=1}, \mbox{T}_{t=2}=\mbox{C}_{t=2}, \mbox{T}_{t=3}\ne\mbox{C}_{t=3}(T群とG群は第3時点で差異が出た)
  •  \mbox{T}_{t=1}=\mbox{C}_{t=1}, \mbox{T}_{t=2}\ne\mbox{C}_{t=2}, \mbox{T}_{t=3}\ne\mbox{C}_{t=3}(T群とG群は第2時点から差異)
  •  \mbox{T}_{t=1}\ne\mbox{C}_{t=1}, \mbox{T}_{t=2}\ne\mbox{C}_{t=2}, \mbox{T}_{t=3}\ne\mbox{C}_{t=3}(T群とG群は第1時点から差異)

の4パターンのモデルを評価することができます。平均だけではなく分散も異なること、平均が同じでも分散が異なれば \mbox{T}_{t}\ne\mbox{C}_{t}となるとしましょう。また、どちらの群でも個体の状態が一定でないこと、 T_{t=1}\ne T_{t=2}\ne T_{t=3},C_{t=1}\ne C_{t=2}\ne C_{t=3}を仮定します。
欠損値の入り方は想像がつかなかったのですが、MCARとしておきます。原理的にはMARでも大丈夫ですし、欠損値が入るパターンが分かっていれば、推定モデルを工夫すればMNARにも対応できます。

2. データセット

真のモデルは以下の表の通りです。

パラメーター t=1 t=2 t=3
C  \mu_{Ct} 0.0 0.0 0.0
C  \mbox{VAR}(\epsilon_{Ct}) 1.0 1.0 1.0
T  \mu_{Tt} 0.0 1.0 2.0
T  \mbox{VAR}(\epsilon_{Tt}) 1.0 1.1 1.2

t=2時点からT群とC群に差があるとしています。
この設定を元に、架空のデータを生成します。練習感あふれてしまいますが、真のパラメーターが分かることから、推定結果が妥当なものか検討できるので、手法の確認には実データを使うより良いとも言えます。

set.seed(604)
rm(list=ls())

n <- 120 # 各時点のサンプルサイズ
t <- 1:3 # 時点
g <- c("C", "T") # 群

dimnames <- list(g, sprintf("t=%d", t))

# 以下のmu, sigma, betaは真の値になるので、本当は分からない
mu <- matrix(c(0+0*t, 0+1*(t-1)), ncol=length(t), byrow=TRUE, dimnames=dimnames)
sigma <- matrix(c(1+0*t, 1+0.1*(t-1)), ncol=length(t), byrow=TRUE, dimnames=dimnames)
beta <- 1

x <- rep(runif(n), length(t))

df01 <- data.frame(
    id=rep(1:n, length(t)),
    time=rep(t, each=n),
    group=rep(rep(g, each=n/length(g)), length(t)),
    x=x,
    z=x/2 + rnorm(n, sd=1),
    score=numeric(n*length(t))
)

rm(x)

for(j in t){
    for(i in 1:length(g)){
        s <- n*(j-1)+n/length(g)*(i-1) + 1
        e <- n*(j-1)+n/length(g)*i
        df01$score[s:e] <- rnorm(n/length(g), mean=mu[i, j], sd=sigma[i, j]) + beta*df01$x[s:e]
    }
}

# xに欠損値をつくる
df01$x[df01$id %in% sample(1:n, n/2)] <- NA

こんな感じのデータセットが出来上がります。

id time group x z score
1 1 C 0.03264839 -0.39328775 0.1545386
2 1 C 0.07093954 -0.07395381 -0.3282276
3 1 C 0.96846472 0.70290536 1.6737583
4 1 C NA -0.79364909 1.0485097
5 1 C 0.40517678 -0.50145817 0.7247011
6 1 C 0.09162524 -1.43459105 0.5523630

3. 推定コード

ベイズファクターの計算の都合上、推定はこれでもかと言うぐらいループが繰り返されます。

library(MCMCpack)

# ある点の確率密度を得る(Parzen Window)
get_density <- function(m, samp, n=200){
    max <- max(samp)
    min <- min(samp)
    range <- (max - min)/n
    sum(samp > (m - range) & samp < (m + range))/length(samp)/2/range
}

# Chib (1995)に習って、MCMCを繰り返して確率密度を計算する
get_posterior_density <- function(objf, init_theta, post.samp, estimated, log=TRUE){
    # 周辺尤度計算用の確率密度関数
    pdensity <- numeric(length(init_theta))
    len <- length(estimated)
    pdensity[len] <- get_density(estimated[len], post.samp[,len])
    for(i in 1:(len - 1)){
        l <- len - i
        post.samp.cnd <- MCMCmetrop1R(objf, theta.init=init_theta[1:l], force.samp=TRUE, estimated=estimated)
        pdensity[l] <- get_density(estimated[l], post.samp.cnd[,l])
    }
    ifelse(log, sum(log(pdensity)), prod(pdensity))
}

# xの係数と欠損値補定用の計4個のパラメーター
# t*g*2=12個のパラメーター
theta2param <- function(theta){
    theta <- theta[-(1:4)] # xと欠損値ダミーの係数のパラメーターは不要
  mlen <- length(t)*length(g)*2 # 制約なしの長さ
  tlen <- length(theta) # 引数の長さ
  slen <- (mlen - tlen)/2
  if(0 < slen){
    # 制約がついている
    # slen=kならば、t=k時点まではCとTは同じ
    len <- tlen/2
    m <- c(theta[1:length(t)], theta[1:slen])
    s <- c(theta[(1+len):(length(t)+len)], theta[(1+len):(slen+len)])
    if((length(t)+1)<=len){
      m <- c(m, theta[(length(t)+1):len])
      s <- c(s, theta[(length(t)+1+len):(2*len)])
    }
    theta <- c(m, s)
  }
  len <- length(theta)/2
  mu <- matrix(theta[1:len], ncol=length(t), byrow=TRUE, dimnames=dimnames)
  sigma <- matrix(theta[(len+1):(2*len)], ncol=length(t), byrow=TRUE, dimname=dimnames)
  list(mu=mu, sigma=sigma)
}
# theta2param(c(beta, 0, 0, 1, t(mu), t(sigma)))

# 尤度関数
llf.nd <- function(theta){
    p <- theta2param(theta)
    beta <- theta[1]
    alpha <- theta[2:4]
    theta <- theta[-(1:4)]
    len <- length(theta)/2
   # 補定処理の係数の尤度を計算
    xi <- df01$x
    s <- sum(dnorm(df01$x[!is.na(xi)] - alpha[1] - alpha[2]*df01$z[!is.na(xi)], sd=alpha[3], log=TRUE))
    # 線形の補定処理を行なう
    xi[is.na(xi)] <- alpha[1] + alpha[2]*df01$z[is.na(xi)]
    for(j in t){
        for(i in g){
            # i群j時点のscoreを説明するモデルの方の尤度を計算して加算
            s <- s + sum(dnorm(with(df01, { score[group==i & time==j] - beta*xi[group==i & time==j] })
                , mean = p$mu[i, j]
                , sd = p$sigma[i, j]
                , log = TRUE))
        }
    }
    s
}

# 事前分布(期待値には正規分布,偏差にはガンマ分布)
prior.nd <- function(theta){
    beta <- theta[1]
    alpha <- theta[2:4]
    theta <- theta[-(1:4)]
  len <- length(theta)/2
  mu <- 1:len
  sigma <- (1+len):(2*len)
  sum(dnorm(mu, mean=1, sd=1, log=TRUE)
    ,dgamma(sigma, shape=1.5, scale=2, log=TRUE)
        ,dnorm(beta, mean=1, sd=1, log=TRUE)
        ,dnorm(alpha[1], mean=0, sd=1, log=TRUE)
        ,dnorm(alpha[2], mean=0, sd=1, log=TRUE)
    ,dgamma(alpha[3], shape=1.5, scale=2, log=TRUE))
}

# 正規分布モデルの目的関数
objf.nd <- function(theta, estimated=theta){
    # thetaの長さが短ければ、周辺尤度の計算過程なので、estimatedの値で補完する
    len <- length(estimated)
    d <- len - length(theta)
    if(0 < d){
        p <- (len-d+1):len
        theta[p] <- estimated[p]
    }
    llf.nd(theta) + prior.nd(theta)
}

calcPosterior <- function(l=16){
    MCMCmetrop1R(objf.nd, theta.init=rep(1, l))
}

getMerginalLikelihood <- function(post.nd, l){
    mean.nd <- summary(post.nd)$statistics[, "Mean"]
    pdensity.nd <- get_posterior_density(objf.nd, rep(1, l), post.nd, mean.nd)
    llf.nd(mean.nd) - pdensity.nd
}

ptrn <- seq(2*max(t), 0, -2) + 10 # パラメーターの数は16(t=1に変化), 14(t=2に変化), 12(t=3に変化), 10(変化なし)でモデルの違いも意味する

# 並行処理の準備
library(parallel)
library(doParallel)
library(foreach)
cl <- makeCluster(max(1, detectCores() - 1))
registerDoParallel(cl)

# 対数周辺尤度を作る
# モデルごとに(パラメーター数, 事後分布, 周辺尤度)のリストを返す
r <- foreach(i=ptrn, .packages = c("MCMCpack")) %dopar% {
    post.name <- "posterior"
    mll.name <- "marginal likelihood"
    r <- list(i, calcPosterior(i))
    names(r) <- c("nop", post.name)
    r[mll.name] <- getMerginalLikelihood(r[[post.name]], i)
    return(r)
}

# 並行処理の終了
stopCluster(cl)

GetBayesFactor <- function(mll.a, mll.b){
    exp(mll.a - mll.b)
}

printMeans <- function(r){
    cat(sprintf("パラメーターの数:%d\n", r[[1]]))
    cat(sprintf("パラメーターの点推定値(mean):\n"))
    print(theta2param(summary(r[[2]])$statistics[,"Mean"]))
}

len <- length(ptrn)
bf <- matrix(NA, len, len, dimnames=list(c(ptrn), c(ptrn)))
for(i in 1:len){
    rownames(bf)[i] <- sprintf("numerator:t=%d", (max(ptrn) + 2 - r[[i]][[1]])/2)
    for(j in 1:len){
        colnames(bf)[j] <- sprintf("denominator:t=%d", (max(ptrn) + 2 - r[[j]][[1]])/2)
        bf[i, j] <- GetBayesFactor(r[[i]][[3]], r[[j]][[3]])
    }
}
# 最後の行と列は介入群と対照群の変化なし
rownames(bf)[len] <- colnames(bf)[len] <- "no change"
print(bf)

# Bayes Factorsからモデル選択
choiced_model_no <- which(bf[,1]==max(bf[,1]))
choiced_model <- r[[choiced_model_no]]
# 選択されたモデルの点推定値を表示
printMeans(choiced_model)
# 選択されたモデルの要約
summary(choiced_model[[2]])
# 選択されたモデルのパラメーターのプロット
plot(choiced_model[[2]])

長い? — 4つあるモデルを1つの処理で片付けようとしたために、ほんの少し複雑になっていますが、大したことありません。時点ダミー×介入群ダミーを入れた線形混合モデルだと10行以内に終わるとしてでもです。なお、ベイズファクターを用いているので、事前分布は絶対に省略しないでください。ところで記事を書いてから気付いたのですが、推定式に zを含めるのを忘れていますが、まぁ、含めても結果は変わらないと思うので気にしないでください。

4. 推定結果

各時点120とそれなりのサンプルサイズなのですが、ベイズファクターで選択されたモデルの推定結果はt=3で変化となりま・・・と最初書いていたのですが、コードにミスがあったので直したらt=1で変化となりました。真のモデルはt=2で変化のはずなので妙な気も済ますが、t=1とt=2のモデルのベイズ・ファクターを比較すると1.03とほとんど1であり、両者に差が無い事が分かります。t=2のモデルの方が簡素なので、第2時点で変化があったと解釈して良いでしょう。

ベイズファクターで選択されたモデル以外、今回の場合はt=2もチェックした方が良さそうです。今回のコードではprintMeans(r[[2]])で、もしくはMCMCPackを使っているので、summary(r[[2]][[2]])で数字を見ることができます。

なお、set.seed(611)にしておくと、t=2で変化と言う結果になります。例に使うシード値を間違えた感。

*1:"Bayesian methods for missing data Bayes Pharma"と言う文書を参考にしました。

*2:欠損値処理をしなければ、library(lme4); mixed.lmer <- lmer(score ~ x + factor(time)*factor(group) + (1|group), data = df01); print(summary(mixed.lmer))で終わります。多重代入法による欠損値処理も容易です( Rで一般化混合線形モデル×多重代入法 - 餡子付゛録゛ )。

*3:ごく普通に期待値と分散を線形モデルで説明する事も出来ます。

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

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

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

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

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

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

完全情報最尤推定法(FIML)で操作変数のある連立方程式モデルを推定

計量経済学では同時方程式を考えることが多く、そのために一般最小二乗法(OLS)の発展である操作変数法(IVM)とそのバリエーションである二段階最小二乗法(2SLS)*1や三段階最小二乗法(3SLS)、モーメント法の発展である一般化モーメント法(GMM)がよく使われているのですが、これらは尤度函数を用いないのでベイズ推定にできません。頻度主義者がベイジアンになるのは難しいです。
しかし、完全情報最尤推定法(FIML)と言うのがあり、実際のところは簡単に最尤法で同時方程式を推定できます。FIMLはRでは欠損値処理アルゴリズムのように扱われている気がするのですが、教科書的な手法なので確認しておきましょう*2。最尤法だけに不均一分散や欠損値処理の対処も工夫できますし、ベイズ推定しない場合でも、便利に使える局面は多そうです。

理論モデルとOLS推定量に入る同時性バイアス

推定の例をつくるわけですが、需要曲線と供給曲線を考えましょう*3


需要関数: d = \alpha_{10} + \alpha_{11} p + \xi_1
供給関数: s = \alpha_{20} + \alpha_{21} p  + \alpha_{22} z + \xi_2
均衡条件: d=s
 dは需要、 sは供給、 pは価格、 zは天候などの操作変数、 \alphaは係数、 \xiは誤差項です。ここで均衡価格を計算すると、

均衡価格: p^*=\frac{\alpha_{20} - \alpha_{10}}{\alpha_{11} - \alpha_{21}} + \frac{ \alpha_{22}}{\alpha_{11} - \alpha_{21}} z  + \frac{ \xi_2 - \xi_1 }{\alpha_{11} - \alpha_{21}}
となります。観測データは均衡価格 p^*になるわけですが、例えば需要関数の推定で p^*を説明変数に用いると、 p^*の第3項に \xi_1があるため、 COV(p^*, \xi_1) \neq 0となりOLSの前提を満たしません。需要関数の価格の係数の推定量 \hat{\alpha_{01}}_{OLS}=\alpha_{01}+COV(p^*, \xi_1) とバイアスが入ります。

完全情報最尤推定

需要関数と均衡価格式を同時推定することを考えましょう。推定する係数を \beta、誤差項を \muと置きなおし、方程式一本ごとに誤差項があると考え、二本の方程式を行列にまとめて書くと、次のようになります。 nは観測数です。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix}  = \begin{pmatrix}
1 & z_1 & p_1\\
1 & z_2 & p_2\\
\vdots & \vdots & \vdots \\
1 & z_n & p_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \\
\beta_{31}  & 0 \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

前節の係数とは、 \beta_{11}=\alpha_{10} \beta_{12}=(\alpha_{20} - \alpha_{10})/(\alpha_{11} - \alpha_{21})と言った対応関係があります。
右辺の内生変数と外生変数を分離します。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix}  = \begin{pmatrix}
1 & z_1 \\
1 & z_2 \\
\vdots & \vdots \\
1 & z_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \end{pmatrix} + \begin{pmatrix}
p_1\\
p_2\\
\vdots \\
p_n
\end{pmatrix}  \begin{pmatrix}
\beta_{31}  & 0 \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

内生変数を左辺にまとめます。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix} - \begin{pmatrix}
p_1\\
p_2\\
\vdots \\
p_n
\end{pmatrix}  \begin{pmatrix}
\beta_{31}  & 0 \end{pmatrix} = \begin{pmatrix}
1 & z_1 \\
1 & z_2 \\
\vdots & \vdots \\
1 & z_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

左辺第1項と第2項をまとめましょう。

 \begin{pmatrix}
d_1 & p_1 \\
d_2 & p_2 \\
\vdots & \vdots \\
d_n & p_n \end{pmatrix} \begin{pmatrix}
1 & 0 \\
 - \beta_{31}   & 1 \end{pmatrix} = \begin{pmatrix}
1 & z_1 \\
1 & z_2 \\
\vdots & \vdots \\
1 & z_n
\end{pmatrix}  \begin{pmatrix}
\beta_{11}  & \beta_{12}  \\
0  & \beta_{22} \end{pmatrix} +  \begin{pmatrix}
\mu_{11}  & \mu_{12}\\
\mu_{21}  & \mu_{22}\\
\vdots & \vdots \\
\mu_{n1} & \mu_{n2} \end{pmatrix}

この行列をいちいち書いていられないので、書き直します。

 Y\Gamma = W B + U

FIMLでは、この構造方程式を推定します。こう書くとOLSで入るバイアスが入らない理由が分かりづらいですが、この構造方程式の尤度函数は説明変数に内生変数が入らない誘導形 Y = W B \Gamma^{-1} + U \Gamma^{-1} = W\Pi + Vの尤度函数から導かれます*4。また、数理的な特性を使って整理すると以下のようなconcentrated loglikelihood function*5にできることが知られています*6


 \Sigma = \frac{1}{n} (Y\Gamma - WB)^{t} (Y\Gamma - WB)
 \log L(\beta) = \frac{-gn}{2}(\log 2\pi + 1)  + n \log | \Gamma | - \frac{n}{2} | \Sigma |
 gは内生変数(もしくは左辺にまとめられた被説明変)の数で、ここでは2です。 \beta = ( \beta_{11},  \beta_{12},  \beta_{22}, -\beta_{31} ) の4つのパラメーターを上の対数尤度函数で推定すればよいです。

データ

理論モデルに沿った式で、乱数を使ってデータを生成します。

set.seed(1231)
n <- 300
a10 <- 100
a11 <- -1
a20 <- 2
a21 <- 3
a22 <- 4
z <- runif(n, min=0, max=3)
mu <- rnorm(n, mean=0, sd=2)
nu <- rnorm(n, mean=0, sd=1)
p <- (a20 - a10 + a22*z + nu - mu)/(a11 - a21)
d <- a10 + a11*p + mu
s <- a20 + a21*p + a22*z + nu

推定

concentrated loglikelihood functionをoptim関数で推定します。

Y <- cbind(d, p)
W <- cbind(rep(1, n), z)
g <- 2

llf <- function(p){
    B <- matrix(c(p[1], 0, p[2], p[3]), 2, 2)
    G <- matrix(c(1, p[4], 0, 1), 2, 2) # Γ
    E <- Y%*%G - W%*%B
    S <- t(E)%*%E/n # ∑
    -g*n/2*(log(2*pi) + 1) + n*log(det(G)) - n/2*log(det(S))
}

r_optim <- optim(c(1, 0, 0, 0), function(p){ -llf(p) }, method="BFGS")

r_optim$parが推定された \hat{\beta}になりますが、内生変数の係数となる4番目のパラメーターの符号が逆になることに注意してください。なお、optim関数の引数にhessian=TRUEをつけるとヘッセ行列がとれるので、標準誤差を計算することができます*7

OLSとIVとの比較

推定結果の比較をしておきましょう。FIMLはOLSよりはかなりマシ、IVよりは僅かに劣るといった感じでしょうか。

係数 真の値 OLS IV FIML
 a_{10} 100 75.95501 100.11595 100.1228135
 a_{11} -1 0.03949 -1.01101 -1.0113085

なお、最尤法は漸近的に有効だけに、サンプルサイズが小さいとIVとの差が広がる傾向があります。IV推定量は以下のように計算できるので、nを変えて試してみてください。

zm <- matrix(c(rep.int(1, n), c(z)), n, 2)
xm <- matrix(c(rep.int(1, n), c(p)), n, 2)
solve(t(zm) %*% xm) %*% (t(zm) %*% d)

*1:操作変数が1つの場合はIVMと同値なのですが、一応、別なので。

*2:昔読んだGreene (2003)やWooldridge (2002)では最尤法で連立方程式モデルの推定の仕方は書いていなかったのですが、Davidson and MacKinnon (2003)では第12章で詳しく説明されています。

*3:Hayashi (2000)にある例を踏襲しています。

*4:Davidson and MacKinnon (2003)もしくはオンラインで配布している2021年版のpp.503–505とpp.544–545を参照。誘導形はSURと同じ多変量正規分布の尤度函数になるので、 \Gammaによる変形の影響を整理すれば構造形の尤度函数が導出できます。

*5:定訳不明

*6:一階条件の節の一階条件に直接関係のないExercise 12.24を \Sigma \Sigma^{-1}=Iに注意して \Sigma^{-1}の成分を考えて解くことで得ることができる・・・ハズ。

*7:Rと手作業で覚える最尤法 - 餡子付゛録゛

MCMCpackのMCMCmetrop1RでBayes Factorを雑に計算

マルコフ連鎖モンテカルロ法(MCMC)を使いたいベイジアンの統計ユーザーはStanを使っている人が多く、R使いもフロントエンドRStanを通して利用していることがほとんどだと思うのですが、MCMCpackと言うのもあります。このMCMCpackでベイズ因子(Bayes Factor)を使ってモデル選択をしてみましょう。

MCMCpackにはよく使われるらしい特定の共役事前分布*1を置ける線形回帰やポアソン回帰などの関数が用意されており、それらはベイズ因子を簡単に出力できるように設計されています。しかし、今回は正規分布と対数正規分布のどちらが当てはまりがよいかモデル選択をするので使えません。そして、汎用的なベイズ因子の計算関数は無いようです。

詰まった感があるのですが、推定モデルごとに周辺尤度を計算できればベイズ因子は簡単に計算できます。そして、事後分布の推定が終わっていれば、ベイズの定理から周辺尤度を簡単に計算できます。精度を考えなければ。
ベイズの定理の事後確率 P(\hat{\theta}|X) と周辺尤度 \int P(X|\theta) P(\theta) d\thetaを入れ替えてみましょう。

 \int P(X|\theta) P(\theta) d\theta = \frac{ P(X|\hat{\theta})P(\hat{\theta})}{P(\hat{\theta}|X)}

 Xはデータと言うか観測値のベクトル、\thetaはパラメーターのベクトル、P(\cdots)は確率、P(\cdots|\cdots)は条件付き確率です。 \hat{\theta}は任意の \thetaですが、精度の都合でMAP定量など、なるべく大きな値を選びます。 P(X|\hat{\theta}) \hat{\theta}で計算した尤度、 P(\hat{\theta}) \hat{\theta}の事前確率です。

1. データセット

今回は以下の対数正規分布からの乱数を使います。

set.seed(20211204)
n <- 15
y <- rlnorm(n, meanlog=1.2, sdlog=0.8)

真のパラメーターが分かっている方が、推定がうまくいっているか分かりやすいので。

ヒストグラムを描くと、こんな感じです。

2. 対数尤度関数

正規分布と対数正規分布の対数尤度関数は以下になります。

# 正規分布の対数尤度関数
llf_nd <- function(theta){
  if(any(theta < 0.0)){
    return(-Inf)
  }
  sum(dnorm(y, mean=theta[1], sd=theta[2], log=TRUE))
}

# 対数正規分布の対数尤度関数
llf_lnd <- function(theta){
  if(any(theta < 0.0)){
    return(-Inf)
  }
  sum(dlnorm(y, meanlog=theta[1], sdlog=theta[2], log=TRUE))
}

3. MCMCで推定

MCMCpackをインストールしてあったら、簡単にできます。なお、MCMCmetrop1R関数は初期値をoptimでMAP推定してからMCMCを開始するため、最適化アルゴリズムがMAP推定に失敗するとエラーを出します。そしてデフォルト指定のBFGSがエラーを出してきたので、L-BFGS-Bに代えました。

init_theta <- c(mean(y), sd(y))
library(MCMCpack)

# 正規分布のモデルの事前分布は(0 7)と(0 5)の一様分布
prior.nd <- function(theta){
    dunif(theta[1], min=0, max=7, log=TRUE) + dunif(theta[2], min=0, max=5, log=TRUE)
}
post.samp.nd <- MCMCmetrop1R(function(theta){
  llf_nd(theta) + prior.nd(theta)
}, theta.init=log(init_theta), optim.method="L-BFGS-B")

# 対数正規分布のモデルの事前分布は(0 log(7))と(0 log(5))の一様分布
prior.lnd <- function(theta){
    dunif(theta[1], min=0, max=log(7), log=TRUE) + dunif(theta[2], min=0, max=log(5), log=TRUE)
}
post.samp.lnd <- MCMCmetrop1R(function(theta){
  llf_lnd(theta) + prior.lnd(theta)
}, theta.init=log(init_theta))

4. 中央値から周辺尤度を計算

今回は手軽に中央値を使います。ノンパラメトリック確率密度関数の推定法…とかっこよく言いたいところですが、MCMCを繰り返して精度をあげるChib (1995)のような事はしていないので雑です。

# 中央値を得る
medians.nd <- summary(post.samp.nd)$quantiles[, "50%"]
medians.lnd <- summary(post.samp.lnd)$quantiles[, "50%"]

# パラメーターyに対応する事後確率を求める
get_posterior_probability <- function(post.samp, y, nos=100){
    p <- length(y) # 次元
    nr <- length(post.samp[, 1])
    count <- function(a){
        v <- 1
        flag <- rep(TRUE, nr)
        for(column in 1:p){
            s <- post.samp[, column]
            range <- max(s) - min(s)
            # y±range*aに含まれる点を選ぶ
            flag <- flag & s <= y[column] + range*a & s >= y[column] - range*a
            v <- v * (2*range*a)
        }
        c("nflag"=sum(flag), "volume"=v)
    }
    # 区間に入る点がnos個になるように、区間を調整する
    l <- 0
    h <- 1
    m <- 0.5
    nz <- h
    for(i in 1:64){
        m <- (l + h)/2
        n <- count(m)[["nflag"]] - nos
        if(n == 0){
            break;
        } else if(n < 0){
            l <- m
        } else {
            h <- m
            nz <- m
        }
    }
    r <- count(ifelse(0<=n, m, nz))
    # r["nflag"]/length(post.samp)は確率(累積分布)になるので、区間の大きさr["volume"]で割って(平均)確率密度にする
    r[["nflag"]]/length(post.samp)/r[["volume"]]
}

# 中央値の事後確率
posterior_p.nd <- get_posterior_probability(post.samp.nd, medians.nd)
posterior_p.lnd <- get_posterior_probability(post.samp.lnd, medians.lnd)

# 対数周辺尤度を作る
mll.nd <- llf_nd(medians.nd)-log(posterior_p.nd) + prior.nd(medians.nd)
mll.lnd <- llf_lnd(medians.lnd)-log(posterior_p.lnd) + prior.lnd(medians.lnd)

5. ベイズ因子を計算

割り算と言うか、引き算して指数をとって終わりです。

sprintf("Bayes Factor: %.3f", exp(mll.lnd-mll.nd))


8.967でSubstantialと言えるような値になりますが、ヒストグラムの見た目にあった結果ではないでしょうか。なお、シャピロ・ウィルク検定をかけると10%有意で正規分布が棄却されて、対数正規分布は棄却されません。

shapiro.test(x=y)$p.value
shapiro.test(x=log(y))$p.value

6. 注意事項

ほぼ思いつきで書いた胡散臭い方法なので、実用の参考にはしないでください*2

*1: MCMCregress関数だと、b0, B0, c0, d0, sigma.mu, sigma.varで事前分布になる多変量正規分布、逆ガンマ分布のパラメーターを指定できます。

*2: Chib (1995)の方法で計算しなおしたら7.154になりました。8.967は目安として使えなくもないですが、論文に書くと怒られると思います。

Rの主成分回帰で多重共線性を制御

一般線形回帰(OLS)でコスト分析をしたら係数がマイナスの項目が出てきた、工数をかけるほど費用が減っていく・・・と言う嘆きを見かけたのですが、説明変数の相関が高いと推定される係数の標準誤差が大きくなって推定量がおかしくなる、多重共線性による問題だと考えられます*1

こういう問題を解決する素朴な方法として、主成分回帰(Principal Component Regression)で多重共線性を制御する技があるので紹介します。試してみると望んだ結果が得られるかも知れないので、学部生で卒論に使う回帰分析の係数が気に入らない人にp-hacking的にお勧めです。なお、べたべたとコードを書いていきますが、plsパッケージで手軽に実行でき、さらに同様の目的に使えるPLS回帰も試せるので、基本的にはplsパッケージを使う方が賢いです。

さて、まず解決すべき多重共線性をつくってみましょう。最後に比較するときに、真のモデルの切片項と係数はそれぞれ1であることに注意してください。

set.seed(20201115)
n <- 150 # 十分なサンプルサイズ(e.g. 1500)にすると多重共線性が問題なくなる
x1 <- runif(n)
x2 <- x1 + rnorm(n, sd=0.1)
y <- 1 + x1 + x2 + rnorm(n)

VIFは8.02675です。
計算結果を格納する行列をつくります。

beta <- matrix(0, 3, 2)
se <- matrix(0, 3, 2)
colnames(beta) <- colnames(se) <- c("OLS", "PCR")
rownames(beta) <- rownames(se) <- c("Const.", "x1", "x2")

比較のためにOLSで回帰します。

lm_r <- lm(y ~ x1 + x2)
beta[, "OLS"] <- coef(summary(lm_r))[, 1]
se[, "OLS"] <- coef(summary(lm_r))[, 2]

主成分分析を行います。scale=TRUEを忘れないようにしましょう。

r_prcomp <- prcomp(~ x1 + x2, scale=TRUE)

寄与率を低い主成分を省略し、回帰を行います。今回は、summary(r_prcomp)してCumulative Proportionが0.9を超えているので、第1主成分だけで行います*2。回帰に用いる主成分の数と説明変数の数が一致すると、OLSと同じ予測力になります。

npc <- 1 # 回帰に用いる主成分の数
lm_prc_r <- lm(y ~ predict(r_prcomp)[,1:npc])

第1主成分を、例えば平均全体工数のような抽象的な概念に対応するとして分析を進めてもよいのですが、そのような概念をあまり説明したくない場合は、主成分の係数を、元の説明変数の係数に分解しましょう。第1主成分の回帰係数と第1主成分の固有ベクトル(説明変数を標準化した値の係数)から、元の説明変数の回帰係数を求めます*3

beta["x1", "PCR"] <- sum(r_prcomp$rotation[1,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x1))
beta["x2", "PCR"] <- sum(r_prcomp$rotation[2,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x2))
beta["Const.", "PCR"] <- coef(lm_prc_r)[1] - beta["x1", "PCR"]*mean(x1) - beta["x2", "PCR"]*mean(x2)

続けて分散の加法公式を使って、上の計算手順をもとに、元の説明変数の擬似的な標準誤差を求めます。

se["x1", "PCR"] <- sqrt(sum((r_prcomp$rotation[1, ][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x1) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2))
se["x2", "PCR"] <- sqrt(sum((r_prcomp$rotation[2, ][1:npc]*coef(lm_prc_r)[2:(npc + 1)]/sd(x2) * coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2))
se["Const.", "PCR"] <- sqrt(
    coef(summary(lm_prc_r))[, 2][1]^2
    + sum((r_prcomp$rotation[1,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]*mean(x1)/sd(x1)*coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2)
    + sum((r_prcomp$rotation[2,][1:npc]*coef(lm_prc_r)[2:(npc + 1)]*mean(x2)/sd(x2)*coef(summary(lm_prc_r))[, 2][2:(npc + 1)])^2)
)

OLSとPCRの推定量と残差平方和を比較しておきましょう。

# 推定量を比較する
cat("\nestimated coeffcients:\n")
print(beta)

cat("\nstandard errors:\n")
print(se)

# 残差平方和を比較する
cat(
    sprintf("\nOLSの残差平方和: %.3f\n主成分回帰の残差平方和: %.3f\n", 
        sum((y - predict(lm_r))^2), 
        sum((y - model.matrix(~ x1 + x2) %*% beta[, 2])^2)
    )
)

係数 真のモデル OLS PCR
Const. 1 1.10562969 1.0806319
x1 1 -0.01407916 0.9303908
x2 1 1.80575174 0.8960958
標準誤差 OLS PCR
Const. 0.1775728 0.09022349
x1 0.8366178 0.05424751
x2 0.8057794 0.05224790

OLSの残差平方和は143.823、PCRの残差平方和は145.11と予測力は僅かに悪化しますが、係数を見るとPCRの推定量の方が真のモデルに近い事がわかります。

*1:さらにomitted variable biasで、見えない成分を代理してしまっている可能性もあります。例えば、常連客との取引では低コストになる仕組みがあり、常連客によく提供するサービス項目の工数であったりすると、常連客のコスト低減効果を代理することもありえます。

*2:省略する成分は、偶発的に生じた見かけ上の、もしくは想定外だが従属変数への影響が大きくない構造(e.g. 上述の常連客のコスト低減効果)になります。

*3: C_{pcr}を主成分回帰の切片工、 \beta_{pcr}を主成分回帰の回帰式の係数、 r_iを主成分の説明変数 x_iの係数、主成分回帰の誤差項を \etaとすると、主成分回帰は  y = C_{pcr} + \beta_{pcr} \bigg( \sum^2_{i=1} r_i \frac{x_i - \bar{x_i}}{sd(x_i)} \bigg) + \eta と書けるので、式を x_iの項で分けて整理すると、OLSと同様の説明変数とその係数の式に直せます。