マルチレベルモデリングChapter4:モデルの比較方法
UPDATE: 2022-11-16 21:55:29
はじめに
縦断データの分析I―変化についてのマルチレベルモデリング―を利用してマルチレベルモデリングの勉強内容をまとめたもの。下記のサポートサイトにサンプルデータなどが保存されている。
分析内容について
分析内容は、前回に引き続き、青年期のアルコール摂取量の変化に関する分析。縦断データの分析I―変化についてのマルチレベルモデリング―のp75に詳細は書いてあるとおり、青年期のアルコール摂取量の変化に関する分析をここでもお借りする。82人(n)を対象に14歳,15歳,16歳の3回(age)の計測タイミングがある。alcuseはアルコール摂取に関する頻度の合成スコアで、予測変数として、友達の飲酒割合に関するpeerと親がアルコール依存症かを表すcoaが用意されている。
この分析の目的は、親がアルコール依存症であったり、自分の周囲の友だちが飲酒していれば、自分のアルコール摂取量も14歳,15歳,16歳の時間経過とともに、アルコール摂取量が増加していくのではないか、という仮説を検証すること。
library(tidyverse)
library(nlme)
library(DT)
alcohol1 <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/alcohol1_pp.txt", header=T, sep=",")
alcohol1 <- alcohol1 %>% mutate(
mean_peer = mean(peer),
flg_peer = ifelse(peer < mean_peer, "low peer", "high peer"),
flg_peer = factor(flg_peer, levels = c("low peer", "high peer"))
)
#Model A B
model.a <- lme(alcuse ~ 1, alcohol1, random = ~1 |id, method = "ML")
model.b <- lme(alcuse ~ age_14 , data = alcohol1, random= ~ age_14 | id, method = "ML")
datatable(alcohol1 %>% mutate_if(is.numeric, round, digit = 2))モデルC
この分析の主たる目的は、親のアルコール依存症(COA)の影響であり、周囲の飲酒環境(PEER)は統制変数である。このモデルは初期値と変化率の予測変数として、親のアルコール依存症(COA)が含まれているモデル。
\[ \begin{eqnarray} Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} COA_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{11} COA_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{10} TIME_{ij} + \gamma_{11} COA_{i} × TIME_{ij} + (\epsilon_{ij} + \zeta_{0i} + \zeta_{1i} × TIME_{ij}) \end{eqnarray} \]
model.c <- lme(
fixed = alcuse ~ coa*age_14,
random = ~ age_14 | id,
data = alcohol1,
method = "ML")
summary(model.c)## Linear mixed-effects model fit by maximum likelihood
## Data: alcohol1
## AIC BIC logLik
## 637.2026 665.2453 -310.6013
##
## Random effects:
## Formula: ~age_14 | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.6982690 (Intr)
## age_14 0.3880684 -0.219
## Residual 0.5807688
##
## Fixed effects: alcuse ~ coa * age_14
## Value Std.Error DF t-value p-value
## (Intercept) 0.3159517 0.13177099 162 2.397734 0.0176
## coa 0.7432120 0.19616697 80 3.788671 0.0003
## age_14 0.2929552 0.08492089 162 3.449742 0.0007
## coa:age_14 -0.0494299 0.12642140 162 -0.390993 0.6963
## Correlation:
## (Intr) coa age_14
## coa -0.672
## age_14 -0.460 0.309
## coa:age_14 0.309 -0.460 -0.672
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.5480534 -0.3879877 -0.1057547 0.3602391 2.3960739
##
## Number of Observations: 246
## Number of Groups: 82- 親がアルコール依存症ではない個人の初期値\(\gamma_{00}\)は0.32
- 親がアルコール依存症ではない個人とそうではない個人の初期値の差\(\gamma_{01}\)は0.74
- 親がアルコール依存症ではない個人の変化率\(\gamma_{10}\)は0.29
- 親がアルコール依存症ではない個人とそうではない個人の変化率の差\(\gamma_{11}\)は-0.05。ただ0の可能性がある。
モデルCからわかることは、親がアルコール依存症の個人は、そうではない個人よりも初期値は高いが、経年変化率については差がない、と解釈できる。また、\(\hat{\sigma}_{0}^{2}\)は前回のモデルBと比較して21％減少している。つまり、COAは有効な変数であったことがわかる。
list(
modelB = VarCorr(model.b),
modelC = VarCorr(model.c)
)## $modelB
## id = pdLogChol(age_14)
## Variance StdDev Corr
## (Intercept) 0.6243558 0.7901619 (Intr)
## age_14 0.1512021 0.3888472 -0.223
## Residual 0.3372926 0.5807690
##
## $modelC
## id = pdLogChol(age_14)
## Variance StdDev Corr
## (Intercept) 0.4875797 0.6982690 (Intr)
## age_14 0.1505971 0.3880684 -0.219
## Residual 0.3372924 0.5807688モデルD
このモデルは初期値と変化率の予測変数として、親のアルコール依存症(COA)と周囲の飲酒環境(PEER)が含まれているモデル。つまり、初期値と変化率へのCOAの影響を推定するものの、PEERで統制されているCOAの影響を推定している。
\[ \begin{eqnarray} Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{11} COA_{i} + \gamma_{12} PEER_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} PEER_{i} + \gamma_{10} TIME_{ij} + \gamma_{11} COA_{i} × TIME_{ij} + \gamma_{12} PEER_{i} × TIME_{ij} + (\epsilon_{ij} + \zeta_{0i} + \zeta_{1i} × TIME_{ij}) \end{eqnarray} \]
model.d <- lme(
fixed = alcuse ~ coa * age_14 + peer * age_14,
random = ~ age_14 | id,
data = alcohol1,
method = "ML")
summary(model.d)## Linear mixed-effects model fit by maximum likelihood
## Data: alcohol1
## AIC BIC logLik
## 608.6906 643.744 -294.3453
##
## Random effects:
## Formula: ~age_14 | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.4908199 (Intr)
## age_14 0.3729816 -0.033
## Residual 0.5807670
##
## Fixed effects: alcuse ~ coa * age_14 + peer * age_14
## Value Std.Error DF t-value p-value
## (Intercept) -0.3165142 0.14990060 161 -2.111494 0.0363
## coa 0.5791651 0.16450410 79 3.520673 0.0007
## age_14 0.4294286 0.11510206 161 3.730851 0.0003
## peer 0.6942956 0.11291813 79 6.148664 0.0000
## coa:age_14 -0.0140319 0.12631544 161 -0.111086 0.9117
## age_14:peer -0.1498150 0.08670485 161 -1.727873 0.0859
## Correlation:
## (Intr) coa age_14 peer c:g_14
## coa -0.371
## age_14 -0.436 0.162
## peer -0.686 -0.162 0.299
## coa:age_14 0.162 -0.436 -0.371 0.071
## age_14:peer 0.299 0.071 -0.686 -0.436 -0.162
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.59554954 -0.40004985 -0.07769103 0.46002869 2.29373877
##
## Number of Observations: 246
## Number of Groups: 82パラメタの係数解釈が少し複雑なので注意が必要。\(\hat{\gamma}_{00}\)、\(\hat{\gamma}_{10}\)は、親がアルコール依存症ではない個人の初期値と変化率ではなく、親がアルコール依存症ではない個人の一部、PEERの値も0の子どもの初期値と変化率。つまり、COA=0,PEER=0のときの初期値と変化率が\(\hat{\gamma}_{00}\)、\(\hat{\gamma}_{10}\)となる。\(\hat{\gamma}_{01}\)、\(\hat{\gamma}_{11}\)は、PEERで統制した後のCORの影響と解釈できる。
COA=0,PEER=0のときの初期値は負になることから、親がアルコール依存症ではなく、周囲の友人がアルコールを摂取しない個人は0ではない量のアルコールを摂取している、とは解釈できなくないが、alcuseは0以下を取らないので解釈には注意が必要。
モデルDは、PEERの効果を統制した上で、親がアルコール依存症ではない個人とそうではない個人の初期値の推定された差は0.579で、親がアルコール依存症ではない個人とそうではない個人の変化率の差は0(-0.014)かもしれない、と解釈できる。残差に関しても、\(\hat{\sigma}_{0}^{2}\)、\(\hat{\sigma}_{1}^{2}\)の両方が減少しており、PEERが有効な変数であることがわかる。
list(
modelC = VarCorr(model.c),
modelD = VarCorr(model.d)
)## $modelC
## id = pdLogChol(age_14)
## Variance StdDev Corr
## (Intercept) 0.4875797 0.6982690 (Intr)
## age_14 0.1505971 0.3880684 -0.219
## Residual 0.3372924 0.5807688
##
## $modelD
## id = pdLogChol(age_14)
## Variance StdDev Corr
## (Intercept) 0.2409042 0.4908199 (Intr)
## age_14 0.1391153 0.3729816 -0.033
## Residual 0.3372903 0.5807670モデルE
このモデルは初期値と変化率の予測変数として、親のアルコール依存症(COA)と周囲の飲酒環境(PEER)が含まれているモデル。ただ、初期値にはPEERとCOAが含まれる一方で、変化率にはPEERのみでCOAは含まれていない。
\[ \begin{eqnarray} Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{12} PEER_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} PEER_{i} + \gamma_{10} TIME_{ij} + \gamma_{12} PEER_{i} × TIME_{ij} + (\epsilon_{ij} + \zeta_{0i} + \zeta_{1i} × TIME_{ij}) \end{eqnarray} \]
model.e <- lme(
fixed = alcuse ~ coa + peer * age_14,
random = ~ age_14 | id,
data = alcohol1,
method = "ML")
summary(model.e)## Linear mixed-effects model fit by maximum likelihood
## Data: alcohol1
## AIC BIC logLik
## 606.7033 638.2513 -294.3516
##
## Random effects:
## Formula: ~age_14 | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.4908386 (Intr)
## age_14 0.3730384 -0.034
## Residual 0.5807690
##
## Fixed effects: alcuse ~ coa + peer * age_14
## Value Std.Error DF t-value p-value
## (Intercept) -0.3138215 0.14762318 162 -2.125828 0.0350
## coa 0.5711970 0.14773607 79 3.866334 0.0002
## peer 0.6951827 0.11240366 79 6.184698 0.0000
## age_14 0.4246867 0.10667949 162 3.980959 0.0001
## peer:age_14 -0.1513771 0.08538525 162 -1.772872 0.0781
## Correlation:
## (Intr) coa peer age_14
## coa -0.338
## peer -0.709 -0.146
## age_14 -0.410 0.000 0.351
## peer:age_14 0.334 0.000 -0.431 -0.814
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.59554431 -0.40414085 -0.08351916 0.45549511 2.29975165
##
## Number of Observations: 246
## Number of Groups: 82PEERの影響を統制したCOA、つまり、親がアルコール依存症ではない個人とそうではない個人の初期値は0.57で、COAの影響を統制したPEERの1ポイントの差は0.695高くなる。平均的な変化率は-0.151低くなる。
モデルEは、PEERの効果を統制した上で、親がアルコール依存症の個人は、そうではない個人よりも初期値が高いが、経年変化率については差がない、と解釈できる。残差に関しても、\(\hat{\sigma}_{0}^{2}, \hat{\sigma}_{1}^{2},
\hat{\sigma}_{\epsilon}^{2}\)は、あまり変化しておらず、変化率に対するCOAはあってもなくても同じということがわかる。
list(
modelD = VarCorr(model.d),
modelE = VarCorr(model.e)
)## $modelD
## id = pdLogChol(age_14)
## Variance StdDev Corr
## (Intercept) 0.2409042 0.4908199 (Intr)
## age_14 0.1391153 0.3729816 -0.033
## Residual 0.3372903 0.5807670
##
## $modelE
## id = pdLogChol(age_14)
## Variance StdDev Corr
## (Intercept) 0.2409226 0.4908386 (Intr)
## age_14 0.1391576 0.3730384 -0.034
## Residual 0.3372926 0.5807690モデルサマリー
| Dependent variable: | |||||||
| alcuse | |||||||
| A | B | C | D | E | F_c | G_c | |
| Intercept(\(\pi_{0i},\gamma_{00}\)) | 0.922*** | 0.651*** | 0.316** | -0.317** | -0.314** | 0.394*** | 0.651*** |
| (0.096) | (0.106) | (0.132) | (0.150) | (0.148) | (0.105) | (0.081) | |
| coa(\(\pi_{0i},\gamma_{01}\)) | 0.743*** | 0.579*** | 0.571*** | 0.571*** | 0.571*** | ||
| (0.196) | (0.165) | (0.148) | (0.148) | (0.148) | |||
| peer(\(\pi_{0i},\gamma_{02}\)) | 0.694*** | 0.695*** | 0.695*** | 0.695*** | |||
| (0.113) | (0.112) | (0.112) | (0.112) | ||||
| age_14(\(\pi_{1i},\gamma_{10}\)) | 0.271*** | 0.293*** | 0.429*** | 0.425*** | 0.271*** | 0.271*** | |
| (0.063) | (0.085) | (0.115) | (0.107) | (0.062) | (0.062) | ||
| coa:age_14(\(\pi_{1i},\gamma_{11}\)) | -0.049 | -0.014 | |||||
| (0.126) | (0.126) | ||||||
| peer:age_14(\(\pi_{1i},\gamma_{12}\)) | -0.150* | -0.151* | -0.151* | -0.151* | |||
| (0.087) | (0.085) | (0.085) | (0.085) | ||||
| Observations | 246 | 246 | 246 | 246 | 246 | 246 | 246 |
| Log Likelihood | -335 | -318 | -310 | -294 | -294 | -294 | -294 |
| Akaike Inf. Crit. | 676 | 648 | 637 | 608 | 606 | 606 | 606 |
| Bayesian Inf. Crit. | 686 | 669 | 665 | 643 | 638 | 638 | 638 |
| Note: | p<0.1; p<0.05; p<0.01 | ||||||
モデルの可視化
書籍に従って、モデルCを可視化しておく。前回のノートに可視化の方法をまとめているので細かい話はそちらを参照のこと。このモデルは、親がアルコール依存症でない平均的な個人は、初期値が0.316、変化率が0.293の変化の軌跡を持ち、親がアルコール依存症の平均的な個人は、初期値が1.059、変化率が0.244の変化の軌跡を持つ。
\[ \begin{eqnarray} Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} COA_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{11} COA_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{10} TIME_{ij} + \gamma_{11} COA_{i} × TIME_{ij} + (\epsilon_{ij} + \zeta_{0i} + \zeta_{1i} × TIME_{ij}) \\ &=& 0.315 + 0.743 COA_{i} + 0.292 TIME_{ij} -0.049 COA_{i} × TIME_{ij} \end{eqnarray} \]
初期値の高さに違いはあるものの、変化率に関してはあまり差はない。
fixef.c <- fixef(model.c)
df_fit_plt_c <-
expand_grid(age_14 = 0:2, coa = 0:1) %>%
arrange(coa) %>%
mutate(alcuse =
fixef.c[[1]] +
fixef.c[[2]] * coa +
fixef.c[[3]] * age_14 +
fixef.c[[4]] * age_14 * coa
)
ggplot(df_fit_plt_c, aes(age_14, alcuse, col = as.character(coa))) +
geom_path(size = 1) +
geom_text(aes(y = alcuse + 0.1, label = round(alcuse,2))) +
scale_x_continuous(breaks = c(0, 1, 2), label = c("14", "15", "16")) +
scale_color_discrete(name = "coa") +
scale_y_continuous(breaks = seq(0, 2, 0.5), limits = c(0, 2)) +
xlab("age") +
ggtitle("Model C for the effects of COA") +
theme_bw()
モデルEも可視化するが、モデルEのPEERは連続変数なので、典型的な個人を考える必要がある。取りうる全てのPEERで可視化することはもちろん可能ではあるが、線が大ぎて傾向がわからないため、典型的な個人の値を考えてモデルを可視化することが一般的。連続変数の典型的な値のとり方は下記の通り。
- 興味のある値を設定する
- パーセタイルを利用する
- 平均値に集約する
- 平均値±標準偏差×0.5(or 1)の値を利用する
ここでも書籍に従い、PEERは平均値±標準偏差×0.5の値を採用することにする。
\[ \begin{eqnarray} Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{12} PEER_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} PEER_{i} + \gamma_{10} TIME_{ij} + \gamma_{12} PEER_{i} × TIME_{ij} + (\epsilon_{ij} + \zeta_{0i} + \zeta_{1i} × TIME_{ij}) \\ &=& -0.313 + 0.571 COA_{i} + 0.695 PEER_{i} + 0.424 TIME_{ij} + -0.151 PEER_{i} × TIME_{ij} \end{eqnarray} \]
fixef.e <- fixef(model.e)
peer_prototype <- mean(alcohol1$peer) + c(-1,1) * sd(alcohol1$peer)/2
df_fit_plt_e <-
expand_grid(age_14 = 0:2, coa = 0:1, peer = peer_prototype) %>%
mutate(peer_flg = rep(c("LowPeer", "HighPeer"), times = n()/2),
type = paste0("coa=",coa," & ",peer_flg)) %>%
mutate(alcuse =
fixef.e[[1]] +
fixef.e[[2]] * coa +
fixef.e[[3]] * peer +
fixef.e[[4]] * age_14 +
fixef.e[[5]] * age_14 * peer
)
ggplot(df_fit_plt_e, aes(age_14, alcuse, col = as.character(type))) +
geom_path(size = 1) +
geom_text(aes(y = alcuse + 0.1, label = round(alcuse,2))) +
scale_x_continuous(breaks = c(0, 1, 2), label = c("14", "15", "16")) +
scale_color_manual(values = c("#ec6d71", "#a22041", "#a0d8ef", "#19448e"), name = "coa") +
scale_y_continuous(breaks = seq(0, 2, 0.5), limits = c(0, 2)) +
xlab("age") +
ggtitle("Model E for the controlled effects of COA") +
theme_bw()
親がアルコール依存症である初期値は一貫して高い位置にあり、HighPeerだと時間経過の変化が緩やかである一方、LowPeerだと時間経過の変化がHighPeerに比べるときつい。組み合わせて解釈すると、初期値において、親がアルコール依存症ではなく、周囲の友人がアルコールを摂取するケースcoa=0 & HighPeerと親がアルコール依存症で、周囲の友人がアルコールを摂取しないケースcoa=1 & LowPeerを比べると、そこまで大きな違いはないことがわかる。