UPDATE: 2022-11-26 14:59:20

はじめに

縦断データの分析I―変化についてのマルチレベルモデリング―を利用してマルチレベルモデリングの勉強内容をまとめたもの。下記のサポートサイトにサンプルデータなどが保存されている。

ここでの目的は、時変の予測変数を扱うマルチレベルモデルについてまとめておく。

時変の予測変数

時変の予測変数とは、時不変な変数とは異なり、名前の通り、観測タイミングによって値が変わる変数のこと。ここでは、失業が抑うつ症状に与える影響を調べた研究のデータを利用する。

今回のデータは、254人の調査対象(id)がおり、失業後すぐの0-2ヶ月後(months)、3-8ヶ月後、10−16ヶ月後に面接を行うことにより、うつ病自己評価尺度(cesd)を観測している。うつ病自己評価尺度は0から80までの範囲で値を取り、面接時点で失業しているかどうかは、失業時をunemp=1として管理され、最初の面接時点では失業(unemp=1)となる。

library(tidyverse)
library(broom)
library(nlme)
library(DT)
library(patchwork)
library(stargazer)

unemployment <- read.table("https://stats.idre.ucla.edu/stat/r/examples/alda/data/unemployment_pp.txt", header=T, sep=",")%>% 
  mutate(interview = case_when(
    months < 3 ~ 1,
    months < 8 ~ 2,
    TRUE ~ 3)
    )

datatable(unemployment %>% mutate_if(is.numeric, round, 2))

このデータの特徴は、測定回数や測定間隔も異なれば、その時点での失業状態もバラバラという点である。下記は、測定回数と失業状態の組み合わせを計算したもので、[1-0-0]は最初は失業していて、仕事を得ている状態を表す。

df_tmp <- unemployment %>%
  select(-months, -cesd) %>%
  mutate(interview = str_c("interview_", interview)) %>%
  pivot_wider(id_cols = id, names_from = interview, values_from = unemp)

df_stat <- df_tmp %>% 
  group_by(interview_1, interview_2, interview_3) %>%
  count() %>% 
  ungroup()

df_stat
## # A tibble: 11 × 4
##    interview_1 interview_2 interview_3     n
##          <int>       <int>       <int> <int>
##  1           1           0           0    55
##  2           1           0           1    19
##  3           1           0          NA     4
##  4           1           1           0    41
##  5           1           1           1    78
##  6           1           1          NA    22
##  7           1          NA           0     2
##  8           1          NA           1     4
##  9           1          NA          NA    27
## 10          NA           1           0     1
## 11          NA           1           1     1

緑は失業状態、赤は就業状態、グレーはNULLを表す。このように時変な変数は、測定タイミングに応じて、異なる値をとっている。

df_stat %>% 
  select(-n) %>% 
  mutate(rowid = row_number()) %>% 
  pivot_longer(cols = interview_1:interview_3,
               names_to = "interview",
               values_to = "value") %>% 
  mutate(value = as.character(value)) %>% 
  ggplot(., aes(interview, rowid, fill = value)) + 
  geom_tile() + 
  theme_bw()

もう少し調べると、132人は全ての面接時点で失業しており、62人は最初の面接後は就業している([1-0-0],[1-0-NA],[1-NA-0],[NA-1-0])。このように就業パターンが多様であることがわかる。

df_stat %>% 
  rowwise() %>% 
  mutate(min = min(interview_1, interview_2, interview_3, na.rm = TRUE),
         max = max(interview_1, interview_2, interview_3, na.rm = TRUE),
         all_unemp = min == max) %>% 
  group_by(all_unemp) %>% 
  summarise(total = sum(n))
## # A tibble: 2 × 2
##   all_unemp total
##   <lgl>     <int>
## 1 FALSE       122
## 2 TRUE        132

さらに理解を深めるために、各パターンに応じた個人成長プロットも可視化しておく。

