事例：米国の160の学校に通う7185人の生徒の数学の成績に関するデータがある。 従属変数として、生徒の数学の成績、個人レベルの独立変数として（１）性別、（２）人種的マイノリティかどうか、（３）家庭の社会経済的地位、学校レベルの独立変数として（４）学校における差別の程度を設定し、各変数の効果を検証したい ## データの読み込みと準備
# データ取得のためパッケージを読み込む
library(nlme)
# nlmeパッケージの中に入っているサンプルデータを取り出す
dat.i <- data.frame(MathAchieve) # 個人レベルの属性
dat.s <- data.frame(MathAchSchool) # 学校レベルの属性
# データの中身を見てみる
head(dat.i) # 個人レベル（School変数で学校を識別できることが分かる）## School Minority Sex SES MathAch MEANSES
## 1 1224 No Female -1.528 5.876 -0.428
## 2 1224 No Female -0.588 19.708 -0.428
## 3 1224 No Male -0.528 20.349 -0.428
## 4 1224 No Male -0.668 8.781 -0.428
## 5 1224 No Male -0.158 17.898 -0.428
## 6 1224 No Male 0.022 4.583 -0.428head(dat.s) # 学校レベル## School Size Sector PRACAD DISCLIM HIMINTY MEANSES
## 1224 1224 842 Public 0.35 1.597 0 -0.428
## 1288 1288 1855 Public 0.27 0.174 0 0.128
## 1296 1296 1719 Public 0.32 -0.137 1 -0.420
## 1308 1308 716 Catholic 0.96 -0.622 0 0.534
## 1317 1317 455 Catholic 0.95 -1.694 1 0.351
## 1358 1358 1430 Public 0.25 1.535 0 -0.014# 両者のデータを、school変数をキーとして統合する
dat <- merge(x = dat.i,
y = dat.s,
by = "School")
head(dat) # 個人レベルのデータに学校情報の列が統合される## School Minority Sex SES MathAch MEANSES.x Size Sector PRACAD DISCLIM
## 1 1224 No Female -1.528 5.876 -0.428 842 Public 0.35 1.597
## 2 1224 No Female -0.588 19.708 -0.428 842 Public 0.35 1.597
## 3 1224 No Male -0.528 20.349 -0.428 842 Public 0.35 1.597
## 4 1224 No Male -0.668 8.781 -0.428 842 Public 0.35 1.597
## 5 1224 No Male -0.158 17.898 -0.428 842 Public 0.35 1.597
## 6 1224 No Male 0.022 4.583 -0.428 842 Public 0.35 1.597
## HIMINTY MEANSES.y
## 1 0 -0.428
## 2 0 -0.428
## 3 0 -0.428
## 4 0 -0.428
## 5 0 -0.428
## 6 0 -0.428教科書とは異なり、lme4パッケージを使用する
library(lme4) # マルチレベル分析用のパッケージ
library(texreg) # 結果を並べて見るのに便利なパッケージ学校のランダム切片のみのモデルを作る
out1.0 <- lmer(MathAch ~ 1|School, data = dat)
summary(out1.0) # 結果## Linear mixed model fit by REML ['lmerMod']
## Formula: MathAch ~ 1 | School
## Data: dat
##
## REML criterion at convergence: 47116.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0631 -0.7539 0.0267 0.7606 2.7426
##
## Random effects:
## Groups Name Variance Std.Dev.
## School (Intercept) 8.614 2.935
## Residual 39.148 6.257
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 12.6370 0.2444 51.71screenreg(out1.0) # texregパッケージのこちらの関数を用いると、より見やすい結果が見やすく表示される）##
## ======================================
## Model 1
## --------------------------------------
## (Intercept) 12.64 ***
## (0.24)
## --------------------------------------
## AIC 47122.79
## BIC 47143.43
## Log Likelihood -23558.40
## Num. obs. 7185
## Num. groups: School 160
## Var: School (Intercept) 8.61
## Var: Residual 39.15
## ======================================
## *** p < 0.001; ** p < 0.01; * p < 0.05covmat <- data.frame(summary(out1.0)$varcor) # データフレームに変換する
covmat # 中身を見てみる## grp var1 var2 vcov sdcor
## 1 School (Intercept) <NA> 8.614025 2.934966
## 2 Residual <NA> <NA> 39.148322 6.256862covmat$sdcor[1]^2/sum(covmat$sdcor^2) # Schoolの違いによってどの程度分散が説明されるか## [1] 0.1803518【参考】なお、級内相関は関数を使って直接計算することもできる
library(misty)
multilevel.icc(x = dat$MathAch, cluster = dat$School)out1.1 <- lmer(MathAch ~
Sex + SES + Minority + DISCLIM +
(1 + Minority|School),
data = dat)
summary(out1.1)## Linear mixed model fit by REML ['lmerMod']
## Formula: MathAch ~ Sex + SES + Minority + DISCLIM + (1 + Minority | School)
## Data: dat
##
## REML criterion at convergence: 46337.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2750 -0.7204 0.0366 0.7596 2.9994
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## School (Intercept) 2.413 1.553
## MinorityYes 1.624 1.274 -0.20
## Residual 35.724 5.977
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.1100 0.1758 80.256
## SexFemale -1.3164 0.1609 -8.182
## SES 2.0461 0.1054 19.415
## MinorityYes -3.0341 0.2333 -13.005
## DISCLIM -1.1496 0.1498 -7.673
##
## Correlation of Fixed Effects:
## (Intr) SexFml SES MnrtyY
## SexFemale -0.479
## SES -0.091 0.064
## MinorityYes -0.332 0.013 0.176
## DISCLIM 0.025 0.045 0.101 -0.014社会経済的地位の高いこどもが多い学校に行く効果を除くため、グループ平均をしたSES変数を投入して、効果を検証してみる。
# SESの学校ごとの平均を計算する
sesmean <- aggregate(x = dat.i$SES, # データベクトル
by = list(dat.i$School), # グループ
FUN = mean) # 関数
colnames(sesmean) <- c("School", "ses_m") # 変数名をつける
# データを統合する
dat2 <- merge(x = dat, y = sesmean,
by = "School")
# 個々のSESを所属学校の平均SESから引く（グループ平均で中心化）
dat2$ses_c <- dat2$SES - dat2$ses_m
# グループ平均で中心化した社会経済的地位を含めてマルチレベルを行う。
out1.2 <- lmer(MathAch ~
Sex + ses_c + DISCLIM + ses_m + Minority
+ (1 + Minority|School),
data = dat2)
# 結果
summary(out1.2)## Linear mixed model fit by REML ['lmerMod']
## Formula: MathAch ~ Sex + ses_c + DISCLIM + ses_m + Minority + (1 + Minority |
## School)
## Data: dat2
##
## REML criterion at convergence: 46302.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2703 -0.7199 0.0340 0.7593 3.0033
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## School (Intercept) 1.739 1.319
## MinorityYes 1.740 1.319 -0.14
## Residual 35.701 5.975
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.0505 0.1627 86.350
## SexFemale -1.2996 0.1596 -8.142
## ses_c 1.8947 0.1086 17.446
## DISCLIM -0.8151 0.1448 -5.630
## ses_m 4.1611 0.3511 11.850
## MinorityYes -2.8281 0.2360 -11.984
##
## Correlation of Fixed Effects:
## (Intr) SexFml ses_c DISCLI ses_m
## SexFemale -0.515
## ses_c -0.076 0.051
## DISCLIM -0.003 0.062 0.007
## ses_m -0.106 0.064 0.032 0.359
## MinorityYes -0.328 0.019 0.129 0.038 0.196マイノリティであることの効果が学校の差別の程度によって異なるかを調べるため、クロスレベルの交互作用を含めてみる。
out1.3 <- lmer(MathAch ~
Sex + ses_c + Minority + DISCLIM + ses_m + Minority*DISCLIM
+ (1 + Minority|School),
data = dat2)
summary(out1.3)## Linear mixed model fit by REML ['lmerMod']
## Formula: MathAch ~ Sex + ses_c + Minority + DISCLIM + ses_m + Minority *
## DISCLIM + (1 + Minority | School)
## Data: dat2
##
## REML criterion at convergence: 46282.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2297 -0.7212 0.0341 0.7580 3.1200
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## School (Intercept) 1.7773 1.333
## MinorityYes 0.8446 0.919 0.01
## Residual 35.6573 5.971
## Number of obs: 7185, groups: School, 160
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 14.0932 0.1631 86.421
## SexFemale -1.3003 0.1596 -8.148
## ses_c 1.8806 0.1085 17.340
## MinorityYes -2.8533 0.2201 -12.966
## DISCLIM -0.5514 0.1561 -3.533
## ses_m 4.1137 0.3532 11.647
## MinorityYes:DISCLIM -1.0325 0.2163 -4.773
##
## Correlation of Fixed Effects:
## (Intr) SexFml ses_c MnrtyY DISCLI ses_m
## SexFemale -0.515
## ses_c -0.078 0.052
## MinorityYes -0.310 0.020 0.142
## DISCLIM 0.016 0.053 -0.010 0.025
## ses_m -0.097 0.064 0.034 0.201 0.330
## MnY:DISCLIM -0.046 0.012 0.047 0.041 -0.350 0.017モデルを並べて結果を比較してみる。
screenreg(list(out1.0, out1.1, out1.2, out1.3))##
## ===============================================================================================
## Model 1 Model 2 Model 3 Model 4
## -----------------------------------------------------------------------------------------------
## (Intercept) 12.64 *** 14.11 *** 14.05 *** 14.09 ***
## (0.24) (0.18) (0.16) (0.16)
## SexFemale -1.32 *** -1.30 *** -1.30 ***
## (0.16) (0.16) (0.16)
## SES 2.05 ***
## (0.11)
## MinorityYes -3.03 *** -2.83 *** -2.85 ***
## (0.23) (0.24) (0.22)
## DISCLIM -1.15 *** -0.82 *** -0.55 ***
## (0.15) (0.14) (0.16)
## ses_c 1.89 *** 1.88 ***
## (0.11) (0.11)
## ses_m 4.16 *** 4.11 ***
## (0.35) (0.35)
## MinorityYes:DISCLIM -1.03 ***
## (0.22)
## -----------------------------------------------------------------------------------------------
## AIC 47122.79 46355.73 46322.59 46304.64
## BIC 47143.43 46417.65 46391.39 46380.31
## Log Likelihood -23558.40 -23168.86 -23151.30 -23141.32
## Num. obs. 7185 7185 7185 7185
## Num. groups: School 160 160 160 160
## Var: School (Intercept) 8.61 2.41 1.74 1.78
## Var: Residual 39.15 35.72 35.70 35.66
## Var: School MinorityYes 1.62 1.74 0.84
## Cov: School (Intercept) MinorityYes -0.40 -0.25 0.01
## ===============================================================================================
## *** p < 0.001; ** p < 0.01; * p < 0.05最尤法を用いたマルチレベル分析は、モデルが複雑になるほどエラーが起きやすい。モデルを単純化することは一つの方策ではるが、モデルはそのままにしてbrmsというパッケージを用いて、ベイズ推定を用いて計算をさせてみるということも考えられる。ただし、Rでベイズ統計を扱えるようにするためには、以下のようなやや複雑な準備が必要となる。
自分のOSごとにインストール方法が異なる。こちらのページ（https://learnb4ss.github.io/learnB4SS/articles/install-brms.html） を参照し、自身のOSに従ってインストールを行うこと。
以下は、OSがWindowsの場合 の説明であるが、Rtoolsの版が変わっている場合には機能しない場合がある。あくまで参考程度のものとすること。
[Windowsの場合のインストール方法]
手順1：Rtoolsをインストールする(Windowsの場合) このページ（https://github.com/stan-dev/rstan/wiki/Configuring-C---Toolchain-for-Windows）を参考に最新のRToolsをインストール（この文書の執筆時点では、Rtools43）。文書真ん中あたりにRtools installerというリンクがあるので、そちらをクリックする。
手順2: Rtoolsをインストールする
install.packages("rstan")以下の警告が出る場合があるが、無視してよい。
# Warning message:
# In system(paste(CXX, ARGS), ignore.stdout = TRUE, ignore.stderr = TRUE) :
# 'C:/rtools40/usr/mingw_/bin/g++' not foundinstall.packages("brms")library(rstan)## 要求されたパッケージ StanHeaders をロード中です##
## rstan version 2.32.3 (Stan version 2.26.1)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file##
## 次のパッケージを付け加えます: 'rstan'## 以下のオブジェクトは 'package:texreg' からマスクされています:
##
## extractlibrary(brms) # ベイジアン・マルチレベル分析用のパッケージ## 要求されたパッケージ Rcpp をロード中です## Loading 'brms' package (version 2.20.4). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## 次のパッケージを付け加えます: 'brms'## 以下のオブジェクトは 'package:rstan' からマスクされています:
##
## loo## 以下のオブジェクトは 'package:lme4' からマスクされています:
##
## ngrps## 以下のオブジェクトは 'package:stats' からマスクされています:
##
## arlibrary(parallel) # 並列計算用のパッケージはじめに自分のパソコンがいくつ分の並列計算が可能なのかを調べてみる。
library(parallel)
detectCores()## [1] 40ここで、複数の数字が返されれば、brmsを回す際にその分の並列計算が可能であり時間短縮になる。以下では、４個分のコアが見つかった場合を想定し、それをコア数に設定している。
注意：ベイズモデルなのですぐに反応しない。慌てず待つこと。
out1.3brm <- brm(MathAch ~
Sex + ses_c + Minority + DISCLIM + ses_m + Minority*DISCLIM
+ (1 + Minority|School),
data = dat2,
prior = NULL, # 事前確率。NULLとして一様分布にしている
chains = 4, # chainの回数
iter = 1000, # 反復の回数。試行用に少なめに設定、もっと少なくても良い。逆に本番は最低5000は必要
warmup = 500, # 本番は最低2000は必要。
cores = 4) # 並列計算用のコア数 ## Compiling Stan program...## Start sampling## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-esssummary(out1.3brm)## Warning: There were 1 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: MathAch ~ Sex + ses_c + Minority + DISCLIM + ses_m + Minority * DISCLIM + (1 + Minority | School)
## Data: dat2 (Number of observations: 7185)
## Draws: 4 chains, each with iter = 1000; warmup = 500; thin = 1;
## total post-warmup draws = 2000
##
## Group-Level Effects:
## ~School (Number of levels: 160)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 1.35 0.13 1.10 1.60 1.00 912
## sd(MinorityYes) 0.67 0.38 0.04 1.41 1.01 312
## cor(Intercept,MinorityYes) 0.12 0.38 -0.61 0.88 1.01 1027
## Tail_ESS
## sd(Intercept) 1282
## sd(MinorityYes) 844
## cor(Intercept,MinorityYes) 916
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 14.09 0.17 13.78 14.42 1.00 1306 1665
## SexFemale -1.30 0.16 -1.62 -0.99 1.00 3013 1754
## ses_c 1.88 0.11 1.68 2.09 1.00 3790 1336
## MinorityYes -2.85 0.22 -3.28 -2.45 1.00 2201 1090
## DISCLIM -0.55 0.16 -0.86 -0.24 1.00 1163 1591
## ses_m 4.10 0.36 3.42 4.83 1.00 1127 1128
## MinorityYes:DISCLIM -1.04 0.21 -1.47 -0.62 1.00 1749 1636
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 5.98 0.05 5.88 6.08 1.00 3191 1431
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# ベイズモデルなので、メモリを大量に消費する。なるべく作業ごとにメモリを解放したほうがよい。
gc()## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 4280163 228.6 8394674 448.4 5024134 268.4
## Vcells 11394939 87.0 19932257 152.1 19932236 152.1参考文献
永吉希久子，2016，『行動科学の統計学：社会調査のデータ分析』共立出版．