マルチレベルモデリングChapter5:時変な予測変数の中心化
UPDATE: 2022-11-29 21:29:01
はじめに
縦断データの分析I―変化についてのマルチレベルモデリング―を利用してマルチレベルモデリングの勉強内容をまとめたもの。下記のサポートサイトにサンプルデータなどが保存されている。
ここでの目的は、時変の予測変数の中心化についてまとめておく。
時変の予測変数の中心化
男子高校生の中退者の労働市場経験と賃金の関係を追跡したデータを測定時点の異なる場合の例とする。このデータは下記の特徴を持っている。
- 最初の時点で年齢が14歳から17歳までばらついている
- 測定の間隔が1年後であったり、2年後の場合もある
- 調査月もバラバラな状態
- 面接調査のときに2つ以上の仕事に回答できる
- 異なる時期に退学して、異なる時期に就職して、異なる時期に転職できる
lnwは自然対数変換された賃金で、\(e^{2.03}=7.6\)ドルとなる。experは時間的な変数でありlnwの観測値と関連付けられた特定の時点を表す。例えば、id=206は3回の測定回数をもっており、1.87, 2.81, 4.31年経過後を意味している。転職すると、調査があって、レコードが生成されるというイメージ。
アフリカ系を表すblack
、中退までに修了した学年を中心化したhgc、hgcは6年(小学6年)から12年(高校3年)までの値をとり、平均は8.8なので、9をひくことで、平均値に近い意味のある値の周辺に中心化(hgc.9)されている。uerateはその地域での失業率を時変な変数。
library(tidyverse)
library(broom)
library(nlme)
library(DT)
library(patchwork)
library(stargazer)
wages_pp <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/wages_pp.txt", header=T, sep=",")
datatable(wages_pp %>% filter(id %in% c(206, 332, 1028)) %>% mutate_if(is.numeric, round, 2))uerateが時変の変数で、この変数は地元地域の失業率を表す。この変数をモデルに組み込んで予測変数の中心化することの理解を深めていく。
\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{10} Time_{ij} + \gamma_{01} (HGC_{i}-9) + \gamma_{12} BLACK_{i}×TIME_{ij} + \gamma_{20}UERATE_{ij} + [\zeta_{0i} + \zeta_{1i}Time_{ij} + \epsilon_{ij}] \\ \end{eqnarray} \]
この変数をモデルに入れる際に、そのまま利用することもできるし、全平均7から引くことで中心化する方法や特定の定数を使って中心化することもできる。例えば定数7を引いた場合のモデルの\(\gamma_{00}\)の解釈は、失業率が7%の地域に住む人の平均的な対数賃金を表す。そのため、地域の失業率の1%の変化は、賃金の1.2％の減少を表す。
# 7で中心化
model.a <- lme(
lnw ~ hgc.9 + ue.7 + exper + exper:black,
random = ~exper|id,
data = wages_pp,
method = "ML")個人の平均で中心化する方法もある。下記のモデル式をみると違和感があるが、時変の予測変数を1つの変数で表すのではなく、成分を分解することで、パラメタの解釈を柔軟にできる。下記は個人内中心化法で、個人の時変変数の平均値を使って中心化する。この場合、モデルには個人の平均値と偏差を投入する。
datatable(wages_pp %>% filter(id %in% c(206, 332, 1028)) %>% mutate_if(is.numeric, round, 2) %>% select(id, lnw, exper, uerate, ue.mean, ue.person.cen)) 個人内中心化では、個人の平均は時不変、偏差は時変ということになる。その結果、このモデルでは、
- 平均的な失業率が低いほど対数賃金も低い
- 偏差は、各個人平均と比較した各時点での相対的な大きさを表す。パラメタまマイナスなので、相対的に大きくなれば対数賃金も低くなり、相対的に小さくなれば偏差はマイナスなので、対数賃金は増加する。(という解釈であっているはず。)
という2つの側面との関連を明らかにする。
# 個人平均で中心化
model.b <- lme(
lnw ~ hgc.9 + ue.mean + ue.person.cen + exper + exper:black,
random = ~exper|id,
data = wages_pp,
method = "ML")
summary(model.b)## Linear mixed-effects model fit by maximum likelihood
## Data: wages_pp
## AIC BIC logLik
## 4846.978 4914.622 -2413.489
##
## Random effects:
## Formula: ~exper | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 0.22585797 (Intr)
## exper 0.04035578 -0.332
## Residual 0.30789895
##
## Fixed effects: lnw ~ hgc.9 + ue.mean + ue.person.cen + exper + exper:black
## Value Std.Error DF t-value p-value
## (Intercept) 1.8742604 0.029537393 5511 63.45382 0
## hgc.9 0.0401695 0.006353687 885 6.32224 0
## ue.mean -0.0177091 0.003521841 885 -5.02837 0
## ue.person.cen -0.0099015 0.002098270 5511 -4.71890 0
## exper 0.0450568 0.002651046 5511 16.99584 0
## exper:black -0.0188696 0.004478993 5511 -4.21290 0
## Correlation:
## (Intr) hgc.9 ue.men u.prs. exper
## hgc.9 0.058
## ue.mean -0.926 -0.029
## ue.person.cen -0.097 -0.027 -0.013
## exper -0.187 -0.031 -0.030 0.334
## exper:black -0.112 -0.019 0.105 0.018 -0.360
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -4.34617316 -0.51627696 -0.03591668 0.44100757 7.00029530
##
## Number of Observations: 6402
## Number of Groups: 888モデルCでは時点1で中心化されているため、モデルの変数は、中退者が最初に労働市場に入ったときの初期値、この初期を基準として増加分、減少分を表していることになる。
# 時点1で中心化
model.c <- lme(
lnw ~ hgc.9 + ue1 + ue.centert1 + exper + exper:black,
random = ~exper|id,
data = wages_pp,
method = "ML")このように予測変数の中心化方法はいくつかあるが、必ずしも理解しやすくなるわけでもない。そのため、盲目的に中心化する、中心化しない、ではなく、目的に応じて、適切な中心化方法を選択する必要がある。各モデルの結果を下記にまとめておく。
| Dependent variable: | |||
| lnw | |||
| (7で中心化) | (個人内で中心化) | (時点1で中心化) | |
| Constant\(\gamma_{00}\) | 1.749*** | 1.874*** | 1.869*** |
| (0.011) | (0.030) | (0.026) | |
| (hgc-9)\(\gamma_{01}\) | 0.040*** | 0.040*** | 0.040*** |
| (0.006) | (0.006) | (0.006) | |
| uerate\(\gamma_{20}\) | -0.012*** | -0.018*** | -0.016*** |
| (0.002) | (0.004) | (0.003) | |
| (中心値からのUERATEの偏差)\(\gamma_{30}\) | -0.010*** | -0.010*** | |
| (0.002) | (0.002) | ||
| exper\(\gamma_{10}\) | 0.044*** | 0.045*** | 0.045*** |
| (0.003) | (0.003) | (0.003) | |
| exper:black\(\gamma_{12}\) | -0.018*** | -0.019*** | -0.018*** |
| (0.004) | (0.004) | (0.004) | |
| Observations | 6,402 | 6,402 | 6,402 |
| Log Likelihood | -2,415 | -2,413 | -2,412 |
| Akaike Inf. Crit. | 4,848 | 4,846 | 4,845 |
| Bayesian Inf. Crit. | 4,909 | 4,914 | 4,913 |
| Note: | p<0.1; p<0.05; p<0.01 | ||
TIMEの中心化(その2)
ここでは、TIME変数の中心化について考える。まずはサンプルデータを紹介する。このデータはうつ病患者への薬の効果を調査した研究のデータで、毎日1周間、8時、15時、22時に気分を回答したことで経時的なデータとなっている。ただ、最大で個人は21回の回答を行うことができるが、全員が必ずしもちゃんと回答しているわけではない。データの変数の意味は下記の通り。
waveは回数の連番で、1習慣を概念的に21に分割しているdayは何日目かを表すtime.of.dayは1日の各時点を0から1の範囲で表すtimeはtime.of.dayを7日まで累積したもので、ここでは尤もらしいTIME変数で0から6.67までの範囲を持つ。下記の中心化変数と比べるために表記を揃えるにおであれば、0を引いて中心化したとも考えることができる。time333は研究の中心時点である3.33という時点で中心化したものtime667は研究の最終時点である6.67という時点で中心化したもの
medication <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/medication_pp.txt", header=T, sep=",")
datatable(medication[c(1:6, 11, 16:21), c(1:8)] %>% mutate_if(is.numeric, round, 2))ここでは下記のモデルを使って、TIME変数を中心化したことによる解釈の違いを考えていく。\(C\)は中心化に利用した定数を表す。時間を中心化することで、変化率に関する\(\gamma_{10},\gamma_{11}\)は変化せず、切片に関する\(\gamma_{00},\gamma_{01}\)の値が変化する事になり、切片を評価する位置が変化することになる。つまり、初期時点、中間時点、最終時点のどのポイントで切片を評価するのかを、時間を中心化した定数によって変更できる。
\[ \begin{eqnarray} Y_{ij} &=& \pi_{0i} + \pi_{1i}(Time_{ij} - C) + \epsilon_{0i} \\ \pi_{0i} &=& \gamma_{00} + \gamma_{01} TREAT_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{11} TREAT_{i} + \zeta_{1i} \\ \\ Y_{ij} &=& \gamma_{00} + \gamma_{01} TREAT_{i} + \gamma_{10}(Time_{ij} - C) + \gamma_{11} TREAT_{i}×(Time_{ij} - C) + [\zeta_{0i} + \zeta_{1i}(Time_{ij} - C)+ \epsilon_{0i}] \\ \end{eqnarray} \]
モデルAは、時間経過とともに気分が下がる傾向(-2.41)があるが有意ではなく、初期値はランダム割付が行われているため違いがない(-3.1)。交互作用のtreat:timeは、線形変化に有意な影響があるため、変化の軌跡の傾きは異なる。つまり、treat=1の処置群では気分が時間経過とともに上がるが、treat=1
の対照群では気分が時間経過とともに下がる。
#Using time (Model A)
model.a <- lme(
pos ~ treat*time,
random= ~ time|id,
data = medication,
method = "ML")
summary(model.a)## Linear mixed-effects model fit by maximum likelihood
## Data: medication
## AIC BIC logLik
## 12696.45 12737.45 -6340.226
##
## Random effects:
## Formula: ~time | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 45.949935 (Intr)
## time 7.983475 -0.332
## Residual 35.070353
##
## Fixed effects: pos ~ treat * time
## Value Std.Error DF t-value p-value
## (Intercept) 167.46346 9.341307 1176 17.927198 0.0000
## treat -3.10930 12.352456 62 -0.251715 0.8021
## time -2.41813 1.733646 1176 -1.394822 0.1633
## treat:time 5.53681 2.281523 1176 2.426805 0.0154
## Correlation:
## (Intr) treat time
## treat -0.756
## time -0.404 0.305
## treat:time 0.307 -0.408 -0.760
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2708599 -0.4899446 -0.1174939 0.4240840 5.6180036
##
## Number of Observations: 1242
## Number of Groups: 64このモデルを可視化すると、下記のようになる。時間を中心化することで、どの時点で切片を評価するかをコントロールできる。
fixef.a <- fixef(model.a)
expand_grid(treat = c(0,1),
time = seq(from = 0, to = 7, length.out = 21)) %>%
mutate(type = paste0("treat=", treat),
pos =
fixef.a[[1]] +
fixef.a[[2]] * treat +
fixef.a[[3]] * time +
fixef.a[[4]] * time*treat
) %>% ggplot(., aes(time, pos, col = type)) +
geom_path(size = 1) +
geom_text(aes(y = pos + 1, label = round(pos,1))) +
geom_vline(xintercept = c(0, 3.33, 6.67)) +
scale_x_continuous(breaks = seq(0, 21, 1)) +
scale_y_continuous(breaks = seq(0, 200, 10)) +
xlab("time") +
theme_bw()
モデルBは切片を中間時点で評価することになるが、各線との差分を表すtreatは、この時点ではtreat=15.34で有意ではない。
#Using time - 3.33 (Model B)
model.b <- lme(
pos ~ treat * time333,
random= ~ time|id,
data = medication,
method = "ML")
summary(model.b)## Linear mixed-effects model fit by maximum likelihood
## Data: medication
## AIC BIC logLik
## 12696.45 12737.45 -6340.226
##
## Random effects:
## Formula: ~time | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 45.949935 (Intr)
## time 7.983475 -0.332
## Residual 35.070353
##
## Fixed effects: pos ~ treat * time333
## Value Std.Error DF t-value p-value
## (Intercept) 159.40304 8.778656 1176 18.158023 0.0000
## treat 15.34674 11.563284 62 1.327196 0.1893
## time333 -2.41813 1.733646 1176 -1.394822 0.1633
## treat:time333 5.53681 2.281523 1176 2.426805 0.0154
## Correlation:
## (Intr) treat tim333
## treat -0.759
## time333 0.229 -0.174
## treat:time333 -0.174 0.222 -0.760
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2708599 -0.4899446 -0.1174939 0.4240840 5.6180036
##
## Number of Observations: 1242
## Number of Groups: 64モデルCは切片を最終時点で評価することになるが、最終時点ではtreat=33.80で有意となる。これは先程、可視化した図より明らかである。
#Using time - 6.67 (Model C)
model.c <- lme(
pos ~ treat * time667,
random= ~ time|id,
data = medication,
method = "ML")
summary(model.c)## Linear mixed-effects model fit by maximum likelihood
## Data: medication
## AIC BIC logLik
## 12696.45 12737.45 -6340.226
##
## Random effects:
## Formula: ~time | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 45.949935 (Intr)
## time 7.983475 -0.332
## Residual 35.070353
##
## Fixed effects: pos ~ treat * time667
## Value Std.Error DF t-value p-value
## (Intercept) 151.34261 11.561102 1176 13.090673 0.0000
## treat 33.80278 15.182565 62 2.226421 0.0296
## time667 -2.41813 1.733646 1176 -1.394822 0.1633
## treat:time667 5.53681 2.281523 1176 2.426805 0.0154
## Correlation:
## (Intr) treat tim667
## treat -0.761
## time667 0.673 -0.513
## treat:time667 -0.512 0.670 -0.760
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2708599 -0.4899446 -0.1174939 0.4240840 5.6180036
##
## Number of Observations: 1242
## Number of Groups: 64各モデルのサマリを下記の通り、まとめておく。
| Dependent variable: | |||
| pos | |||
| (ModelA) | (ModelB_time333) | (ModelC_time667) | |
| Constant\(\gamma_{00}\) | 167.463*** | 159.403*** | 151.343*** |
| (9.341) | (8.779) | (11.561) | |
| treat\(\gamma_{01}\) | -3.109 | 15.347 | 33.803** |
| (12.352) | (11.563) | (15.183) | |
| time\(\gamma_{10}\) | -2.418 | -2.418 | -2.418 |
| (1.734) | (1.734) | (1.734) | |
| treat:time\(\gamma_{11}\) | 5.537** | 5.537** | 5.537** |
| (2.282) | (2.282) | (2.282) | |
| Observations | 1,242 | 1,242 | 1,242 |
| Log Likelihood | -6,340.226 | -6,340.226 | -6,340.226 |
| Akaike Inf. Crit. | 12,696.450 | 12,696.450 | 12,696.450 |
| Bayesian Inf. Crit. | 12,737.450 | 12,737.450 | 12,737.450 |
| Note: | p<0.1; p<0.05; p<0.01 | ||
TIMEの中心化(その2)
この他にも、時間を中心化する方法は他にもある。最小最大変換に似ている変換を施すことで、time=0のときは初期値を表し、time=6.67のときは最終時点を表すようなモデルを作ることもできる。
\[ \begin{eqnarray} Y_{ij} &=& \pi_{0i} \left( \frac{max(TIME) - TIME_{ij}}{max(TIME) - min(TIME)} \right) + \pi_{1i} \left( \frac{TIME_{ij} - min(TIME)}{max(TIME) - min(TIME)} \right) + \epsilon_{ij} \\ &=& \pi_{0i} \left( \frac{6.67 - TIME_{ij}}{6.67 - 0} \right) + \pi_{1i} \left( \frac{TIME_{ij} - 0}{6.67 - 0} \right) + \epsilon_{ij} \\ &=& \pi_{0i} \left( \frac{6.67 - TIME_{ij}}{6.67} \right) + \pi_{1i} \left( \frac{TIME_{ij}}{6.67} \right) + \epsilon_{ij} \end{eqnarray} \]
これはレベル2サブモデルとしても表現することができ、
\[ \begin{eqnarray} \pi_{0i} &=& \gamma_{00} + \gamma_{01} TREAT_{i} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \gamma_{11} TREAT_{i} + \zeta_{1i} \\ \end{eqnarray} \]
モデルを実行すると下記が得られる。対照群の初期値を167.46と推定し、初期時点での差分を- 3.11と推定する。これはモデルAの推定結果と同じであり、最終時点では、対照群の切片を151.34と推定し、差分を33.80と推定するが、これはモデルCの推定結果と同じである。
\[ \begin{eqnarray} \pi_{0i} &=& 167.46 - 3.11 × TREAT_{i} \\ \pi_{1i} &=& 151.34 + 33.80 × TREAT_{i} \\ \end{eqnarray} \]
サポートサイトにはRでの実装はのっていないが、モデルは下記の通り実行できる。
# モデル式はこれであっているのか
model.d <- lme(
pos ~ 0 + initial + final + initial:treat + final:treat,
random= ~ 0 + initial + final|id,
data = medication,
method = "ML")
summary(model.d)## Linear mixed-effects model fit by maximum likelihood
## Data: medication
## AIC BIC logLik
## 12696.45 12737.45 -6340.226
##
## Random effects:
## Formula: ~0 + initial + final | id
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## initial 45.94994 initil
## final 57.64098 0.491
## Residual 35.07035
##
## Fixed effects: pos ~ 0 + initial + final + initial:treat + final:treat
## Value Std.Error DF t-value p-value
## initial 167.46346 9.341307 1175 17.927198 0.0000
## final 151.34261 11.561102 1175 13.090673 0.0000
## initial:treat -3.10930 12.352456 1175 -0.251715 0.8013
## final:treat 33.80278 15.182565 1175 2.226421 0.0262
## Correlation:
## initil final intl:t
## final 0.404
## initial:treat -0.756 -0.306
## final:treat -0.308 -0.761 0.405
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.2708599 -0.4899446 -0.1174939 0.4240840 5.6180036
##
## Number of Observations: 1242
## Number of Groups: 64