unemployment %>% 
  left_join(., 
            df_tmp %>% 
              mutate_at(vars(matches("interview")), as.character) %>% 
              mutate(pattern = paste0(interview_1, "-", interview_2, "-", interview_3)) %>% 
              select(id, pattern),
            by = c("id" = "id")
  ) %>% 
  mutate(id = factor(id)) %>% 
  ggplot(., aes(months, cesd, group = id)) + 
  geom_point(size = 1) + 
  geom_line(size = 1, alpha = 0.1) + 
  scale_x_continuous(breaks = seq(0, 16, 1)) + 
  ylim(0, 80) + 
  facet_wrap( ~ pattern, nrow = 4) + 
  theme_bw() + 
  theme(legend.position = "none")

時変な変数を扱うモデル

まずは無条件成長モデルを当てはめてから時変変数を扱うマルチレベルモデルを当てはめていく。各モデルの誤差項は記載を省略したいので、ここでまとめておく。

\[ \begin{eqnarray} \begin{bmatrix} \zeta_{0i} \\ \zeta_{1i} \end{bmatrix} &\sim& N \begin{bmatrix} \sigma_{0}^{2} & \sigma_{01} \\ \sigma_{10} & \sigma_{1}^{2} \end{bmatrix},\quad \epsilon_{ij} \sim N(0, \sigma_{\epsilon}^{2}) \end{eqnarray} \]

無条件成長モデルは下記のモデル。

\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{10} TIME_{ij} + [\zeta_{0i} + \zeta_{1i} TIME_{ij} + \epsilon_{ij}] \end{eqnarray} \]

このモデルを当てはめた結果、失業初日(time=0)だとcesdは17.6点をとり、変化率は-0.42なので、時間経過とともに減少していくことがわかる。分散成分は有意であることからも、説明されていない部分が大きく、改善の余地があることがわかる。

model.a <- lme(cesd ~ months, 
               data = unemployment, 
               random = ~ months|id, 
               method = "ML")
list(
  summary(model.a),
  VarCorr(model.a)
)
## [[1]]
## Linear mixed-effects model fit by maximum likelihood
##   Data: unemployment 
##        AIC      BIC    logLik
##   5145.137 5172.217 -2566.569
## 
## Random effects:
##  Formula: ~months | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 9.3192519 (Intr)
## months      0.5958468 -0.551
## Residual    8.2975997       
## 
## Fixed effects:  cesd ~ months 
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 17.669363 0.7767173 419 22.748768       0
## months      -0.421994 0.0831023 419 -5.078006       0
##  Correlation: 
##        (Intr)
## months -0.632
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.9978128 -0.5415801 -0.1608097  0.4454676  3.5519839 
## 
## Number of Observations: 674
## Number of Groups: 254 
## 
## [[2]]
## id = pdLogChol(months) 
##             Variance   StdDev    Corr  
## (Intercept) 86.8484564 9.3192519 (Intr)
## months       0.3550334 0.5958468 -0.551
## Residual    68.8501602 8.2975997

時変変数を追加したモデル

時変変数を追加してモデルを改善しようと思っても、ここで問題になるのが、時変の予測変数unempはどのレベルに組み込めばよいかわからないということ。これを手助けをしてくれるのが合成モデル。unempは時変な予測変数なので添字\(i,j\)を取れるため、主効果として投入してみると、合成モデルは下記のとおりとなる。

\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{10} TIME_{ij}+ \gamma_{20} \color{red}{UNEMP_{ij}} + [\zeta_{0i} + \zeta_{1i} TIME_{ij} + \epsilon_{ij}] \\ \pi_{0i} &=& \gamma_{00} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \zeta_{1i} \\ \pi_{2i} &=& \gamma_{20} \end{eqnarray} \]

この合成モデルのcesdは、失業後の月数monthsと雇用状態unempに影響を受けて、残差\(\zeta_{0i}, \zeta_{1i}, \epsilon_{ij}\)が加わることで決定されると考えている。

  • \(\gamma_{00}\)は、論理的には不可能な値。失業初日(time=0)で雇用されている(unemp=0)であるため。
  • \(\gamma_{10}\)は、unempで統制した変化率の効果
  • \(\gamma_{20}\)は、monthsで統制した雇用(unemp)の効果

