演習の概要

注意事項

  • 慶應義塾大学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を被説明変数、POPDEMP3Dを説明変数とする線形回帰モデルをベースとし、マルチレベル・モデルと比較する。
  • マルチレベル・モデルでは、都県によるランダム効果を推定する。
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

マルチレベル・モデル

マルチレベル・モデルの最尤推定

  • マルチレベル・モデルは、最尤推定法を用いて推定する。パッケージlme4lmer()を用いれば、lm()と同様の手順でモデル推定できる。
  • ここでは、都県毎にランダムなモデルを推定する。
  1. 傾きが固定で切片がランダムなモデル
lph.lme1 <- lmer(LPH~-1+POPD+EMP3D+(1|AREA),
                 data=lph)
summary(lph.lme1)
ranef(lph.lme1)
  1. 傾きがランダムで切片が固定されたモデル
lph.lme2 <- lmer(LPH~1+(-1+POPD+EMP3D|AREA),
                 data=lph)
summary(lph.lme2)
ranef(lph.lme2)
  1. 傾きと切片がランダムなモデル
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回、チェーン数3として設定すると、以下のようになる。
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)