演習の概要

  • 線形回帰モデルのベイズ推定を行う
  • パッケージMCMCpackとbrmsを用いる

注意事項

  • 慶應義塾大学SFCで開講している「ベイズ統計」の授業履修者向けの演習用ページです。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
  • Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

線形回帰モデルのベイズ推定

パッケージのインストール

  • install.packages()でパッケージをインストールし、library()でパッケージを呼び出す
  • MCMCpackはMCMCで様々な回帰モデルをRのみで実行するのに汎用的に用いられるパッケージ。https://cran.r-project.org/web/packages/MCMCpack/MCMCpack.pdf
  • brmsはStanを使って様々な回帰モデルをベイズ推定するのに用いられる。自身でStanコードを書く必要がなく、モデルを指定すると自動的にStanコードが作成され、サンプリングが生成される。https://cran.r-project.org/web/packages/brms/brms.pdf
install.packages("MCMCpack")
install.packages("brms")
install.packages("sp")
install.packages("coda")
library(MCMCpack)
library(brms)
library(sp)
library(coda)

データの読み込み

  • ここでは、『家計調査』2017年9月、第２−６表「年間収入階級別１世帯当たり１か月間の収入と支出」「二人以上の世帯」から「年間年収階級」の18階級について、「消費支出(千円)」を被説明変数、「年間収入(10万円)」と「有業人員(人)」を説明変数とする線形回帰モデル（重回帰モデル）を推定します
  • データは政府統計ポータルサイトe-statから入手可能です
  • 「消費支出(千円)」をCONS1、「年間収入(10万円)」をINC1、「有業人員(人)」をWORKという変数で表す
  • 分析用のデータはここから入手可能です。データを各自の作業用フォルダにダウンロードしてください。
  • RまたはRStudio上で、PC上の作業用フォルダをRのワーキングディレクトリとして設定してください。setwd("作業フォルダのパス")で作業フォルダのパスを指定して設定するか、RStudioの「Session」→「Set Working Directory」→「Choose Directory」で設定できます。
kakei <- read.csv("kakei201709_l18.csv", header=TRUE, sep=",")
  • 最小二乗法を用いて線形回帰モデルを推定してみましょう。lm()を用いて以下のように推定します。
kakei.lm1 <- lm(CONS1~INC1+WORK, data=kakei)
summary(kakei.lm1)
## 
## Call:
## lm(formula = CONS1 ~ INC1 + WORK, data = kakei)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8671 -0.9996 -0.1032  1.3300  1.9492 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 12.03768    1.11306  10.815 1.77e-08 ***
## INC1         0.18704    0.01419  13.181 1.19e-09 ***
## WORK         2.98212    1.19404   2.498   0.0246 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.539 on 15 degrees of freedom
## Multiple R-squared:  0.9793, Adjusted R-squared:  0.9765 
## F-statistic: 354.9 on 2 and 15 DF,  p-value: 2.337e-13

MCMCpackを用いる方法

  • MCMCpackで線形回帰モデルのベイズ推定を行うには、MCMCregress()を用いて回帰式、バーンイン期間(burnin)、事後分布のMCMC期間(mcmc)、チェーン数(chain)などを指定する
kakei.mcmc1 <- MCMCregress(
  CONS1~INC1+WORK, 
  data=kakei, 
  burnin = 1000, 
  mcmc = 1000, 
  thin = 1,
  chain = 3)
summary(kakei.mcmc1)
## 
## Iterations = 1001:2000
## Thinning interval = 1 
## Number of chains = 1 
## Sample size per chain = 1000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean      SD Naive SE Time-series SE
## (Intercept) 12.0744 1.18082 0.037341        0.03734
## INC1         0.1874 0.01499 0.000474        0.00052
## WORK         2.9331 1.24651 0.039418        0.03942
## sigma2       2.7050 1.29458 0.040938        0.04834
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%    50%     75%   97.5%
## (Intercept) 9.7259 11.3071 12.097 12.8168 14.4324
## INC1        0.1571  0.1782  0.188  0.1971  0.2162
## WORK        0.4141  2.1234  2.897  3.6946  5.4244
## sigma2      1.2803  1.8778  2.412  3.1751  5.8529
  • summary()を用いて、事後平均、事後標準偏差、95%革新区間などの要約を示すことができる
  • plot()を用いれば、事後分布とそのサンプリング結果を可視化できる
plot(kakei.mcmc1)

  • codaパッケージを用いて、ベイズ推定結果を診断することができる
  • geweke.diag()はGewekeの診断結果を示す
# Gewekeの診断方法
geweke.diag(kakei.mcmc1)
## 
## Fraction in 1st window = 0.1
## Fraction in 2nd window = 0.5 
## 
## (Intercept)        INC1        WORK      sigma2 
##     -0.5775     -1.0475      1.3016     -1.9118
  • raftery.diag()はRaftery-Lewisの診断結果を示す
# Raftery-Lewisの診断方法
raftery.diag(kakei.mcmc1)
## 
## Quantile (q) = 0.025
## Accuracy (r) = +/- 0.005
## Probability (s) = 0.95 
## 
## You need a sample size of at least 3746 with these values of q, r and s

brmsを用いる方法

  • brmsで線形回帰モデルのベイズ推定を行うには、brm()を用いて回帰式、バーンイン期間（warm-up）、事後分布のMCMC期間(iter)、チェーン数(chain)などを指定する
kakei.brm1 <- brm(CONS1~INC1+WORK, 
                  data=kakei,
                  family=gaussian,
                  iter=2000,
                  warmup=1000,
                  seed=1234,
                  chain=3
                  )
summary(kakei.brm1)

  • MCMCの計算途中経過は以下のように出力される brm1-0.png

  • モデル推定結果は以下のように出力される brm1-1.png

  • plot()を用いて事後分布とMCMCの挙動を示すことができる

plot(kakei.brm1)