ただunempは時変の予測変数なので、時変のパターンが存在する。典型的なパターンとしては、

  • 常に失業中の1-1-1
  • すぐ仕事が見つかる1-0-0
  • しばらくして仕事が見つかる1-1-0
  • 再度失業するパターン1-0-1

が考えられる。

モデルBはunempの主効果を固定効果として追加した下記のモデルをあてはめる。

model.b <- lme(cesd ~ months + unemp, 
               data = unemployment, 
               random = ~ months|id, 
               method = "ML")
list(
  summary(model.b),
  VarCorr(model.b)
)
## [[1]]
## Linear mixed-effects model fit by maximum likelihood
##   Data: unemployment 
##        AIC      BIC    logLik
##   5121.603 5153.196 -2553.802
## 
## Random effects:
##  Formula: ~months | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 9.6705165 (Intr)
## months      0.6816916 -0.591
## Residual    7.8985843       
## 
## Fixed effects:  cesd ~ months + unemp 
##                 Value Std.Error  DF   t-value p-value
## (Intercept) 12.665598 1.2448444 418 10.174443  0.0000
## months      -0.201984 0.0935246 418 -2.159682  0.0314
## unemp        5.111305 0.9910523 418  5.157452  0.0000
##  Correlation: 
##        (Intr) months
## months -0.715       
## unemp  -0.780  0.459
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -1.9609562 -0.5357110 -0.1191789  0.4303293  3.6857137 
## 
## Number of Observations: 674
## Number of Groups: 254 
## 
## [[2]]
## id = pdLogChol(months) 
##             Variance   StdDev    Corr  
## (Intercept) 93.5188890 9.6705165 (Intr)
## months       0.4647034 0.6816916 -0.591
## Residual    62.3876346 7.8985843

常に失業中の1-1-1の場合、unempの効果はないので、線形で減少することになる。これは典型的なパターンを可視化しているので、実際は個人ごとに切片と変化はランダムにばらつくことになる。残りは非連続なパターンなので、unempの効果は各タイミングによって作用することになる。ここでは典型的な軌跡を描くために、5ヶ月と10ヶ月のタイミングを設定している。このモデルでは、unempの効果は定数として現れるので、ギャップが5.1程度存在する。

fixef.b <- fixef(model.b)
df_fit_plt_p1 <- 
  tibble(unemp = 1, months = seq(from = 0, to = 14, by = 1)) %>% 
  mutate(cesd = 
           fixef.b[[1]] + 
           fixef.b[[2]] * months + 
           fixef.b[[3]] * unemp
  )

p1 <- ggplot(df_fit_plt_p1, aes(months, cesd)) + 
  geom_path(size = 1) +
  geom_text(aes(y = cesd + 0.5, label = round(cesd, 0))) +
  scale_x_continuous(breaks = seq(0, 14, 1)) +
  scale_y_continuous(breaks = seq(0, 20, 1), limits = c(5, 20)) +
  xlab("months") + 
  ggtitle("Unemployment Pattern 1-1-1")


df_fit_plt_p2 <- 
  tibble(unemp  = rep(1:0, times = c(6, 10)),
         months = c(seq(0, 5, 1), seq(5, 14, 1))) %>%
  mutate(cesd = 
           fixef.b[[1]] +
           fixef.b[[2]] * months +
           fixef.b[[3]] * unemp)

p2 <- ggplot(df_fit_plt_p2, aes(months, cesd)) + 
  geom_path(size = 1) +
  geom_text(aes(y = cesd + 0.5, label = round(cesd, 0))) +
  scale_x_continuous(breaks = seq(0, 14, 1)) +
  scale_y_continuous(breaks = seq(0, 20, 1), limits = c(5, 20)) +
  xlab("months") + 
  ggtitle("Unemployment Pattern 1-0-0")


df_fit_plt_p3 <- 
  tibble(unemp  = rep(1:0, times = c(11, 5)),
         months = c(seq(0, 10, 1), seq(10, 14, 1))) %>%
  mutate(cesd = 
           fixef.b[[1]] +
           fixef.b[[2]] * months +
           fixef.b[[3]] * unemp)

