d01 <- read.csv("practice.csv")『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』
第13章 階層構造のあるデータの分析
0 全体の章構成
1 マルチレヴェル分析
変数の準備
d01$Site <- factor(floor(d01$id/100)) # グループ変数（調査地点）
d01$SEX <- factor(d01$sex, levels = 1:2, labels = c("男性", "女性"))1-1 級内相関係数(Intra-Class Correlation coefficient)
このページでは，記号表記は尾崎・川端・山田編『Rで学ぶマルチレベルモデル〔入門編〕』（朝倉書店）を参考としつつ多少変更する。 \[\begin{align*} y_{ij} =& \color{green}\beta_{0j}\color{black}+\color{red}r_{ij} \\ &\color{green}\beta_{0j}\color{black} = \color{blue}\gamma_{00}\color{black}+\color{red}u_{0j} \\ y_{ij} =& \color{blue}\gamma_{00}\color{black}+\color{red}u_{0j}\color{black}+\color{red}{r_{ij}} \\ &\color{red}u_{0j}\color{black} \sim\;N(0,\tau_{00})\\ &\color{red}r_{ij}\color{black} \sim\;N(0,\sigma^2)\\ ICC &=\frac{\tau_{00}}{\tau_{00}+\sigma^2}\\ \end{align*}\]
模擬データの作成
グループによって平均水準だけ異なり，グループの内部では同じだけランダムに個人差のある模擬データを作成する。
s01 <- data.frame(group = rep(c(1:10), 30))
s01$outcome <- rnorm(nrow(s01), s01$group, sd = 5) + 10with(s01, boxplot(outcome ~ group,
main = "模擬乱数データの概要"))
変量効果モデルでの分析とICC
summary(lme00 <- lme4::lmer(outcome ~ 1 + (1|group), s01))Linear mixed model fit by REML ['lmerMod']
Formula: outcome ~ 1 + (1 | group)
Data: s01
REML criterion at convergence: 1812.3
Scaled residuals:
Min 1Q Median 3Q Max
-2.79806 -0.73990 0.03094 0.68605 2.93883
Random effects:
Groups Name Variance Std.Dev.
group (Intercept) 8.001 2.829
Residual 22.895 4.785
Number of obs: 300, groups: group, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 15.4005 0.9362 16.45print(std.dev00 <- lme4::VarCorr(lme00), comp = c("Variance", "Std.Dev")) Groups Name Variance Std.Dev.
group (Intercept) 8.0011 2.8286
Residual 22.8953 4.7849 d_std.dev00 <- as.data.frame((std.dev00))
d_std.dev00$vcov[1]/sum(d_std.dev00$vcov)[1] 0.2589657固定効果モデル（分散分析，一般線形モデル）で分析するとこうなる。
anova00 <- anova(aov(outcome ~ group, s01))
anova00R.squared <- anova00$"Sum Sq"[1]/sum(anova00$"Sum Sq")
R.squared[1] 0.2382034ICCと分散説明率が似た様な値になっている。
1-2 マルチレヴェルモデル
ランダム切片モデルの模擬データを作成する。
模擬データセットの生成スクリプト
n_group <- 10
n_within <- 30
group <- c(1:n_group)
# constant <- rnorm(n_group, group, sd = 5) + 10
constant <- group*2 + 10
indep <- NULL
for (i in 1:n_within) {
indep <- c(indep, rep(i, n_group))
}
s02 <- data.frame(group = rep(group, n_within),
constant = rep(constant, n_within),
factor = indep)
s02$outcome <- s02$constant + s02$factor * 0.5 + rnorm(nrow(s02), mean = 0, sd = 2)
s02$Group <- factor(s02$group)散布図とグループを区別しない回帰直線
ggplot(s02, aes(x = factor, y = outcome)) +
geom_point(alpha = .4,
color = "#7CAE00") +
geom_smooth(method = "lm", se = F, linetype = "dotted") +
labs(title = "ランダム切片模擬データ（全体）",
subtitle = sprintf("（切片： %.3f，傾き %.3f）",
lm(outcome ~ factor, s02)$coefficients["(Intercept)"],
lm(outcome ~ factor, s02)$coefficients["factor"]),
caption = "※ 10グループ，各30人",
x = "説明変数",
y = "結果変数") +
theme_classic()
散布図とグループごとの回帰直線
ggplot(s02, aes(x = factor, y = outcome, color = Group)) +
geom_point(alpha = .4,
color = "#7CAE00") +
geom_smooth(method = "lm", se = F, linetype = "dotted") +
labs(title = "ランダム切片模擬データ（回帰直線は切片・傾き共にグループ別）",
# subtitle = sprintf("（切片： %.3f，傾き %.3f）",
# lm(outcome ~ factor, s02)$coefficients["(Intercept)"],
# lm(outcome ~ factor, s02)$coefficients["factor"]),
caption = "※ 10グループ，各30人",
x = "説明変数",
y = "結果変数") +
theme_classic()
切片のみグループ毎の特定の値（固定効果）を取り，傾きは共通であると仮定した一般線形モデルを当てはめると，
summary(reg_s02 <- lm(outcome ~ Group + factor, s02))
Call:
lm(formula = outcome ~ Group + factor, data = s02)
Residuals:
Min 1Q Median 3Q Max
-5.0892 -1.2162 0.0538 1.2120 4.8224
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.15873 0.37977 32.016 < 2e-16 ***
Group2 2.54024 0.46735 5.435 1.16e-07 ***
Group3 4.38422 0.46735 9.381 < 2e-16 ***
Group4 5.96637 0.46735 12.766 < 2e-16 ***
Group5 7.88572 0.46735 16.873 < 2e-16 ***
Group6 9.93238 0.46735 21.253 < 2e-16 ***
Group7 12.27570 0.46735 26.267 < 2e-16 ***
Group8 13.88156 0.46735 29.703 < 2e-16 ***
Group9 16.09770 0.46735 34.445 < 2e-16 ***
Group10 17.83086 0.46735 38.153 < 2e-16 ***
factor 0.48620 0.01207 40.270 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.81 on 289 degrees of freedom
Multiple R-squared: 0.9401, Adjusted R-squared: 0.938
F-statistic: 453.5 on 10 and 289 DF, p-value: < 2.2e-16anova(reg_s02)- factorの係数は，真値は .5 で，推定値はほぼ0.486なので，まずまず正しく推定できている。
- 切片は，真値は12, 14, 16, 18, 20, 22, 24, 26, 28, 30であるが，推定値は12.16, 14.7, 16.54, 18.13, 20.04, 22.09, 24.43, 26.04, 28.26, 29.99である。
10の切片の要約統計量は次の通り。
summary(intercepts); sd(intercepts); var(intercepts) Min. 1st Qu. Median Mean 3rd Qu. Max.
12.16 16.94 21.07 21.24 25.64 29.99 [1] 5.946108[1] 35.35619このデータを，マルチレヴェル分析で切片にのみ変量効果を仮定して計算してみる。
sum_lme01 <- summary(lme01 <- lme4::lmer(outcome ~ 1 + factor + (1|group), s02))
sum_lme01Linear mixed model fit by REML ['lmerMod']
Formula: outcome ~ 1 + factor + (1 | group)
Data: s02
REML criterion at convergence: 1267.1
Scaled residuals:
Min 1Q Median 3Q Max
-2.80349 -0.67757 0.03319 0.66513 2.67921
Random effects:
Groups Name Variance Std.Dev.
group (Intercept) 35.247 5.937
Residual 3.276 1.810
Number of obs: 300, groups: group, 10
Fixed effects:
Estimate Std. Error t value
(Intercept) 21.23820 1.88961 11.24
factor 0.48620 0.01207 40.27
Correlation of Fixed Effects:
(Intr)
factor -0.099- factorの係数0.4861985は上の固定効果モデル0.4861985と一致している。
- ランダム切片の平均21.2382018も固定効果の切片の平均21.2382018と一致している。
- ランダム切片の標準偏差35.2469889958483は，固定効果の切片の分散や標準偏差35.3561946とかなり近い。
- 切片のみにグループ差を仮定する固定効果モデルでは，グループの数だけ切片が具体的な定数として計算される（上の模擬データでは10グループ）。ランダム切片モデルでは，グループ数が幾つであろうと，切片の平均と分散の二つしか計算しない。「グループによってどの程度異なるか」にしか関心が無ければ，全ての値を具体的に計算しなくても平均と分散で十分かもしれないし，グループ数が増えるほどこの事があてはまる。
ランダム切片・ランダム傾きモデル
独立変数は既に集団平均でセンタリングしてあるとすると， \[\begin{align*} y_{ij} =& \color{green}\beta_{0j}\color{black}+\color{green}\beta_{1j}\color{black}x_{ij}+\color{red}r_{ij} \\ &\color{green}\beta_{0j}\color{black} = \color{blue}\gamma_{00}\color{black}+\color{red}u_{0j} \\ &\color{green}\beta_{ij}\color{black} = \color{blue}\gamma_{10}\color{black}+\color{red}u_{1j} \\ y_{ij} =& (\color{blue}\gamma_{00}\color{black}+\color{red}u_{0j}\color{black})+(\color{blue}\gamma_{10}\color{black}+\color{red}u_{1j}\color{black})x_{ij}+\color{red}r_{ij} \\ =& \color{blue}\gamma_{00}\color{black}+\color{blue}\gamma_{10}\color{black}x_{ij}+\color{red}u_{1j}\color{black}x_{ij}+\color{red}u_{0j}\color{black}+\color{red}r_{ij}\\ \end{align*}\]
\[\begin{align*} \color{red}{r_{ij}} &\sim\;N(0,\sigma^2)\\ \left( \begin{array}{c} \color{red}{u_{0j}} \\ \color{red}{u_{1j}} \\ \end{array} \right) &\sim\;N\left(0, \left[ \begin{array}{cc} \tau_{00} & \tau_{01}\\ \tau_{10} & \tau_{11}\\ \end{array} \right] \right) \end{align*}\]
切片と傾きを集団特性で説明
切片が，属する集団\(j\)の特性\(w_{j}\)によって変化するとし，
更に傾きも\(w_{j}\)によって変化すると考えると，
\[\begin{align*} y_{ij} =& \color{green}\beta_{0j}\color{black}+\color{green}\beta_{1j}\color{black}x_{ij}+\color{red}r_{ij} \\ &\color{green}\beta_{0j}\color{black} = \color{blue}\gamma_{00}\color{black}+\color{blue}\gamma_{01}\color{black}w_{j}+\color{red}u_{0j} \\ &\color{green}\beta_{ij}\color{black} = \color{blue}\gamma_{10}\color{black}+\color{blue}\gamma_{11}\color{black}w_{j}+\color{red}u_{1j} \\ y_{ij} =& (\color{blue}\gamma_{00}\color{black}+\color{blue}\gamma_{01}\color{black}w_{j}+\color{red}u_{0j}\color{black})+(\color{blue}\gamma_{10}\color{black}+\color{blue}\gamma_{11}\color{black}w_{j}+\color{red}u_{1j}\color{black})x_{ij}+\color{red}r_{ij} \\ =& \color{blue}\gamma_{00}\color{black}+\color{blue}\gamma_{10}\color{black}x_{ij}+\color{blue}\gamma_{01}\color{black}w_{j}+\color{blue}\gamma_{11}\color{black}w_{j}x_{ij}+\color{red}u_{1j}\color{black}x_{ij}+\color{red}u_{0j}\color{black}+\color{red}r_{ij}\\ \end{align*}\]
\[\begin{align*} \color{red}{r_{ij}} &\sim\;N(0,\sigma^2)\\ \left( \begin{array}{c} \color{red}{u_{0j}} \\ \color{red}{u_{1j}} \\ \end{array} \right) &\sim\;N\left(0, \left[ \begin{array}{cc} \tau_{00} & \tau_{01}\\ \tau_{10} & \tau_{11}\\ \end{array} \right] \right) \end{align*}\]
モデルを作る時はまずは最初の展開前の式を考えて，統計ソフトで指定する際に展開後の式を与える。
2 イヴェントヒストリー分析
2-1 パーソンピリオドデータ
2-2 離散時間ロジスティック回帰モデル
発展1 マルチレヴェルロジットモデル
発展2 統計的因果推論
Footnotes
教育年数はかなり離散的な値を取るので通常は平均と等しいケースは無いだろう。それが気になるならば，例えば各ケースのeduから12を引いてセンタリングすると，高卒の人を基準とする事も出来る。↩︎