マルチレベルモデリングChapter4:予測変数の中心化
UPDATE: 2022-11-16 21:55:02
はじめに
縦断データの分析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"))
)
#Models
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")
model.c <- lme(fixed = alcuse ~ coa*age_14, random = ~ age_14 | id, data = alcohol1, method = "ML")
model.d <- lme(fixed = alcuse ~ coa * age_14 + peer * age_14, random = ~ age_14 | id, data = alcohol1, method = "ML")
model.e <- lme(fixed = alcuse ~ coa + peer * age_14, random = ~ age_14 | id, data = alcohol1, method = "ML")
datatable(alcohol1 %>% mutate_if(is.numeric, round, digit = 2))予測変数の中心化
以前のノートで時間を表す変数の中心化を行うことの意味についてはまとめた。時間を表す変数から特定の定数を引くことで、切片\(\pi_{0i}\)は初期値を表すことになった。予測変数を中心化した場合、モデルの解釈はどのように変わるのであろうか。
前回のノートのモデルEを例に考える。まず、推定された\(\hat{\gamma}_{00}\)は、関連するレベル2の変数が全て0のという条件であれば、レベル1の成長パラメタ\(\pi_{0i}\)を表すが、より多くの予測変数がモデルに含まれ、すべてが0ではない場合や0がそもそも適切ではない場合などに遭遇すると、解釈が難しくなる。
\[ \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} \]
このようなケースで、予測変数を中心化しておくと解釈が簡単になる場合もある。予測変数を平均で引いて中心化すると、レベル2の当てはめられた切片は、初期値の平均的な予測値となる。平均的なとはどういうことか。
例えば、IQテストであれば、基準値は100なので、100が中心化に適した値となる。このように、中心化に用いる定数が現実的に意味を持つとき、または定数として、標本平均を利用するのであれば、それは平均的な影響となるためである。
下記は、x軸にIQスコアをとって、何らかのyとの関係を見たものとする。平均で中心化をするため、平均値の100を引くことになる。つまり、赤色が中心化したもの。中心化しても傾きに変化が出るわではないが、中心化前の平均値100が取るyの値は、中心化後の0がとるyの値と一致している。つまり、予測変数を中心化しておけば、元の変数の平均的な影響として解釈ができるようになり、それがIQスコアのように馴染みがあるものであれば、より直感的に理解しやすくなる。
set.seed(3)
X <- 80:120
Y <- 10 + 3 * X + rnorm(length(X), sd = 10)
X_c <- X - mean(X)
DF <- tibble(Y, X, X_c)
plot(Y ~ X, xlim = c(-50, 150), ylim = c(200, 400), xlab = "X or X_c", xaxt = "n")
axis(1, at = seq(-50, 150, 10))
par(new = TRUE)
plot(Y ~ X_c, xlim = c(-50, 150), ylim = c(200, 400), xlab = "", pch = 2, col = "red")
legend("bottomright", pch = 1:2, col = 1:2, legend = c("Y ~ X", "Y ~ X_c"))
LM <- lm(Y ~ X)
LM_c <- lm(Y ~ X_c)
abline(LM)
abline(LM_c, col = "red", lty=2)
segments(0, 0, 0, 309)
segments(100, 0, 100, 309)
segments(x0 = -50, x1 = 150, y0 = mean(Y), y1 = mean(Y))
予測変数を中心化したモデル
モデルFはモデルEと比べると、PEERが中心化されているC_PEERかどうかの違いしかない。初期値\(\pi_{0i}\)についてはCOAとC_PEERの効果を含み、変化率\(\pi_{1i}\)についてはC_PEERの効果を含むモデル。
\[ \begin{eqnarray} ModelF\\ Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} C\_PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{12} C\_PEER_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} COA_{i} + \gamma_{02} C\_PEER_{i} + \gamma_{10} TIME_{ij} + \gamma_{12} C\_PEER_{i} × TIME_{ij} + (\epsilon_{ij} + \zeta_{0i} + \zeta_{1i} × TIME_{ij}) \end{eqnarray} \]
モデルGは、初期値\(\pi_{0i}\)についてはC_COAとC_PEERの効果を含み、変化率\(\pi_{1i}\)についてはPEERの効果を含むモデル。
\[ \begin{eqnarray} ModelG\\ Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} C\_COA_{i} + \gamma_{02} C\_PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{12} PEER_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} C\_COA_{i} + \gamma_{02} C\_PEER_{i} + \gamma_{10} TIME_{ij} + \gamma_{12} C\_PEER_{i} × TIME_{ij} + (\epsilon_{ij} + \zeta_{0i} + \zeta_{1i} × TIME_{ij}) \end{eqnarray} \]
下記が各モデルの係数をまとめたもので、モデルE,F,Gに注目する。まず、\(\hat{\gamma}_{01}\)、\(\hat{\gamma}_{02}\)、\(\hat{\gamma}_{12}\)はモデルに使用されている変数は中心化されているかどうかの違いしかないが、回帰係数自体は同じで変化がない。標準誤差も同じであり、中心化したからといって分散成分が変化するわけでもない。各レベル2で異なっているのは、切片\(\hat{\gamma}_{00}\)、\(\hat{\gamma}_{10}\)のパラメタと標準誤差。
\[ \begin{eqnarray} Model E \\ Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \color{red}{-0.314} + 0.571 × COA_{i} + 0.695 × PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \color{red}{0.425} - 0.151 × PEER_{i} + \zeta_{1i} \\ \\ ModelF\\ Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \color{red}{0.394} + 0.571 × COA_{i} + 0.695 × C\_PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \color{red}{0.271} - 0.151 × C\_PEER_{i} + \zeta_{1i} \\ \\ ModelG\\ Y_{ij} &=& \pi_{0i} + \pi_{1i} TIME_{ij} + \epsilon_{ij} \\ \pi_{0i} &=& \color{red}{0.651} + 0.571 × C\_COA_{i} + 0.695 × C\_PEER_{i} + \zeta_{0i} \\ \pi_{1i} &=& \color{red}{0.271} - 0.151 × PEER_{i} + \zeta_{1i} \\ \\ \end{eqnarray} \]
モデルEはCOAとPEERは中心化されていないため、友人がアルコールを全く摂取せず(PEER=0)、親がアルコール依存症でない(COA=0)ときの個人の14歳時点のアルコール摂取量を表している(PEER=0 & COA=0)。
モデルFは、PEERは中心化されており、COAは中心化されていないため、PEERが平均的(C_PEER)で、親がアルコール依存症でない(COA=0)個人の14歳時点のアルコール摂取量を表している(PEER=1.08 & COA=0)。切片は、親がアルコール依存症でなく、平均的な個人を表現しており、\(\hat{\gamma}_{00}\)は、そのような子どもの初期値を表し、\(\hat{\gamma}_{10}\)は、そのような子どもの変化率を表す。
モデルGは、、PEERとCOAは中心化されているため、PEERが平均的(C_PEER)で、親がアルコール依存症でもなければ、依存症でもあるような、平均的な親(表現が難しい)をもつ個人の14歳時点のアルコール摂取量を表している。(PEER=1.08 & COA=0.451)。
モデルサマリーは下記の通りである。
| 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 | ||||||
結局のところ、今回のケースであればモデルFが望ましく、パラメタを解釈しやすいモデルになる。親がアルコール依存症でない(COA=0)個人の表す場合に、0という値は意味を持っており、PEERは中心化されているので、平均的なPEERかつ親がアルコール依存症かどうかを表現しやすいモデルであるためである。
ちなみに、モデルF、Gを実行するときは下記のスクリプトで行う。
#Model F
model.f <- lme(
fixed = alcuse ~ coa + cpeer * age_14 ,
random = ~ age_14 | id,
data = alcohol1,
method = "ML")
summary(model.f)## 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.4908388 (Intr)
## age_14 0.3730382 -0.034
## Residual 0.5807690
##
## Fixed effects: alcuse ~ coa + cpeer * age_14
## Value Std.Error DF t-value p-value
## (Intercept) 0.3938745 0.10460705 162 3.765277 0.0002
## coa 0.5711970 0.14773611 79 3.866333 0.0002
## cpeer 0.6951827 0.11240367 79 6.184697 0.0000
## age_14 0.2705847 0.06189979 162 4.371335 0.0000
## cpeer:age_14 -0.1513771 0.08538524 162 -1.772872 0.0781
## Correlation:
## (Intr) coa cpeer age_14
## coa -0.637
## cpeer 0.094 -0.146
## age_14 -0.336 0.000 0.000
## cpeer:age_14 0.000 0.000 -0.431 0.001
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.59554437 -0.40414068 -0.08351888 0.45549525 2.29975187
##
## Number of Observations: 246
## Number of Groups: 82#Model G
model.g <- lme(
fixed = alcuse ~ ccoa + cpeer * age_14,
random = ~ age_14 | id,
data = alcohol1,
method = "ML")
summary(model.g)## 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.4908391 (Intr)
## age_14 0.3730379 -0.034
## Residual 0.5807690
##
## Fixed effects: alcuse ~ ccoa + cpeer * age_14
## Value Std.Error DF t-value p-value
## (Intercept) 0.6514843 0.08060975 162 8.081954 0.0000
## ccoa 0.5711970 0.14773617 79 3.866331 0.0002
## cpeer 0.6951827 0.11240371 79 6.184695 0.0000
## age_14 0.2705847 0.06189976 162 4.371337 0.0000
## cpeer:age_14 -0.1513771 0.08538520 162 -1.772873 0.0781
## Correlation:
## (Intr) ccoa cpeer age_14
## ccoa 0.000
## cpeer 0.001 -0.146
## age_14 -0.436 0.000 0.000
## cpeer:age_14 0.000 0.000 -0.431 0.001
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.59554456 -0.40414040 -0.08351834 0.45549543 2.29975209
##
## Number of Observations: 246
## Number of Groups: 82