p3 <- ggplot(df_fit_plt_p3, aes(months, cesd)) + 
  geom_path(size = 1) +
  geom_text(aes(y = cesd + 0.5, label = round(cesd, 0))) +
  scale_x_continuous(breaks = seq(0, 14, 1)) +
  scale_y_continuous(breaks = seq(0, 20, 1), limits = c(5, 20)) +
  xlab("months") + 
  ggtitle("Unemployment Pattern 1-1-0")


df_fit_plt_p4 <- 
  tibble(unemp  = rep(c(1, 0, 1), times = c(6, 6, 5)),
         months = c(seq(0, 5, 1), seq(5, 10, 1), seq(10, 14, 1))) %>%
  mutate(cesd = 
           fixef.b[[1]] +
           fixef.b[[2]] * months +
           fixef.b[[3]] * unemp)

p4 <- ggplot(df_fit_plt_p4, aes(months, cesd)) + 
  geom_path(size = 1) +
  geom_text(aes(y = cesd + 0.5, label = round(cesd, 0))) +
  scale_x_continuous(breaks = seq(0, 14, 1)) +
  scale_y_continuous(breaks = seq(0, 20, 1), limits = c(5, 20)) +
  xlab("months") + 
  ggtitle("Unemployment Pattern 1-0-1") 
  

(p1 | p2) / (p3 | p4) & theme_bw()

モデルAとモデルBを比べると、monthsの効果は0.42から0.20と小さくなり、unempによって調整されていることがわかる。モデル指標や分散もモデルBのほうが良くなっていることがわかる。モデルBを使って極端は対比パターンを可視化すると、よりモデルへの理解が深まる。極端なパターンとして、ここでは常に雇用されなかった個人と、4ヶ月目から雇用された個人を可視化している。個人は、雇用状況の各タイミングによって、2本の線を上下することになる。

職を得て面接可能になるのは3.5ヶ月目からであるため、書籍ではこれ以前の部分を点線で可視化しているが、ここでは点線で可視化していない、また、3.5ヶ月目ではなく、4ヶ月目として変更している。

tibble(unemp = rep(c(1,0), times = c(15, 11)), 
       months = c(seq(0, 14, 1), seq(4, 14, 1))) %>%
  mutate(type = paste0("unemp=", unemp),
         cesd = 
           fixef.b[[1]] + 
           fixef.b[[2]] * months + 
           fixef.b[[3]] * unemp
         ) %>% 
  ggplot(., aes(months, cesd, col = type)) + 
  geom_path(size = 1) +
  geom_text(aes(y = cesd + 0.5, label = round(cesd,1))) +
  scale_x_continuous(breaks = seq(0, 14, 1)) +
  scale_y_continuous(breaks = seq(0, 20, 1), limits = c(5, 20)) +
  xlab("months") + 
  ggtitle("Main effects of unemp and time, Unemployment Pattern 1-1-1 vs. x-0-0") + 
  theme_bw()

現状、このモデルBのunempのレベル2サブモデルには残差がなく固定されているが、

\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{10} TIME_{ij}+ \gamma_{20} \color{red}{UNEMP_{ij}} + [\zeta_{0i} + \zeta_{1i} TIME_{ij} + \epsilon_{ij}] \\ \pi_{0i} &=& \gamma_{00} + \zeta_{0i} \\ \pi_{1i} &=& \gamma_{10} + \zeta_{1i} \\ \pi_{2i} &=& \gamma_{20} \end{eqnarray} \] 下記のとおり、残差を追加することで、unempの効果を個人でばらつかせることも可能である。

\[ \begin{eqnarray} \pi_{2i} &=& \gamma_{20} + \zeta_{2i} \\ \begin{bmatrix} \zeta_{0i} \\ \zeta_{1i} \\ \zeta_{2i} \end{bmatrix} &\sim& N \left( \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix}, \begin{bmatrix} \sigma_{0}^{2} & \sigma_{01} & \sigma_{02} \\ \sigma_{10} & \sigma_{1}^{2} & \sigma_{12} \\ \sigma_{20} & \sigma_{21}^{2} & \sigma_{2}^{2} \end{bmatrix} \right) ,\quad \epsilon_{ij} \sim N(0, \sigma_{\epsilon}^{2}) \end{eqnarray} \]

