マルチレベル型のデータは、集団と個人レベルの情報からなるが、これは家族のような集団にも応用できる場合がある。例えば成人子を持つ親に、自分自身と自分の子供の最終学歴をたずねたデータがあり、子供の教育年数を親の教育年数で予測したいとしよう。すべての対象者に関して成人子が一人以下ならばデータのレベルは一つだけであり、OLS で推定することができる。しかし、成人子が二人以上いる家族の情報がある程度存在する場合、データのレベルは二つあり、マルチレベル分析するのが適切である。
以下はガルトンの有名な親子の身長のデータであり、子の身長を親の身長で予測する場合、OLS ではなく、マルチレベル分析するのが適切である。
library(mosaicData)
head(Galton)
## family father mother sex height nkids
## 1 1 78.5 67.0 M 73.2 4
## 2 1 78.5 67.0 F 69.2 4
## 3 1 78.5 67.0 F 69.0 4
## 4 1 78.5 67.0 F 69.0 4
## 5 2 75.5 66.5 M 73.5 4
## 6 2 75.5 66.5 M 72.5 4
d0 <- Galton[, c(2,3,5)]
plot(d0)
round(cor(d0), 2)
## father mother height
## father 1.00 0.07 0.28
## mother 0.07 1.00 0.20
## height 0.28 0.20 1.00
library(lme4)
m1 <- lmer(height ~ father + mother + sex + nkids + (1 | family),
data = Galton)
summary(m1)
## Linear mixed model fit by REML ['lmerMod']
## Formula: height ~ father + mother + sex + nkids + (1 | family)
## Data: Galton
##
## REML criterion at convergence: 3897.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2059 -0.6031 0.0134 0.6195 3.7692
##
## Random effects:
## Groups Name Variance Std.Dev.
## family (Intercept) 0.8278 0.9099
## Residual 3.8406 1.9597
## Number of obs: 898, groups: family, 197
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 17.14566 3.66248 4.68
## father 0.39552 0.03846 10.28
## mother 0.30838 0.04206 7.33
## sexM 5.21775 0.13815 37.77
## nkids -0.03795 0.03721 -1.02
##
## Correlation of Fixed Effects:
## (Intr) father mother sexM
## father -0.678
## mother -0.677 -0.078
## sexM -0.071 0.040 0.026
## nkids -0.110 0.118 -0.043 0.044
少しだけ父親のほうが母親よりも相関係数もマルチレベル分析の回帰係数も大きい。父親と母親の係数が等しいかどうかを Wald 検定すると、以下のようになる。
library(aod)
wald.test(vcov(m1), fixef(m1), L = rbind(c(0, 1, -1, 0, 0)))
## Wald test:
## ----------
##
## Chi-squared test:
## X2 = 2.2, df = 1, P(> X2) = 0.14
父母の回帰係数に有意な差はない。
また、個人をレベル 2 、質問項目のタイプをレベル 1 として分析することもできる(と言われている)。例えば、仕事の内容、賃金、労働条件についての満足度を労働者の性別や年齢で予測したいとしよう。もちろん、別々に回帰分析してもいいが、例えば性別の効果が、内容、賃金、労働条件で有意に異なるかどうか検定したいという場合、マルチレベル分析を用いることができる(構造方程式モデルでも同じことができる)。
library(faraway)
head(hsb)
## id gender race ses schtyp prog read write math science socst
## 1 70 male white low public general 57 52 41 47 57
## 2 121 female white middle public vocation 68 59 53 63 61
## 3 86 male white high public general 44 33 54 58 31
## 4 141 male white high public vocation 63 44 47 53 56
## 5 172 male white middle public academic 47 52 57 53 61
## 6 113 male white middle public academic 44 52 51 63 61
d1 <- hsb[, c(2:3, 7:11)]
head(d1)
## gender race read write math science socst
## 1 male white 57 52 41 47 57
## 2 female white 68 59 53 63 61
## 3 male white 44 33 54 58 31
## 4 male white 63 44 47 53 56
## 5 male white 47 52 57 53 61
## 6 male white 44 52 51 63 61
d1.1 <- reshape(d1, varying = list(names(d1)[3:7]), direction = "long")
head(d1.1)
## gender race time read id
## 1.1 male white 1 57 1
## 2.1 female white 1 68 2
## 3.1 male white 1 44 3
## 4.1 male white 1 63 4
## 5.1 male white 1 47 5
## 6.1 male white 1 44 6
d1.1 <- d1.1[order(d1.1$id), ]
names(d1.1)[3:4] <- c("subject", "score")
d1.1$subject <- factor(d1.1$subject, labels = names(hsb)[7:11])
head(d1.1)
## gender race subject score id
## 1.1 male white read 57 1
## 1.2 male white write 52 1
## 1.3 male white math 41 1
## 1.4 male white science 47 1
## 1.5 male white socst 57 1
## 2.1 female white read 68 2
library(lme4)
m1 <- lmer(score ~ subject + gender + race + (1 | id), data = d1.1)
m2 <- update(m1, ~ . + subject : (gender + race))
library(texreg)
screenreg(list(m1, m2))
##
## =======================================================
## Model 1 Model 2
## -------------------------------------------------------
## (Intercept) 46.78 *** 46.44 ***
## (1.81) (2.17)
## subjectwrite 0.54 3.41
## (0.64) (2.00)
## subjectmath 0.41 0.08
## (0.64) (2.00)
## subjectscience -0.38 -4.46 *
## (0.64) (2.00)
## subjectsocst 0.17 3.42
## (0.64) (2.00)
## gendermale -0.37 1.03
## (1.10) (1.35)
## raceasian 7.10 * 5.19
## (2.88) (3.55)
## racehispanic 0.01 -0.33
## (2.33) (2.88)
## racewhite 7.21 *** 7.00 **
## (1.84) (2.26)
## subjectwrite:gendermale -5.74 ***
## (1.25)
## subjectmath:gendermale -0.36
## (1.25)
## subjectscience:gendermale 1.32
## (1.25)
## subjectsocst:gendermale -2.21
## (1.25)
## subjectwrite:raceasian 4.25
## (3.28)
## subjectmath:raceasian 5.39
## (3.28)
## subjectscience:raceasian 3.65
## (3.28)
## subjectsocst:raceasian -3.73
## (3.28)
## subjectwrite:racehispanic -0.51
## (2.66)
## subjectmath:racehispanic 0.87
## (2.66)
## subjectscience:racehispanic 2.46
## (2.66)
## subjectsocst:racehispanic -1.10
## (2.66)
## subjectwrite:racewhite -0.59
## (2.09)
## subjectmath:racewhite 0.14
## (2.09)
## subjectscience:racewhite 4.12 *
## (2.09)
## subjectsocst:racewhite -2.63
## (2.09)
## -------------------------------------------------------
## AIC 6935.42 6858.56
## BIC 6989.41 6991.07
## Log Likelihood -3456.71 -3402.28
## Num. obs. 1000 1000
## Num. groups: id 200 200
## Var: id (Intercept) 50.88 51.35
## Var: Residual 40.54 38.21
## =======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
ちなみにマルチレベル分析では従属変数は必ず個人レベル(レベル 1)の変数でなければならないので注意。
マルチレベル分析が最も応用されているのは、パネル・データの分析だが、パネル・データ分析では、成長曲線モデル、ランダム効果モデル、ハイブリッド・モデルなど別の名前で呼ばれていることが多い。しかし、実質的には同じものを指していっている。
パネル・データでは、個人ごとにある変数の変化を見ることができる。それは収入でも体重でも仕事に関する価値観でもよいが、これらが個人の属性によってどのように異なるのかをマルチレベル・モデルで予測できる。このようなモデルは成長曲線モデルとも言われる。
パネル・データはマルチレベル構造をしているが、個人がレベル 2 で時点がレベル 1 と見なされる。時点 \(t\) における個人 \(i\) の従属変数の値を \(Y_{ti}\) とすると、切片のみのランダム切片モデルは、 \[ \begin{array}{lll} Y_{ti} & = \ \beta_{0i} + R_{ti} & \\ \beta_{0i} & = \gamma_{00} + U_{0i} \end{array} \] である。このモデルでは従属変数の時間とともに変化しない個人差は \(U_{0i}\) で表現され、時間による変動は \(R_{ti}\) で表現される。
成長曲線モデルにおいては、時間とともに変化しない変数(time constant variable, 時定変数と呼んでおく)と時間とともに変化する変数(time varying variable, 時変変数とよんでおく)。の区別が重要である。時定変数は性別や出生コーホートのように時間とともには変化しない変数で、レベル 2 の変数である。時変変数は、年齢や収入のように時間とともに変化する変数であり、レベル 1 に属する。
あとは通常のマルチレベル・モデルと同じように、各自の関心や仮説に応じてモデルを組んで行けばよく、とりたてて通常のマルチレベル・モデルと違いがあるわけではない。とはいえ、いくつか考えられるモデルを、以下の架空のデータを例に解説していこう。この data2 というデータフレームは 100人の人に関して 5年間、毎年 BMI (Body Mass Index, 肥満の簡易指標) を調べた結果であり、性別、教育年数、年齢も調べられている。
head(data2, 12)
## id.person sex edu time age bmi
## 1.1 1 male 14 1 59 21.54290
## 1.2 1 male 14 2 60 21.19478
## 1.3 1 male 14 3 61 25.48479
## 1.4 1 male 14 4 62 23.23149
## 1.5 1 male 14 5 63 26.97935
## 2.1 2 female 11 1 41 23.16405
## 2.2 2 female 11 2 42 19.76011
## 2.3 2 female 11 3 43 23.29545
## 2.4 2 female 11 4 44 23.70071
## 2.5 2 female 11 5 45 22.73477
## 3.1 3 male 10 1 32 24.72924
## 3.2 3 male 10 2 33 18.22938
このデータは、1行がある個人のある年の情報に対応している。こういったデータをパーソン・イヤー・データということがある。また、ロング形式のパネル・データとも言う。成長曲線モデルをマルチレベル分析のソフトウェアで推定するには、まずデータをこのようなロング形式で用意する必要がある。
年齢とともに BMI は変化すると考えられるが、このような変化の軌跡が、性別や学歴によってどう異なるのかを調べたいとしよう。既存の研究から BMI は年齢とともに上昇するが、その後減少に転じることがわかっているとしよう。その場合、検討すべきモデルは、独立変数として年齢と年齢の二乗を投入したモデルであろう。すなわち、個人 \(i\) の時点 \(t\) における年齢を \(X_{ti}\) とすると、 \[ \begin{array}{lll} Y_{ti} & = \ \beta_{0i} + \beta_1 X_{ti} + \beta_2 X_{ti}^2 + R_{ti} & \\ \beta_{0i} & = \gamma_{00} + U_{0i} \end{array} \] であろう。これは切片に関しては個人差があるものの、年齢による変化はすべての人に関して同じであると仮定したモデルである。切片のみのモデルとこの年齢のみを投入したモデルを推定してみよう。
library(lme4)
lg0 <- lmer(bmi ~ 1 + (1 | id.person), data = data2)
ci0 <- confint(lg0, oldNames = FALSE)
ci0
## 2.5 % 97.5 %
## sd_(Intercept)|id.person 1.755649 2.600367
## sigma 2.799787 3.216224
## (Intercept) 21.512429 22.512951
lg1 <- lmer(bmi ~ age + I(age^2) + (1 | id.person), data = data2)
library(texreg)
screenreg(list(lg0, lg1), custom.model.names = c("Model 0", "Model 1"))
##
## ======================================================
## Model 0 Model 1
## ------------------------------------------------------
## (Intercept) 22.01 *** 8.72 ***
## (0.25) (1.76)
## age 0.47 ***
## (0.08)
## I(age^2) -0.00 ***
## (0.00)
## ------------------------------------------------------
## AIC 2650.09 2567.71
## BIC 2662.73 2588.79
## Log Likelihood -1322.04 -1278.86
## Num. obs. 500 500
## Num. groups: id.person 100 100
## Var: id.person (Intercept) 4.66 0.51
## Var: Residual 8.98 8.94
## ======================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
切片のみのモデルでランダム効果の 95% 信頼区間の推定結果を見ると、グループ・レベル(つまりこの場合は個人レベル)の標準偏差は 1.8~2.6, 時点レベルの標準偏差が 2.8~3.2 なので、個人間の分散も個人内の分散もそれなりにあることがわかる。次に年齢と年齢の二乗を投入した Model 1 の推定結果を見ると、確かに曲線的な関係が見られる。次にこのような年齢による変化に個人差があるかどうかを見るために、ランダム係数モデルを推定する。すなわち、 \[ \begin{array}{lll} Y_{ti} & = \ \beta_{0i} + \beta_{1i} X_{ti} + \beta_2 X_{ti}^2 + R_{it} & \\ \beta_{0i} & = \gamma_{00} + U_{0i} & \\ \beta_{1i} & = \gamma_{10} + U_{1i} \end{array} \] というモデルを推定してみる。
data2$age10 <- (data2$age - 25)/10 # 年齢を25歳でセンタリングして10歳単位に。
lg2 <- lmer(bmi ~ age10 + I(age10^2) + (age10 | id.person), data = data2)
anova(lg1, lg2, refit=FALSE) # さきほどのように区間推定してもよい
## Data: data2
## Models:
## lg1: bmi ~ age + I(age^2) + (1 | id.person)
## lg2: bmi ~ age10 + I(age10^2) + (age10 | id.person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## lg1 5 2567.7 2588.8 -1278.9 2557.7
## lg2 7 2553.8 2583.3 -1269.9 2539.8 17.922 2 0.0001283 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
screenreg(list(lg1, lg2), custom.coef.names = c(NA, NA, NA, "age", "I(age^2)"))
##
## ============================================================
## Model 1 Model 2
## ------------------------------------------------------------
## (Intercept) 8.72 *** 18.21 ***
## (1.76) (0.32)
## age 0.47 *** 2.96 ***
## (0.08) (0.37)
## I(age^2) -0.00 *** -0.35 ***
## (0.00) (0.09)
## ------------------------------------------------------------
## AIC 2567.71 2553.79
## BIC 2588.79 2583.29
## Log Likelihood -1278.86 -1269.89
## Num. obs. 500 500
## Num. groups: id.person 100 100
## Var: id.person (Intercept) 0.51 0.00
## Var: Residual 8.94 8.83
## Var: id.person age10 0.11
## Cov: id.person (Intercept) age10 0.00
## ============================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
確かに年齢の効果の分散は有意である。年齢を25歳でセンタリングしているのに深い意味は無いが、センタリングは必要である。年齢を 10歳単位に変換しているのは、独立変数のスケールが違いすぎると lmer() が警告を発してくるからである。
次に年齢の二乗の効果に関しても個人差が無いか調べてみよう。推定するモデルは \[ \begin{array}{lll} Y_{ti} & = \ \beta_{0i} + \beta_{1i} X_{ti} + \beta_2 X_{ti}^2 + R_{it} & \\ \beta_{0i} & = \gamma_{00} + U_{0i} & \\ \beta_{1i} & = \gamma_{10} + U_{1i} & \\ \beta_{2i} & = \gamma_{20} + U_{2i} & \end{array} \] である。以下がこのモデルの推定。
lg3 <- lmer(bmi ~ age10 + I(age10^2) + (age10 + I(age10^2) | id.person), data = data2)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge: degenerate Hessian with 1 negative
## eigenvalues
anova(lg2, lg3, refit=FALSE) # さきほどのように区間推定してもよい
## Data: data2
## Models:
## lg2: bmi ~ age10 + I(age10^2) + (age10 | id.person)
## lg3: bmi ~ age10 + I(age10^2) + (age10 + I(age10^2) | id.person)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## lg2 7 2553.8 2583.3 -1269.9 2539.8
## lg3 10 2559.7 2601.9 -1269.9 2539.7 0.0787 3 0.9943
lmer() が警告を発しており、推定結果が収束しなかったようである。尤度比検定の結果も有意ではないので、年齢の二乗の傾きにかんしてはランダム効果を仮定しないほうがよさそうである。
次に個人レベルの変数として性別 (\(Z_{1i}\)) と教育年数 (\(Z_{2i}\))を投入しよう。まずは切片のみ。推定するモデルは、以下のとおり。 \[ \begin{array}{lll} Y_{ti} & = \ \beta_{0i} + \beta_{1i} X_{ti} + \beta_2 X_{ti}^2 + R_{it} & \\ \beta_{0i} & = \gamma_{00} + Z_{1i} + Z_{2i}+ U_{0i} & \\ \beta_{1i} & = \gamma_{10} + U_{1i} \end{array} \]
data2$edu.c <- data2$edu - mean(data2$edu)
lg4 <- lmer(bmi ~ age10 + I(age10^2) + sex + edu.c + (age10 | id.person), data = data2)
screenreg(list(lg2, lg4), custom.model.names = c("Model 2", "Model 3"))
##
## ============================================================
## Model 2 Model 3
## ------------------------------------------------------------
## (Intercept) 18.21 *** 18.36 ***
## (0.32) (0.33)
## age10 2.96 *** 3.04 ***
## (0.37) (0.34)
## I(age10^2) -0.35 *** -0.36 ***
## (0.09) (0.07)
## sexfemale -0.55 *
## (0.27)
## edu.c -0.37 ***
## (0.05)
## ------------------------------------------------------------
## AIC 2553.79 2520.40
## BIC 2583.29 2558.33
## Log Likelihood -1269.89 -1251.20
## Num. obs. 500 500
## Num. groups: id.person 100 100
## Var: id.person (Intercept) 0.00 0.00
## Var: id.person age10 0.11 0.00
## Cov: id.person (Intercept) age10 0.00 -0.00
## Var: Residual 8.83 8.59
## ============================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05
交差レベル交互作用効果も推定してみよう。
lg5 <- lmer(bmi ~ I(age10^2) + age10 * (sex + edu.c) + (age10 | id.person), data = data2)
summary(lg5)
## Linear mixed model fit by REML ['lmerMod']
## Formula: bmi ~ I(age10^2) + age10 * (sex + edu.c) + (age10 | id.person)
## Data: data2
##
## REML criterion at convergence: 2504.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8720 -0.7140 -0.0272 0.6076 3.4500
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id.person (Intercept) 7.946e-14 2.819e-07
## age10 3.749e-14 1.936e-07 -1.00
## Residual 8.560e+00 2.926e+00
## Number of obs: 500, groups: id.person, 100
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 18.19850 0.36964 49.23
## I(age10^2) -0.36264 0.07473 -4.85
## age10 3.16601 0.34952 9.06
## sexfemale -0.19597 0.47131 -0.42
## edu.c -0.23598 0.09439 -2.50
## age10:sexfemale -0.21550 0.20522 -1.05
## age10:edu.c -0.07328 0.04064 -1.80
##
## Correlation of Fixed Effects:
## (Intr) I(10^2 age10 sexfml edu.c ag10:s
## I(age10^2) 0.529
## age10 -0.754 -0.919
## sexfemale -0.549 0.043 0.176
## edu.c -0.066 -0.019 0.058 0.140
## age10:sxfml 0.473 0.013 -0.287 -0.821 -0.164
## age10:edu.c 0.097 0.032 -0.110 -0.159 -0.829 0.268
最後にハイブリッド・モデルと呼ばれるモデルを推定してみよう。これは時変変数をグループ平均センタリングして、グループ平均も時定変数として投入するようなモデルである。 \[ \begin{array}{lll} Y_{ti} & = \ \beta_{0i} + \beta_{1i} (X_{ti} - \bar X_{ti}) + \beta_2 X_{ti}^2 + R_{it} & \\ \beta_{0i} & = \gamma_{00} + Z_{1i} + Z_{2i} + \bar X_{\cdot i} + \bar X_{\cdot i}^2 + U_{0i} & \\ \beta_{1i} & = \gamma_{10} + Z_{1i} + Z_{2i} + \bar X_{\cdot i} + U_{1i} \end{array} \] こうすると、グループ平均センタリングしたほうの効果は、個人の時間とともに変化しないすべての属性を統制した効果として解釈でき、いっぽうグループ平均のほうは、個人間の変動に対する効果とみなせる。
この場合、対象者の調査期間における平均年齢をモデルに投入しているわけであるから、要するに出生コーホートの効果を検討しているのと同じである。
data2$age10.m <- ave(data2$age10, data2$id.person)
data2$age10.gc <- data2$age10 - data2$age10.m # グループ平均センタリング
data2$age10.m <- data2$age10.m - mean(data2$age10.m) # 全体平均センタリング
lg6 <- lmer(bmi ~ I(age10.gc^2) + age10.gc * (sex + edu.c + age10.m ) + I(age10.m^2) + (age10 | id.person), data = data2)
summary(lg6)
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## bmi ~ I(age10.gc^2) + age10.gc * (sex + edu.c + age10.m) + I(age10.m^2) +
## (age10 | id.person)
## Data: data2
##
## REML criterion at convergence: 2481.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9135 -0.6623 -0.0498 0.6269 3.4190
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## id.person (Intercept) 0.000e+00 0.000000
## age10 8.725e-05 0.009341 NaN
## Residual 8.520e+00 2.918938
## Number of obs: 500, groups: id.person, 100
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 23.17760 0.28431 81.52
## I(age10.gc^2) -12.55417 7.80119 -1.61
## age10.gc 1.22686 1.31715 0.93
## sexfemale -0.55650 0.26699 -2.08
## edu.c -0.37400 0.05233 -7.15
## age10.m 1.65503 0.10294 16.08
## I(age10.m^2) -0.35951 0.07635 -4.71
## age10.gc:sexfemale -0.55685 1.87920 -0.30
## age10.gc:edu.c -0.89659 0.36997 -2.42
## age10.gc:age10.m -0.58056 0.69656 -0.83
##
## Correlation of Fixed Effects:
## (Intr) I(g10.g^2) ag10.g sexfml edu.c ag10.m I(g10.m^2)
## I(g10.gc^2) -0.549
## age10.gc 0.000 0.000
## sexfemale -0.514 0.000 0.000
## edu.c -0.091 0.000 0.000 0.181
## age10.m 0.177 0.000 0.000 -0.085 -0.066
## I(ag10.m^2) -0.519 0.000 0.000 0.095 0.013 -0.289
## ag10.gc:sxf 0.000 0.000 -0.713 0.000 0.000 0.000 0.000
## ag10.gc:d.c 0.000 0.000 -0.129 0.000 0.000 0.000 0.000
## ag10.gc:10. 0.000 0.000 0.043 0.000 0.000 0.000 0.000
## ag10.: a10.:.
## I(g10.gc^2)
## age10.gc
## sexfemale
## edu.c
## age10.m
## I(ag10.m^2)
## ag10.gc:sxf
## ag10.gc:d.c 0.180
## ag10.gc:10. -0.060 -0.065
## convergence code: 0
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
この結果を見ると、個人内の年齢の効果と個人間の年齢の効果には違いがあるようである。
親子関係の分析をするような場合、父親や母親に関する情報と、子供に関する情報を得て、例えば、父の価値観と子の価値観の相関を分析する場合がある。仮に世帯をサンプリングし、世帯成員すべての価値観を調べたとすると、子供が複数いる世帯はそれぞれの子供の価値観が得られる。これもマルチレベル型のデータで、世帯がグループで子供が個人レベルの単位、父の価値観はグループレベルの変数、子供の価値観は個人レベルの変数となる。こういった分析はときどき必要になる。
今まではレベルが2つのデータをあつかってきたが、レベルは 3 つあっても 4 つあっても構わない。例えば、レベル 1 が個人、レベル 2 が市区町村、レベル 3 が都道府県といった構造のデータが考えられる。ランダム効果に関してはすでに述べたような仮定が置かれるが、モデルが複雑になるほど推定に時間がかかったり、収束しなくなったりしやすくなるので、注意。
R で 3レベルのランダム切片モデルを推定する場合、
lmer(固定効果のモデル式 + (1 | レベル2 のグループID) + (1 | レベル3 のグループID), …)
という指定が必要。
3レベルモデルの変形に交差分類モデルがある。通常の3レベルモデルは、レベル1のデータがレベル2のグループに分類され、レベル2のグループがレベル3のグループにさらに分類される、という構造を持つ。しかし、そのような階層関係が成り立たない場合もある。例えば、企業とコーホート(出生年)をグループ・レベルに設定したい場合、同じ出生年の人が複数の企業に存在するので、二種類のグループに階層関係が成り立たない場合もある。
このようなモデルを交差分類モデルと呼ぶ。R の場合、交差分類モデルも 3レベルモデルと同じ仮定のもとに同じように分析できる。
被説明変数が、二値変数だったり、順序をあらわすカテゴリだったり、度数(0位上の整数)だったりする場合、いわゆる一般化線形モデル (Generalized Linear Model: GLM) が用いられる。GLM に関してもマルチレベル分析が可能である。R で推定する場合、lme4パッケージの glmer() 関数を使うのが簡単である。二項ロジスティック回帰分析なら
glmer(… , family = binomial)
… の部分は lmer() 関数の場合と同じように引数を指定すればよい。ポアソン回帰(や対数線形モデル)なら、
glmer(…, family = poisson)
である。順序ロジスティック回帰分析の場合、ordinal パッケージの clmm() 関数が使える。
clmm(lmerと同じ引数)
最後に負二項回帰分析は lme4パッケージの glmer.nb() で計算できるが、開発途中のものなので、問題があるかもしれないらしい。
glmer.nb(lmerと同じ引数)
以下は計算例。
レベル1のランダム効果に系列相関を仮定したり、分散の不均一性を仮定したりする必要がある場合、nlme パッケージの lme() 関数が使える。また、