統計解析R演習10
演習の概要
- マルチレベル・モデル
注意事項
- 慶應義塾大学SFCで開講している「統計解析」の授業履修者向けの演習用ページです。
- 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.6.0を使っています。
- 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を説明変数とする線形回帰モデルをベースとし、マルチレベル・モデルと比較する。- マルチレベル・モデルでは、都県によるランダム効果を推定する。
lph <- read.csv("lph.csv", header=TRUE, sep=",")
head(lph)線形回帰モデル（重回帰モデル）の最小二乗推定
- まずは線形回帰モデルを
lm()を用いて推定する
lph.lm1 <- lm(LPH~POPD+EMP3D,data=lph)
summary(lph.lm1)- 単回帰モデルと比較する
lph.lm2 <- lm(LPH~POPD,data=lph)
summary(lph.lm2)anova()の適用
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)- BICの比較
BIC(lph.lm1, lph.lm2)- ベイズファクターの計算
BF_BIC = exp((BIC(lph.lm1) - BIC(lph.lm2))/2)
BF_BIC[1] 5.44654e-93マルチレベル・モデル
マルチレベル・モデルの最尤推定
- マルチレベル・モデルは、最尤推定法を用いて推定する。パッケージ
lme4のlmer()を用いれば、lm()と同様の手順でモデル推定できる。 - ここでは、都県毎にランダムなモデルを推定する。
- 傾きが固定で切片がランダムなモデル
lph.lme1 <- lmer(LPH~-1+POPD+EMP3D+(1|AREA),
data=lph)
summary(lph.lme1)
ranef(lph.lme1)- 傾きがランダムで切片が固定されたモデル
lph.lme2 <- lmer(LPH~1+(-1+POPD+EMP3D|AREA),
data=lph)
summary(lph.lme2)
ranef(lph.lme2)- 傾きと切片がランダムなモデル
lph.lme3 <- lmer(LPH~1+POPD+EMP3D+(1+POPD+EMP3D|AREA),
data=lph)
summary(lph.lme3)
ranef(lph.lme3)- モデルの適合度を
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 + POPD + EMP3D + (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.00 1 1
lph.lme3 10 2527.5 2566.4 -1253.8 2507.5 62.12 5 4.429e-12 ***
---
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)There were 111 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.85 4.01 9.22 24.84 817 1.00
sd(POPD) 1.06 0.35 0.53 1.86 610 1.00
sd(EMP3D) 1.17 0.68 0.25 2.83 664 1.01
cor(Intercept,POPD) -0.81 0.16 -0.98 -0.37 1520 1.00
cor(Intercept,EMP3D) 0.03 0.37 -0.61 0.77 1400 1.00
cor(POPD,EMP3D) 0.01 0.40 -0.75 0.75 1501 1.00
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept 8.37 4.78 -1.09 17.97 441 1.00
POPD 1.22 0.43 0.40 2.12 618 1.00
EMP3D 1.53 0.85 -0.00 3.42 1094 1.00
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma 7.19 0.28 6.69 7.79 2989 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)
LS0tCnRpdGxlOiAi57Wx6KiI6Kej5p6QUua8lOe/kjEwIgojYXV0aG9yOiAiVG9tb3l1a2kgRnVydXRhbmkiCiNkYXRlOiAiMjAxOS84LzE0IgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdAogIGh0bWxfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQprbml0cjo6b3B0c19jaHVuayRzZXQoY2FjaGUgPSBUUlVFKQprbml0cjo6b3B0c19jaHVuayRzZXQoZXJyb3IgPSBUUlVFKQpgYGAKCiMjIOa8lOe/kuOBruamguimgQotIOODnuODq+ODgeODrOODmeODq+ODu+ODouODh+ODqwoKIyMjIOazqOaEj+S6i+mghQotIOaFtuaHiee+qeWhvuWkp+WtplNGQ+OBp+mWi+ism+OBl+OBpuOBhOOCi+OAjOe1seioiOino+aekOOAjeOBruaOiOalreWxpeS/ruiAheWQkeOBkeOBrua8lOe/kueUqOODmuODvOOCuOOBp+OBmeOAggotIOW/heOBmuOBl+OCguWFqOOBpuOBruODkOODvOOCuOODp+ODs+OBrlLjgoRPU+OBp+WLleS9nOeiuuiqjeOCkuihjOOBo+OBpuOBhOOBvuOBm+OCk+OAguOBk+OBrua8lOe/kueUqOODmuODvOOCuOOCkuS9nOaIkOOBl+OBpuOBhOOCi+autemajuOBp+OBr+OAgVIzLjYuMOOCkuS9v+OBo+OBpuOBhOOBvuOBmeOAggotIFLjga7mm7TmlrDjgarjganjgavjgojjgorjgIFS44Kz44O844OJ44Gq44Gp44KS5LqI5ZGK54Sh44GX44Gr5aSJ5pu044GZ44KL5aC05ZCI44GM44GC44KK44G+44GZ44CCIAoKIyMg5YiG5p6Q44Gu5rqW5YKZCiMjIyDjg5Hjg4PjgrHjg7zjgrjjga7jgqTjg7Pjgrnjg4jjg7zjg6sKYGBge3IgaW5zdGFsbCBwYWNrYWdlcywgZXZhbD1GQUxTRX0KaW5zdGFsbC5wYWNrYWdlcygibG1lNCIpCmluc3RhbGwucGFja2FnZXMoImxtZXJUZXN0IikKaW5zdGFsbC5wYWNrYWdlcygiYnJtcyIpCmluc3RhbGwucGFja2FnZXMoInNwIikKaW5zdGFsbC5wYWNrYWdlcygidGlkeXZlcnNlIikKaW5zdGFsbC5wYWNrYWdlcygic2pQbG90IikKaW5zdGFsbC5wYWNrYWdlcygic2ptaXNjIikKaW5zdGFsbC5wYWNrYWdlcygic2psYWJlbGxlZCIpCmBgYAoKYGBge3IgbGlicmFyeSwgZXZhbD1GQUxTRX0KbGlicmFyeShsbWU0KQpsaWJyYXJ5KGxtZXJUZXN0KQpsaWJyYXJ5KGJybXMpCmxpYnJhcnkoc3ApCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHNqUGxvdCkKbGlicmFyeShzam1pc2MpCmxpYnJhcnkoc2psYWJlbGxlZCkKYGBgCgojIyMg44OH44O844K/44Gu5rqW5YKZCi0g6aaW6YO95ZyP44Gu5L2P5a6F5Zyw5Zyw5L6h44Gu44OH44O844K/44KS44KC44Gh44GE44G+44GZ44CCW+OBk+OBk10oaHR0cDovL3dlYi5zZmMua2Vpby5hYy5qcC9+bWF1bnovRFNCMTkvRFNCZGF0YS56aXAp44GL44KJ44OA44Km44Oz44Ot44O844OJ44Gn44GN44G+44GZ44CCCi0g44CM5L2P5a6F5Zyw5Zyw5L6h77yI5LiH5YaGL+OOoe+8ieOAjeOCkmBMUEhg44CB44CM5bi45L2P5Lq65Y+j5a+G5bqm77yI5Y2D5Lq6L2ttXjIp44CN44KSYFBPUERg44CB44CM56ys5LiJ5qyh55Sj5qWt5b6T5qWt5Lq65Y+j77yI5Y2D5Lq6L2ttXjIp44CN44KSYEVNUDNEYOOBqOihqOOBmeOAguOBk+OBruS7luOBq+OAgemDveecjOWQjUFSRUFg44CB6Z2i56mNYFNg44CB6YO955yM44GU44Go44Go44Gu44K744Kw44Oh44Oz44OIYHNlZ2DjgYzjgYLjgovjgIIKLSBgTFBIYOOCkuiiq+iqrOaYjuWkieaVsOOAgWBQT1BEYOOBqGBFTVAzRGDjgpLoqqzmmI7lpInmlbDjgajjgZnjgovnt5rlvaLlm57luLDjg6Ljg4fjg6vjgpLjg5njg7zjgrnjgajjgZfjgIHjg57jg6vjg4Hjg6zjg5njg6vjg7vjg6Ljg4fjg6vjgajmr5TovIPjgZnjgovjgIIKLSDjg57jg6vjg4Hjg6zjg5njg6vjg7vjg6Ljg4fjg6vjgafjga/jgIHpg73nnIzjgavjgojjgovjg6njg7Pjg4Djg6DlirnmnpzjgpLmjqjlrprjgZnjgovjgIIKCmBgYHtyIHJlYWRscGgxLCBldmFsPUZBTFNFfQpscGggPC0gcmVhZC5jc3YoImxwaC5jc3YiLCBoZWFkZXI9VFJVRSwgc2VwPSIsIikKaGVhZChscGgpCmBgYAoKIyMjIOe3muW9ouWbnuW4sOODouODh+ODq++8iOmHjeWbnuW4sOODouODh+ODq++8ieOBruacgOWwj+S6jOS5l+aOqOWumgotIOOBvuOBmuOBr+e3muW9ouWbnuW4sOODouODh+ODq+OCkmBsbSgpYOOCkueUqOOBhOOBpuaOqOWumuOBmeOCiwpgYGB7ciBscGgubG0xLTEsIGV2YWw9RkFMU0V9CmxwaC5sbTEgPC0gbG0oTFBIflBPUEQrRU1QM0QsZGF0YT1scGgpCnN1bW1hcnkobHBoLmxtMSkKYGBgCi0g5Y2Y5Zue5biw44Oi44OH44Or44Go5q+U6LyD44GZ44KLCmBgYHtyIGxwaC5sbTItMSwgZXZhbD1GQUxTRX0KbHBoLmxtMiA8LSBsbShMUEh+UE9QRCxkYXRhPWxwaCkKc3VtbWFyeShscGgubG0yKQpgYGAKLSBgYW5vdmEoKWDjga7pgannlKgKYGBge3IgbHBoLmxtLmFub3ZhfQphbm92YShscGgubG0xLCBscGgubG0yKQpgYGAKLSBBSUPjga7mr5TovIMKYGBge3IgbHBoLmxtLkFJQ30KQUlDKGxwaC5sbTEsIGxwaC5sbTIpCmBgYAotIEJJQ+OBruavlOi8gwpgYGB7ciBscGgubG0uQklDfQpCSUMobHBoLmxtMSwgbHBoLmxtMikKYGBgCi0g44OZ44Kk44K644OV44Kh44Kv44K/44O844Gu6KiI566XCmBgYHtyIGxwaC5sbS5CRn0KQkZfQklDID0gZXhwKChCSUMobHBoLmxtMSkgLSBCSUMobHBoLmxtMikpLzIpCkJGX0JJQwpgYGAKCiMjIOODnuODq+ODgeODrOODmeODq+ODu+ODouODh+ODqwojIyMg44Oe44Or44OB44Os44OZ44Or44O744Oi44OH44Or44Gu5pyA5bCk5o6o5a6aCi0g44Oe44Or44OB44Os44OZ44Or44O744Oi44OH44Or44Gv44CB5pyA5bCk5o6o5a6a5rOV44KS55So44GE44Gm5o6o5a6a44GZ44KL44CC44OR44OD44Kx44O844K4YGxtZTRg44GuYGxtZXIoKWDjgpLnlKjjgYTjgozjgbDjgIFgbG0oKWDjgajlkIzmp5jjga7miYvpoIbjgafjg6Ljg4fjg6vmjqjlrprjgafjgY3jgovjgIIKLSDjgZPjgZPjgafjga/jgIHpg73nnIzmr47jgavjg6njg7Pjg4Djg6Djgarjg6Ljg4fjg6vjgpLmjqjlrprjgZnjgovjgIIKCjEuIOWCvuOBjeOBjOWbuuWumuOBp+WIh+eJh+OBjOODqeODs+ODgOODoOOBquODouODh+ODqwpgYGB7ciBscGgubG1lMS0xLCBldmFsPUZBTFNFfQpscGgubG1lMSA8LSBsbWVyKExQSH4tMStQT1BEK0VNUDNEKygxfEFSRUEpLAogICAgICAgICAgICAgICAgIGRhdGE9bHBoKQpzdW1tYXJ5KGxwaC5sbWUxKQpyYW5lZihscGgubG1lMSkKYGBgCgoyLiDlgr7jgY3jgYzjg6njg7Pjg4Djg6DjgafliIfniYfjgYzlm7rlrprjgZXjgozjgZ/jg6Ljg4fjg6sKYGBge3IgbHBoLmxtZTItMSwgZXZhbD1GQUxTRX0KbHBoLmxtZTIgPC0gbG1lcihMUEh+MSsoLTErUE9QRCtFTVAzRHxBUkVBKSwKICAgICAgICAgICAgICAgICBkYXRhPWxwaCkKc3VtbWFyeShscGgubG1lMikKcmFuZWYobHBoLmxtZTIpCmBgYAoKMy4g5YK+44GN44Go5YiH54mH44GM44Op44Oz44OA44Og44Gq44Oi44OH44OrCmBgYHtyIGxwaC5sbWUzLTEsIGV2YWw9RkFMU0V9CmxwaC5sbWUzIDwtIGxtZXIoTFBIfjErUE9QRCtFTVAzRCsoMStQT1BEK0VNUDNEfEFSRUEpLAogICAgICAgICAgICAgICAgIGRhdGE9bHBoKQpzdW1tYXJ5KGxwaC5sbWUzKQpyYW5lZihscGgubG1lMykKYGBgCgotIOODouODh+ODq+OBrumBqeWQiOW6puOCkmBhbm92YSgpYOOBp+avlOi8g+OBmeOCiwpgYGB7ciBscG4uYW5vdmF9CmFub3ZhKGxwaC5sbWUxLCBscGgubG1lMiwgbHBoLmxtZTMpCmBgYAoKIyMjIOODnuODq+ODgeODrOODmeODq+ODouODh+ODq+OBruODmeOCpOOCuuaOqOWumgotIOODnuODq+ODgeODrOODmeODq+ODouODh+ODq+OBruODmeOCpOOCuuaOqOWumuOBr+OAgWBicm1zYOODkeODg+OCseODvOOCuOOBrmBicm0oKWDjgpLnlKjjgYTjgabmjqjlrprjgafjgY3jgovjgIIKLSDjgZPjgZPjgafjga/jgIHkuIrov7Djga5scGgubG1lM++8iOWCvuOBjeOBqOWIh+eJh+OBjOOBqOOCguOBq+ODqeODs+ODgOODoOOBquODouODh+ODq++8ieOCkuODmeOCpOOCuuaOqOWumuOBl+OCiOOBhuOAggotIOOBvuOBmuOAgeS6i+WJjeWIhuW4g+OCkuaMh+WumuOBm+OBmuOAgU1DTUPmnJ/plpPjgpIyMDAw5Zue44CBYnVybmlu77yI44Km44Kp44O844Og44Ki44OD44OX77yJ5pyf6ZaT44KSMTAwMOWbnuOAgeODgeOCp+ODvOODs+aVsO+8k+OBqOOBl+OBpuioreWumuOBmeOCi+OBqOOAgeS7peS4i+OBruOCiOOBhuOBq+OBquOCi+OAggoKYGBge3IgbHBoLmJybTMtMSwgZXZhbD1GQUxTRX0KbHBoLmJybTMgPC0gYnJtKExQSH4xK1BPUEQrRU1QM0QrKDErUE9QRCtFTVAzRHxBUkVBKSwKICAgICAgICAgICAgICAgICBkYXRhPWxwaCwKICAgICAgICAgICAgICAgIHByaW9yPU5VTEwsCiAgICAgICAgICAgICAgICBpdGVyPTIwMDAsCiAgICAgICAgICAgICAgICB3YXJtdXA9MTAwMCwKICAgICAgICAgICAgICAgIGNoYWlucz0zCiAgICAgICAgICAgICAgICApCmBgYAoKLSDmjqjlrprntZDmnpzjga/ku6XkuIvjga7jgojjgYbjgavjgarjgovjgIIKYGBge3IgbHBoLmJybTMuc3VtbWFyeX0Kc3VtbWFyeShscGguYnJtMykKYGBgCgotIOS6i+W+jOWIhuW4g+OCkuWbs+ekuuOBmeOCi+OAggpgYGB7ciBscGguYnJtMy5wbG90fQpwbG90KGxwaC5icm0zKQpgYGAKCi0g44OR44Op44Oh44O844K/44Gu5LqL5b6M5bmz5Z2H44Go5LqL5b6M5qiZ5rqW5YGP5beu44KS5Zuz56S644GZ44KL44CCCmBgYHtyIGxwaC5icm0zLnN0YW5wbG90fQpzdGFucGxvdChscGguYnJtMykKYGBgCgo=