ただ、追加できるからといって、追加するべきかどうかは理論的に妥当であるかどうかによる。つまり、unempの効果がcesdに与える効果が個人でランダムに変化するかどうか、である。また、個人で3回の観測しかないと、分散を推定するのに十分とは言えないため、境界制約を含めパラメタの推定ができなくなる可能性もある。そのため、時変な予測変数をレベル2サブモデルでランダムに変動させることは推奨されない。

交互作用を追加したモデル

ここではunempとmonthsの主効果のみを含むモデルBを改良して、モデルCではunempと線形なmonthsの交互作用を追加した下記のモデルをあてはめる。

\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{10} TIME_{ij} + \gamma_{20} UNEMP_{ij} + \gamma_{30} UNEMP_{ij} × TIME_{ij} + [\zeta_{0i} + \zeta_{1i} TIME_{ij} + \epsilon_{ij}] \end{eqnarray} \]

モデルを当てはめた結果を見ると、交互作用のmonths:unempは-0.46で有意であった。つまり、失業状態がスコアに与える影響が時間経過とともに変動し、時間経過に伴うスコアへの変化が失業状態によって異なる、ということになる。

model.c <- lme(cesd ~ months*unemp, 
               data = unemployment, 
               random = ~ months|id, 
               method = "ML")
list(
  summary(model.c),
  VarCorr(model.c)
)
## [[1]]
## Linear mixed-effects model fit by maximum likelihood
##   Data: unemployment 
##        AIC      BIC    logLik
##   5119.047 5155.153 -2551.523
## 
## Random effects:
##  Formula: ~months | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev    Corr  
## (Intercept) 9.6805552 (Intr)
## months      0.6717203 -0.596
## Residual    7.8759887       
## 
## Fixed effects:  cesd ~ months * unemp 
##                  Value Std.Error  DF   t-value p-value
## (Intercept)   9.616744 1.8949397 417  5.074961  0.0000
## months        0.162036 0.1942395 417  0.834207  0.4046
## unemp         8.529059 1.8834713 417  4.528372  0.0000
## months:unemp -0.465222 0.2178620 417 -2.135398  0.0333
##  Correlation: 
##              (Intr) months unemp 
## months       -0.888              
## unemp        -0.911  0.863       
## months:unemp  0.755 -0.878 -0.852
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.0039251 -0.5237204 -0.1548188  0.4308929  3.6960920 
## 
## Number of Observations: 674
## Number of Groups: 254 
## 
## [[2]]
## id = pdLogChol(months) 
##             Variance   StdDev    Corr  
## (Intercept) 93.7131494 9.6805552 (Intr)
## months       0.4512082 0.6717203 -0.596
## Residual    62.0311982 7.8759887

言葉では理解しにくいので可視化してみると、予期していないパターンの直線となっていることがわかる。つまり、雇用されている場合(unemp=0)においては、cesdが上昇する可能性があることをモデルは示している。monthsの効果は0.16で有意ではなく、この場合、0である可能性もありうる。この結果より、再雇用者の直線の傾きを0にし、再雇用者は時間経過とともに変化しないと制約をかけることが望ましいとも考えられる。

書籍では、このように再雇用者のcesdの上昇は望ましくないと考えているが、実際、雇用されたらされたらで、cesdが上昇することは有り得そうな話である。仕事をしていることで、様々な要因(特に人間関係や職務評価など)によって、うつ病に追い込まれるケースは多そうなので。

