ベイズ統計R演習6
演習の概要
- マルチレベル・モデルの最尤推定とベイズ推定を行う
- パッケージ
lme4とbrmsを用いる
注意事項
- 慶應義塾大学SFCで開講している「ベイズ統計」の授業履修者向けの演習用ページです。
- 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
- Rの更新などにより、Rコードなどを予告無しに変更する場合があります。
分析の準備
パッケージのインストール
install.packages("lme4")
install.packages("lmerTest")
install.packages("brms")
install.packages("sp")
install.packages("tidyverse")
install.packages("sjPlot")
install.packages("sjmisc")
install.packages("sjlabelled")library(lme4)
library(lmerTest)
library(brms)
library(sp)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(sjlabelled)データの準備
- 首都圏の住宅地地価のデータをもちいます。ここからダウンロードできます。
- 「住宅地地価（万円/㎡）」を
LPH、「常住人口密度（千人/km2)」をPOPD、「第三次産業従業人口（千人/km2)」をEMP3Dと表す。この他に、都県名AREA、面積S、都県ごととのセグメントseg`がある。 LPHを被説明変数、POPDとEMP3Dを説明変数とする線形回帰モデルをベースとし、マルチレベル・モデルと比較する。- マルチレベル・モデルでは、都県によるランダム効果を推定する。
lpj <- read.csv("lph.csv", header=TRUE, sep=",")
head(lph)## ID JCODE Easting Northing LPH POPD EMPD EMP3D AREA S seg
## 1 1 8201 54029.151 41292.37 6.96 1.425 0.840 0.731 Ibaraki 184 1
## 2 2 8202 70243.203 69620.67 5.23 2.060 1.009 0.621 Ibaraki 97 1
## 3 3 8203 32224.406 11260.67 4.02 1.427 0.779 0.595 Ibaraki 101 1
## 4 4 8204 -5202.486 19960.62 4.75 1.278 0.540 0.333 Ibaraki 114 1
## 5 5 8205 34447.897 25470.16 2.80 0.604 0.221 0.149 Ibaraki 136 1
## 6 6 8207 2762.423 28982.98 3.87 0.832 0.343 0.205 Ibaraki 63 1線形回帰モデル（重回帰モデル）の最小二乗推定
- まずは線形回帰モデルを
lm()を用いて推定する
lph.lm1 <- lm(LPH~POPD+EMP3D,data=lph)
summary(lph.lm1)##
## Call:
## lm(formula = LPH ~ POPD + EMP3D, data = lph)
##
## Residuals:
## Min 1Q Median 3Q Max
## -85.110 -1.873 -0.581 1.230 67.060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.55112 0.56647 4.504 9.05e-06 ***
## POPD 1.68162 0.10611 15.848 < 2e-16 ***
## EMP3D 2.24666 0.07841 28.654 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.188 on 359 degrees of freedom
## Multiple R-squared: 0.8289, Adjusted R-squared: 0.8279
## F-statistic: 869.6 on 2 and 359 DF, p-value: < 2.2e-16- 単回帰モデルと比較する
lph.lm2 <- lm(LPH~POPD,data=lph)
summary(lph.lm2)##
## Call:
## lm(formula = LPH ~ POPD, data = lph)
##
## Residuals:
## Min 1Q Median 3Q Max
## -28.990 -2.180 -1.087 0.300 215.036
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.9532 1.0249 1.906 0.0575 .
## POPD 2.9311 0.1751 16.736 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.82 on 360 degrees of freedom
## Multiple R-squared: 0.4376, Adjusted R-squared: 0.436
## F-statistic: 280.1 on 1 and 360 DF, p-value: < 2.2e-16anova()の適用
anova(lph.lm1, lph.lm2)## Analysis of Variance Table
##
## Model 1: LPH ~ POPD + EMP3D
## Model 2: LPH ~ POPD
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 359 24068
## 2 360 79116 -1 -55047 821.07 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- AICの比較
AIC(lph.lm1, lph.lm2)## df AIC
## lph.lm1 4 2554.631
## lph.lm2 3 2983.413- BICの比較
BIC(lph.lm1, lph.lm2)## df BIC
## lph.lm1 4 2570.197
## lph.lm2 3 2995.088- ベイズファクターの計算
BF_BIC = exp((BIC(lph.lm1) - BIC(lph.lm2))/2)
BF_BIC## [1] 5.44654e-93線形回帰モデルのベイズ推定
brmsパッケージのbrm()を用いて線形回帰モデルをベイズ推定する
lph.brm1 <- brm(LPH~POPD+EMP3D,
data=lph,
prior=NULL,
iter=2000,
warmup=1000,
chains=3,
save_all_pars = TRUE, #ベイズファクターを計算する際に必要
seed=1234
)- 推定結果は以下のようになる。
summary(lph.brm1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: LPH ~ POPD + EMP3D
## Data: lph (Number of observations: 362)
## Samples: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 3000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 2.54 0.57 1.42 3.67 3395 1.00
## POPD 1.68 0.11 1.47 1.89 2694 1.00
## EMP3D 2.25 0.08 2.09 2.39 2423 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 8.20 0.30 7.63 8.85 2878 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).- 事後分布を図示する。
plot(lph.brm1)
- パラメータの事後平均と事後標準偏差を図示する。
stanplot(lph.brm1)
- 次に、単回帰モデルをベイズ推定する
lph.brm2 <- brm(LPH~POPD,
data=lph,
prior=NULL,
iter=2000,
warmup=1000,
chains=3,
save_all_pars = TRUE, #ベイズファクターを計算する際に必要
seed=1234
)- 推定結果は以下のようになる。
summary(lph.brm2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: LPH ~ POPD
## Data: lph (Number of observations: 362)
## Samples: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 3000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 1.92 1.05 -0.17 3.93 3183 1.00
## POPD 2.93 0.18 2.58 3.29 3479 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 14.83 0.56 13.82 15.99 2814 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).- 事後分布を図示する。
plot(lph.brm2)
- パラメータの事後平均と事後標準偏差を図示する。
stanplot(lph.brm2)
lph.brm2とlph.brm2のベイズファクターを計算する
brms::bayes_factor(lph.brm1,lph.brm2)## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 1
## Iteration: 2
## Iteration: 3
## Iteration: 4
## Iteration: 5## Estimated Bayes factor in favor of bridge1 over bridge2: 230532886857468028871920528800930260131844698434920290150177761044211256165880601476870438912.00000マルチレベル・モデル
マルチレベル・モデルの最尤推定
- マルチレベル・モデルは、最尤推定法を用いて推定する。パッケージ
lme4のlmer()を用いれば、lm()と同様の手順でモデル推定できる。 - ここでは、都県毎にランダムなモデルを推定する。
- 傾きが固定で切片がランダムなモデル
lph.lme1 <- lmer(LPH~-1+POPD+EMP3D+(1|AREA),
data=lph)
summary(lph.lme1)
ranef(lph.lme1)## Linear mixed model fit by REML ['lmerMod']
## Formula: LPH ~ -1 + POPD + EMP3D + (1 | AREA)
## Data: lph
##
## REML criterion at convergence: 2558.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -10.4909 -0.1517 -0.0033 0.1748 8.5265
##
## Random effects:
## Groups Name Variance Std.Dev.
## AREA (Intercept) 24.05 4.905
## Residual 63.09 7.943
## Number of obs: 362, groups: AREA, 12
##
## Fixed effects:
## Estimate Std. Error t value
## POPD 1.43112 0.14163 10.10
## EMP3D 2.11984 0.08504 24.93
##
## Correlation of Fixed Effects:
## POPD
## EMP3D -0.063## $AREA
## (Intercept)
## Chiba 1.855681
## Chiba-shi 1.468511
## Gunma 1.258787
## Ibaraki 1.695953
## Kanagawa 5.318527
## Kawasaki-shi 3.127584
## Saitama 3.601225
## Saitama-shi 2.378541
## Tochigi 2.176931
## Tokyo-ku 11.355791
## Tokyo-tama 6.438640
## Yokohama-shi 1.563193
##
## with conditional variances for "AREA"- 傾きがランダムで切片が固定されたモデル
lph.lme2 <- lmer(LPH~1+(-1+POPD+EMP3D|AREA),
data=lph)
summary(lph.lme2)
ranef(lph.lme2)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.00719226
## (tol = 0.002, component 1)## Linear mixed model fit by REML ['lmerMod']
## Formula: LPH ~ 1 + (-1 + POPD + EMP3D | AREA)
## Data: lph
##
## REML criterion at convergence: 2568.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.0021 -0.2694 -0.0666 0.1464 8.2451
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## AREA POPD 2.237 1.496
## EMP3D 2.672 1.635 0.90
## Residual 62.898 7.931
## Number of obs: 362, groups: AREA, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.6045 0.6368 5.66
## convergence code: 0
## Model failed to converge with max|grad| = 0.00719226 (tol = 0.002, component 1)## $AREA
## POPD EMP3D
## Chiba 1.50570585 1.69189687
## Chiba-shi 1.24226936 1.30625115
## Gunma 0.02764093 0.03557282
## Ibaraki 0.30016081 0.32412166
## Kanagawa 1.78984019 1.85846410
## Kawasaki-shi 1.49290423 1.64264607
## Saitama 1.57472333 1.63091140
## Saitama-shi 1.39486463 1.35770973
## Tochigi 0.26021382 0.27302613
## Tokyo-ku 1.76208097 2.25151556
## Tokyo-tama 1.95883877 2.07184964
## Yokohama-shi 1.52280394 0.78354923
##
## with conditional variances for "AREA"- 傾きと切片がランダムなモデル
lph.lme3 <- lmer(LPH~1+POPD+EMP3D+(1+POPD+EMP3D|AREA),
data=lph)
summary(lph.lme3)
ranef(lph.lme3)## boundary (singular) fit: see ?isSingular## Linear mixed model fit by REML ['lmerMod']
## Formula: LPH ~ 1 + POPD + EMP3D + (1 + POPD + EMP3D | AREA)
## Data: lph
##
## REML criterion at convergence: 2505.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -11.4296 -0.1739 -0.0318 0.1387 7.3704
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## AREA (Intercept) 294.023 17.1471
## POPD 1.053 1.0262 -0.97
## EMP3D 0.495 0.7035 -0.08 0.30
## Residual 50.651 7.1170
## Number of obs: 362, groups: AREA, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.3246 5.1105 2.020
## POPD 1.1855 0.3293 3.600
## EMP3D 1.4387 0.3258 4.416
##
## Correlation of Fixed Effects:
## (Intr) POPD
## POPD -0.932
## EMP3D -0.111 0.257
## convergence code: 0
## boundary (singular) fit: see ?isSingular## $AREA
## (Intercept) POPD EMP3D
## Chiba -9.119506 0.68795922 0.50360221
## Chiba-shi -7.303220 0.39373251 -0.07465558
## Gunma -8.997025 0.52075878 0.01659625
## Ibaraki -8.580667 0.51551193 0.07314475
## Kanagawa -3.506601 0.17611270 -0.07517377
## Kawasaki-shi 3.918439 -0.32942099 -0.31921078
## Saitama -6.695836 0.42201635 0.11709827
## Saitama-shi -2.163649 -0.01584167 -0.42491810
## Tochigi -7.952898 0.46892605 0.04082446
## Tokyo-ku 50.184929 -2.76264706 0.33950269
## Tokyo-tama -6.986157 0.68324930 0.86076195
## Yokohama-shi 7.202190 -0.76035713 -1.05757234
##
## with conditional variances for "AREA"- モデルの適合度を
anova()で比較する
anova(lph.lme1, lph.lme2, lph.lme3)## refitting model(s) with ML (instead of REML)## Data: lph
## Models:
## lph.lme1: LPH ~ -1 + POPD + EMP3D + (1 | AREA)
## lph.lme2: LPH ~ 1 + (-1 + POPD + EMP3D | AREA)
## lph.lme3: LPH ~ -1 + (1 + POPD + EMP3D | AREA)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## lph.lme1 4 2561.7 2577.2 -1276.8 2553.7
## lph.lme2 5 2579.6 2599.1 -1284.8 2569.6 0.000 1 1
## lph.lme3 7 2551.7 2578.9 -1268.8 2537.7 31.965 2 1.145e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1マルチレベルモデルのベイズ推定
- マルチレベルモデルのベイズ推定は、
brmsパッケージのbrm()を用いて推定できる。 - ここでは、上述のlph.lme3（傾きと切片がともにランダムなモデル）をベイズ推定しよう。
- まず、事前分布を指定せず、MCMC期間を2000回、burnin（ウォームアップ）期間を1000回、チェーン数３として設定すると、以下のようになる。
lph.brm3 <- brm(LPH~1+POPD+EMP3D+(1+POPD+EMP3D|AREA),
data=lph,
prior=NULL,
iter=2000,
warmup=1000,
chains=3
)- 推定結果は以下のようになる。
summary(lph.brm3)## Warning: There were 249 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
## See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: LPH ~ 1 + POPD + EMP3D + (1 + POPD + EMP3D | AREA)
## Data: lph (Number of observations: 362)
## Samples: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 3000
##
## Group-Level Effects:
## ~AREA (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 15.92 4.05 9.24 25.61 771 1.01
## sd(POPD) 1.08 0.35 0.54 1.86 791 1.00
## sd(EMP3D) 1.12 0.61 0.25 2.71 679 1.00
## cor(Intercept,POPD) -0.81 0.16 -0.98 -0.37 1061 1.00
## cor(Intercept,EMP3D) 0.03 0.37 -0.62 0.79 1180 1.00
## cor(POPD,EMP3D) 0.03 0.40 -0.73 0.73 1680 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 8.48 4.69 -0.59 17.60 333 1.01
## POPD 1.21 0.44 0.36 2.08 498 1.00
## EMP3D 1.53 0.82 0.05 3.39 920 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 7.18 0.28 6.66 7.77 1985 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).- 事後分布を図示する。
plot(lph.brm3)
- パラメータの事後平均と事後標準偏差を図示する。
stanplot(lph.brm3)
- 次に、事前分布を指定して、MCMC期間を2000回、burnin（ウォームアップ）期間を1000回、チェーン数３として設定すると、以下のようになる。
lph.brm4 <- brm(LPH~1+POPD+EMP3D+(1+POPD+EMP3D|AREA),
data=lph,
prior=c(set_prior("normal(0,10)", class = "b")),
iter=2000,
warmup=1000,
chains=3
)- 推定結果は以下のようになる。
summary(lph.brm4)## Warning: There were 118 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
## See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: LPH ~ 1 + POPD + EMP3D + (1 + POPD + EMP3D | AREA)
## Data: lph (Number of observations: 362)
## Samples: 3 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 3000
##
## Group-Level Effects:
## ~AREA (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 15.80 4.14 9.14 25.45 624 1.00
## sd(POPD) 1.07 0.31 0.56 1.78 763 1.01
## sd(EMP3D) 1.00 0.58 0.22 2.33 641 1.01
## cor(Intercept,POPD) -0.81 0.15 -0.98 -0.40 1406 1.00
## cor(Intercept,EMP3D) 0.10 0.35 -0.55 0.81 1398 1.00
## cor(POPD,EMP3D) -0.01 0.39 -0.76 0.69 1474 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 9.09 4.62 -0.14 18.62 480 1.00
## POPD 1.20 0.38 0.46 1.95 740 1.00
## EMP3D 1.32 0.57 0.31 2.56 1647 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 7.19 0.28 6.67 7.79 2350 1.00
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).- 事後分布を図示する。
plot(lph.brm4)
- パラメータの事後平均と事後標準偏差を図示する。
stanplot(lph.brm4)