fixef.c <- fixef(model.c)
tibble(unemp = rep(c(1,0), times = c(15, 11)), 
       months = c(seq(0, 14, 1), seq(4, 14, 1))) %>%
  mutate(type = paste0("unemp=", unemp),
         cesd = 
           fixef.c[[1]] + 
           fixef.c[[2]] * months + 
           fixef.c[[3]] * unemp + 
           fixef.c[[4]] * unemp*months 
  ) %>%   ggplot(., aes(months, cesd, col = type)) + 
  geom_path(size = 1) +
  geom_text(aes(y = cesd + 0.5, label = round(cesd,1))) +
  scale_x_continuous(breaks = seq(0, 14, 1)) +
  scale_y_continuous(breaks = seq(0, 20, 1), limits = c(5, 20)) +
  xlab("months") + 
  ggtitle("Main effects of unemp and time, Unemployment Pattern 1-1-1 vs. x-0-0") + 
  theme_bw()

モデルに雇用されている場合(unemp=0)に、傾きが水平になる制約を入れたいのであれば、\(\gamma_{10} TIME_{ij}\)をモデルから抜けば良いように見えるが、つまり、雇用されている場合(unemp=0)には、unempの効果しか作用しなくなるが、これは構造的な部分と確率的な部分が一致しておらず、モデルとして成立していないので、当てはめることができない。

\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{20} UNEMP_{ij} + \gamma_{30} UNEMP_{ij} × TIME_{ij} + [\zeta_{0i} + \zeta_{1i} TIME_{ij} + \epsilon_{ij}] \end{eqnarray} \]

\(TIME_{ij}\)のランダム効果\(\zeta_{1i}\)が含まれるが主効果\(\gamma_{10}\)は含まない、交互作用\(\gamma_{30}\)を含む一方で対応するランダム効果\(\zeta_{3i}\)は含まないため、モデルとして成立していない。そのため、モデルを当てはまえるためには、下記のように効果とランダム要素を整える必要がある。

\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{20} UNEMP_{ij} + \gamma_{30} UNEMP_{ij} × TIME_{ij} + [\zeta_{0i} + \zeta_{3i} UNEMP_{ij} × TIME_{ij} + \epsilon_{ij}] \end{eqnarray} \]

ただ、このモデルは雇用されている個人の傾きを成約するために作ったモデルでした。\(\zeta_{0i}\)を含めることで、雇用された個人の軌跡がランダムにばらつくことを認めながら、\(\zeta_{2i}\)がないので、雇用されていない個人の切片の増分\(\gamma_{20}\)がランダムにばらつくことを許していない。つまり、再雇用者の水平の水準がランダムに変動することを認めながら、この水準の増加がランダムにばらつくことは許さない、という非常に厳しい制約を課していることになる。そのため、このモデルの当てはまりは非常に悪くなってしまう。

制約つき交互作用を追加したモデル

これを検証するためにモデルDを考える。モデルDはunempが固定効果とランダム効果を持ち、各固定効果にランダム効果を持つことを許容している。

\[ \begin{eqnarray} Y_{ij} &=& \gamma_{00} + \gamma_{20} UNEMP_{ij} + \gamma_{30} UNEMP_{ij} × TIME_{ij} + [\zeta_{0i} + \zeta_{2i} UNEMP_{ij}+ \zeta_{3i} UNEMP_{ij} × TIME_{ij} + \epsilon_{ij}] \end{eqnarray} \]

モデルDは、サポートサイトのコードでは収束しないため、少し調整する必要がある。ただ、この設定では不要な調整も含まれている可能性がある。

# https://groups.google.com/g/davis-rug/c/sbcjeMalCgc?pli=1
# 調整しないと収束しない 
# Error in lme.formula(cesd ~ unemp + unemp:months, unemployment, random = ~unemp +  : 
#   nlminb problem, convergence error code = 1
#   message = iteration limit reached without convergence (10)

control.list <-
  lmeControl(
    maxIter = 500,
    msMaxIter = 500,
    msMaxEval = 500,
    tolerance = 0.1,
    msTol = 0.1,
    sing.tol = 1e-20
  )

model.d <-
  lme(
    cesd ~ unemp + unemp:months,
    random =  ~ unemp + unemp:months | id,
    data = unemployment,
    control = control.list
  )

list(
  summary(model.d),
  VarCorr(model.d)
)
## [[1]]
## Linear mixed-effects model fit by REML
##   Data: unemployment 
##        AIC      BIC    logLik
##   5115.578 5160.666 -2547.789
## 
## Random effects:
##  Formula: ~unemp + unemp:months | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##              StdDev    Corr         
## (Intercept)  6.7644554 (Intr) unemp 
## unemp        6.7554117  0.137       
## unemp:months 0.8746322  0.117 -0.968
## Residual     7.6889727              
## 
## Fixed effects:  cesd ~ unemp + unemp:months 
##                  Value Std.Error  DF   t-value p-value
## (Intercept)  11.197231 0.7924315 418 14.130219  0.0000
## unemp         6.924770 0.9326717 418  7.424659  0.0000
## unemp:months -0.303075 0.1124915 418 -2.694201  0.0073
##  Correlation: 
##              (Intr) unemp 
## unemp        -0.564       
## unemp:months -0.073 -0.444
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.1282215 -0.4899561 -0.1572776  0.4260313  3.7742403 
## 
## Number of Observations: 674
## Number of Groups: 254 
## 
## [[2]]
## id = pdLogChol(unemp + unemp:months) 
##              Variance   StdDev    Corr         
## (Intercept)  45.7578567 6.7644554 (Intr) unemp 
## unemp        45.6355876 6.7554117  0.137       
## unemp:months  0.7649815 0.8746322  0.117 -0.968
## Residual     59.1203010 7.6889727

職を失った直後(month=0)では、失業者のcesdは18.12(=11.20 + 6.92)で、月ごとに-0.30ポイントづつ低下していき、再び雇用されると 6.92ポイント下がることになる。

fixef.d <- fixef(model.d)
tibble(unemp = rep(c(1,0), times = c(15, 11)), 
       months = c(seq(0, 14, 1), seq(4, 14, 1))) %>%
  mutate(type = paste0("unemp=", unemp),
         cesd = 
           fixef.d[[1]] + 
           fixef.d[[2]] * unemp + 
           fixef.d[[3]] * unemp*months 
  ) %>% 
  ggplot(., aes(months, cesd, col = type)) + 
  geom_path(size = 1) +
  geom_text(aes(y = cesd + 0.5, label = round(cesd,1))) +
  scale_x_continuous(breaks = seq(0, 14, 1)) +
  scale_y_continuous(breaks = seq(0, 20, 1), limits = c(5, 20)) +
  xlab("months") + 
  ggtitle("Constraining the effects time among the re-employed,\nUnemployment Pattern 1-1-1 vs. x-0-0") + 
  theme_bw()

このモデルDは、モデルCと比べて、AIC、BICもあまり変わらないため、研究者らの仮説をより良く表現されている。ここまでに当てはめたモデルのサマリを下記にまとめておく。

Dependent variable:
cesd
ModelA ModelB ModelC ModelD
Constant\(\gamma_{00}\) 17.669*** 12.666*** 9.617*** 11.19***
(0.777) (1.245) (1.895) (0.792)
months\(\gamma_{10}\) -0.422*** -0.202** 0.162
(0.083) (0.094) (0.194)
unemp\(\gamma_{20}\) 5.111*** 8.529*** 6.92***
(0.991) (1.883) (0.932)
months:unemp\(\gamma_{30}\) -0.465** -0.30**
(0.218) (0.112)
個人内\(\sigma^{2}_{\epsilon}\) 68.85 62.39 62.03 59.12
切片\(\sigma^{2}_{0}\) 86.85 93.52 93.71 45.75
変化率\(\sigma^{2}_{1}\) 0.36 0.46 0.45
UNEMP\(\sigma^{2}_{2}\) 45.63
UNEMP:TIME\(\sigma^{2}_{3}\) 0.76
Observations 674 674 674 674
Log Likelihood -2,566 -2,553 -2,551 -2,547
Akaike Inf. Crit. 5,145 5,121 5,119 5,115
Bayesian Inf. Crit. 5,172 5,153 5,155 5,160
Note: p<0.1; p<0.05; p<0.01