1 Introduction

The data set is from a psychiatric study focusing on the clinical response over time in 66 depressed patients. Depression is often classified into two types: non-endogenous depression, which is associated with some tragic life event such as the death of a close friend or family member, or endogenous depression, which is not a result of any specific event. Antidepressant medications are held to be more effective for endogenous depression. In this sample, 29 patients were diagnosed as non-endogenous, and 37 patients were endogenous.

Following a placebo period of 1 week, patients received 225 mg/day doses of imipramine for 4 weeks. Patients were rated with the Hamilton depression (HD) rating scale twice during the baseline placebo week (at the start and end of this week), as well as at the end of each of the 4 treatment weeks of the study.

The number of patients with all measures at each of the weeks fluctuated. Only 46 had complete data at all time points.

For the 17-item version, score on the Hamilton Scale goes from 0 to 54. One suggestion is that scores between 0 and 6 indicate a normal person with regard to depression, scores between 7 and 17 indicate mild depression, scores between 18 and 24 indicate moderate depression, and scores over 24 indicate severe depression.

2 Data

# data input
dta <- read.csv("C:/Users/openw/Documents/riesby_data.csv")
dta <- dta[, c("ID", "Endog", "week", "HamD")]
names(dta) <- c("ID", "Endogenous", "Week", "Score")

head(dta)

   ID Endogenous Week Score
1 101          0    0    26
2 101          0    1    22
3 101          0    2    18
4 101          0    3     7
5 101          0    4     4
6 101          0    5     3

str(dta)

'data.frame':   396 obs. of  4 variables:
 $ ID        : int  101 101 101 101 101 101 103 103 103 103 ...
 $ Endogenous: int  0 0 0 0 0 0 0 0 0 0 ...
 $ Week      : int  0 1 2 3 4 5 0 1 2 3 ...
 $ Score     : int  26 22 18 7 4 3 33 24 15 24 ...

# verify the number of endogenous cases
aggregate(Endogenous ~ ID, data=dta, mean)[,2] |> sum()

[1] 37

# make a copy of the original data in the long format
dtaL <- dta %>%                                        
 dplyr::filter(!is.na(Score)) %>%          # remove NA or missing scores
 dplyr::mutate(ID = as.factor(ID),
               Endogenous = as.factor(Endogenous) %>% 
               forcats::fct_recode("N" = "0",
                                   "Y" = "1")) %>%
 dplyr::arrange(ID, Week)                     # sort observations

# long to wide format
dtaw <- dtaL %>%     
 tidyr::pivot_wider(names_from = Week,
                    names_prefix = "Score_",
                    values_from = Score)

3 Summary statistics

dtaw %>% 
  furniture::table1("Endogenous Depression" = Endogenous, 
                    align="c",
                    output = "markdown")

	Mean/Count (SD/%)
	n = 66
Endogenous Depression
N	29 (43.9%)
Y	37 (56.1%)

# number of patients by week
dtaL %>% 
  dplyr::group_by(Week) %>% 
  dplyr::tally() %>%
  knitr::kable()

Week	n
0	61
1	63
2	65
3	65
4	63
5	58

3.1 Missing data

# which indices are true and return array indices
which(is.na(dtaw), arr.ind=TRUE)

      row col
 [1,]   9   3
 [2,]  10   3
 [3,]  12   3
 [4,]  40   3
 [5,]  65   3
 [6,]  24   4
 [7,]  39   4
 [8,]  66   4
 [9,]  35   5
[10,]  18   6
[11,]   6   7
[12,]  23   7
[13,]  66   7
[14,]   5   8
[15,]   8   8
[16,]  22   8
[17,]  26   8
[18,]  30   8
[19,]  43   8
[20,]  50   8
[21,]  61   8

-有遺漏值的儲存格的位置

# compare original data set with missing removed 
dtaw %>%
     dplyr::setdiff(., na.omit(.))

# A tibble: 20 × 8
   ID    Endogenous Score_0 Score_1 Score_2 Score_3 Score_4 Score_5
   <fct> <fct>        <int>   <int>   <int>   <int>   <int>   <int>
 1 106   Y               21      25      23      18      20      NA
 2 107   Y               21      21      16      19      NA       6
 3 113   N               21      23      19      23      23      NA
 4 114   N               NA      17      11      13       7       7
 5 115   Y               NA      16      16      16      16      11
 6 118   Y               NA      26      18      18      14      11
 7 304   Y               21      27      29      NA      12      24
 8 310   Y               24      19      11       7       6      NA
 9 311   Y               20      16      21      17      NA      15
10 312   Y               17      NA      18      17      17       6
11 315   Y               27      21      17      13       5      NA
12 322   Y               28      21      25      32      34      NA
13 334   N               31      25      NA       7       8      11
14 339   Y               27      NA      14      12      11      12
15 344   Y               NA      21      12      13      13      18
16 347   Y               18      15      14      10       8      NA
17 354   Y               28      27      27      26      23      NA
18 604   N               27      27      13       5       7      NA
19 609   Y               NA      25      22      14      15       2
20 610   Y               34      NA      33      23      NA      11

-有遺漏值的項目外，還有受試者編號跟其他分數的情形。

dtaw %>% 
  naniar::gg_miss_upset(sets = c("Score_5_NA",
                                 "Score_4_NA",
                                 "Score_3_NA",
                                 "Score_2_NA",
                                 "Score_1_NA",
                                 "Score_0_NA"),
                        keep.order = TRUE)

-每週都有病人沒有完成HD量表，第一週和第四週皆有3名病人沒有完成，其中第五週最多，有8名病人。

dtaw %>%                
  dplyr::group_by(Endogenous) %>% 
  furniture::table1("Baseline" = Score_0,
                    "Week 1" = Score_1,
                    "Week 2" = Score_2,
                    "Week 3" = Score_3,
                    "Week 4" = Score_4,
                    "Week 5" = Score_5,
                    total = TRUE,
                    test = TRUE,
                    na.rm = FALSE,   # default: COMPLETE CASES ONLY
                    digits = 2,
                    align = "c",
                    output = "markdown",
                    caption = "Weekly Hamilton Depression Scores by Diagnosis for all participants")

Weekly Hamilton Depression Scores by Diagnosis for all participants
	Total	N	Y	P-Value
	n = 66	n = 29	n = 37
Baseline				0.301
	23.44 (4.53)	22.79 (4.12)	24.00 (4.85)
Week 1				0.033
	21.84 (4.70)	20.48 (3.83)	23.00 (5.10)
Week 2				0.095
	18.31 (5.49)	17.00 (4.35)	19.30 (6.08)
Week 3				0.23
	16.42 (6.42)	15.34 (6.17)	17.28 (6.56)
Week 4				0.298
	13.62 (6.97)	12.62 (6.72)	14.47 (7.17)
Week 5				0.48
	11.95 (7.22)	11.22 (6.34)	12.58 (7.96)

-整體而言，隨著時間增加，HD量表的分數下降。其中不是Endogenous depression的病人分數比Endogenous depression的病人分數低，在接受治療後的第二週分數有達到顯著差異。

3.1.1 Problem

Summarize Hamilton depression scores across time by depression type for patients with complete data for all 6 weeks

complete.cases(dtaw)

 [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE
[13]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE
[25]  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE
[37]  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
[49]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[61] FALSE  TRUE  TRUE  TRUE FALSE FALSE

rm.dtaw <- dtaw[complete.cases(dtaw), ]

library(tidyr)

#pivot the data frame into a long format
rm.dtaw %>% 
  pivot_longer(cols=c('Score_0', 'Score_1', 'Score_2', 'Score_3', 'Score_4', 'Score_5'),
                         names_to='Week',
                         values_to='Score')

# A tibble: 276 × 4
   ID    Endogenous Week    Score
   <fct> <fct>      <chr>   <int>
 1 101   N          Score_0    26
 2 101   N          Score_1    22
 3 101   N          Score_2    18
 4 101   N          Score_3     7
 5 101   N          Score_4     4
 6 101   N          Score_5     3
 7 103   N          Score_0    33
 8 103   N          Score_1    24
 9 103   N          Score_2    15
10 103   N          Score_3    24
# … with 266 more rows

#Save as data
rm.dtaw_L <- rm.dtaw %>% 
  pivot_longer(cols=c('Score_0', 'Score_1', 'Score_2', 'Score_3', 'Score_4', 'Score_5'),
                         names_to='Week',
                         values_to='Score')

data_sum_all <- rm.dtaw_L %>% 
  dplyr::group_by(Endogenous,Week) %>%             # specify the groups
  dplyr::summarise(Score_n    = n(),                # count of valid scores
                   Score_mean = mean(Score),        # mean score
                   Score_sd   = sd(Score),          # standard deviation of scores
                   Score_sem  = Score_sd / sqrt(Score_n)) # se for the mean

data_sum_all |> knitr::kable()

Endogenous	Week	Score_n	Score_mean	Score_sd	Score_sem
N	Score_0	25	22.3600	3.90384	0.780769
N	Score_1	25	20.0800	3.68465	0.736931
N	Score_2	25	17.3200	4.33705	0.867410
N	Score_3	25	15.8800	5.84038	1.168075
N	Score_4	25	12.8400	6.68132	1.336263
N	Score_5	25	11.4000	6.53835	1.307670
Y	Score_0	21	24.0952	4.86729	1.062129
Y	Score_1	21	23.9048	5.46722	1.193045
Y	Score_2	21	18.9524	6.01229	1.311989
Y	Score_3	21	17.4762	6.86017	1.497011
Y	Score_4	21	14.1905	6.98297	1.523810
Y	Score_5	21	13.0476	8.72626	1.904226

-有完成每次量表的病人總共有46人，其中25人不是Endogenous depression,21人是Endogenous depression。

3.2 Covariances and correlations

dtaw %>% 
  dplyr::select(starts_with("Score_")) %>%  
  cov(use = "pairwise.complete.obs")  %>%  
  round(3)

        Score_0 Score_1 Score_2 Score_3 Score_4 Score_5
Score_0  20.551  10.115  10.139  10.086   7.191   6.278
Score_1  10.115  22.071  12.277  12.550  10.264   7.720
Score_2  10.139  12.277  30.091  25.126  24.626  18.384
Score_3  10.086  12.550  25.126  41.153  37.339  23.992
Score_4   7.191  10.264  24.626  37.339  48.594  30.513
Score_5   6.278   7.720  18.384  23.992  30.513  52.120

3.2.1 Problem

Compute the covariance matrix of the scores over the weeks using complete cases only

rm.dtaw %>% 
  dplyr::select(starts_with("Score_")) %>%  
  cov(use = "pairwise.complete.obs")  %>%  
  round(3)

        Score_0 Score_1 Score_2 Score_3 Score_4 Score_5
Score_0  19.421  10.716   9.523  12.350   9.062   7.376
Score_1  10.716  24.236  12.545  15.930  11.592   8.471
Score_2   9.523  12.545  26.773  23.848  23.858  20.657
Score_3  12.350  15.930  23.848  39.755  33.316  29.728
Score_4   9.062  11.592  23.858  33.316  45.943  37.107
Score_5   7.376   8.471  20.657  29.728  37.107  57.332

3.2.2 Problem

Compute the correlation matrix of the scores over the weeks as long as pair-wise cases are available

dta_score<-rm.dtaw[,c('Score_0','Score_1','Score_2','Score_3','Score_4','Score_5')]
round(cor(dta_score),3)

        Score_0 Score_1 Score_2 Score_3 Score_4 Score_5
Score_0   1.000   0.494   0.418   0.444   0.303   0.221
Score_1   0.494   1.000   0.492   0.513   0.347   0.227
Score_2   0.418   0.492   1.000   0.731   0.680   0.527
Score_3   0.444   0.513   0.731   1.000   0.780   0.623
Score_4   0.303   0.347   0.680   0.780   1.000   0.723
Score_5   0.221   0.227   0.527   0.623   0.723   1.000

dtaw %>% 
  dplyr::select(starts_with("Score_")) %>% # just the outcome(s)
  cor(use = "pairwise.complete.obs") %>%   # correlation matrix
  corrplot::corrplot.mixed(upper = "ellipse")

dtaw %>% 
  dplyr::filter(Endogenous == "Y") %>%    # for the Endogenous group
  dplyr::select(starts_with("Score_")) %>%
  cor(use = "pairwise.complete.obs") %>%   
  corrplot::corrplot.mixed(upper = "ellipse")

3.2.3 Problem

Draw the correlation plot for scores over the weeks for the non-endogenous group

dtaw %>% 
  dplyr::filter(Endogenous == "N") %>%    # for the Endogenous group
  dplyr::select(starts_with("Score_")) %>%
  cor(use = "pairwise.complete.obs") %>%   
  corrplot::corrplot.mixed(upper = "ellipse")

4 Visualization

# Use theme_set() to completely override the current theme.
old <- theme_set(theme_minimal())

ggplot(dtaL, aes(x = Week, y = Score)) +
  geom_line(aes(group = ID), alpha=.4, size=rel(.4)) +
  facet_grid( ~ Endogenous) +
  scale_y_continuous(limits=c(0, 42), breaks=seq(0, 42, by=6))+
  labs(x="Time (weeks since baseline)",
       y="Hamilton Depression Score",
       subtitle="Endogenous depression")

-在不同組別的每位病人HD量表分數表現，整體來說，隨著時間增加，HD量表分數會逐漸下降，其中Endogenous depression病人的分數比不是Endogenous depression病人的分數來得低。

ggplot(dtaL, aes(x = Week, y = Score, color=Endogenous)) +
  stat_summary(fun.data="mean_se", position=position_dodge(width=0.2)) +
  stat_summary(aes(group=Endogenous, linetype=Endogenous), 
               fun = mean,
               geom = "line",
               position=position_dodge(width=.2)) +
  scale_linetype_manual(values = c("solid", "longdash")) +
  scale_color_manual(values=c('black','gray50'))+
  labs(x="Time (weeks since baseline)",
       y="Mean Hamilton Depression Score")+
  theme(legend.position = "bottom",
        legend.key.width = unit(2, "cm"))

-在不同組別的HD量表分數表現，整體來說，隨著時間增加，HD量表分數會逐漸下降，其中Endogenous depression病人的分數比不是Endogenous depression病人的分數來得低。

ggplot(dtaL, aes(x = Week, y = Score)) +
  geom_line(aes(group = ID), col='gray', alpha=.5, size=rel(.5)) +
  facet_grid( ~ Endogenous) +
  geom_smooth(method="loess", 
              se=FALSE, 
              formula=y~x,
              color="orange") +  
  geom_smooth(method = "lm", 
              formula=y~x,
              se = FALSE, 
              color="skyblue")+
  scale_y_continuous(limits=c(0, 42), breaks=seq(0, 42, by=6))+
  labs(x="Time (weeks since baseline)",
       y="Hamilton Depression Score")

-以區域性迴歸(橘線)及線性迴歸(藍線)做出的迴歸模型

ggplot(dtaL, aes(x = Week,
             y = Score,
             group = Endogenous,
             linetype = Endogenous)) +
  geom_smooth(method = "loess",
              formula=y~x,
              color = "black",
              alpha = .2) +
  scale_linetype_manual(values = c("solid", "longdash")) +
  labs(x = "Time (weeks since baseline)",
       y = "Hamilton Depression Scale",
       linetype = "Endogenous Depression")+
  theme(legend.position = c(1, 1),
        legend.justification = c(1.1, 1.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(1.5, "cm"))

-在區域迴歸分析中，不同組的分數信賴區間範圍有相互重疊，推測即使不是Endogenous depression病人，HD量表的分數也可能會跟是Endogenous depression病人分數相同。

5 Models

5.1 Random intercepts only

# 1 is included whether you type it in or not
m0 <- lme4::lmer(Score ~ 1 + Week + (1 | ID), data=dtaL)

texreg::knitreg(m0, 
                single.row = TRUE,
                stars = numeric(0),
                caption = "Random Intercepts Model",
                caption.above = TRUE,
                custom.note = "Model fit w/ REML")

Random Intercepts Model
	Model 1
(Intercept)	23.55 (0.64)
Week	-2.38 (0.14)
AIC	2294.73
BIC	2310.43
Log Likelihood	-1143.36
Num. obs.	375
Num. groups: ID	66
Var: ID (Intercept)	16.45
Var: Residual	19.10
Model fit w/ REML

-在治療開始前，全體平均HD量表分數為23.55。時間每增加1單位，每位病人的HD量表分數會下降2.38。ICC為0.462(=16.45/(16.45+19.1))，表示在相同時間的兩個分數之間的相似性有46.2%

5.1.1 Fitted values

dtaL %>% 
  dplyr::mutate(pred_fixed = predict(m0, re.form = NA)) %>% # fixed effects only
  dplyr::mutate(pred_wrand = predict(m0)) %>%               # fixed and random effects together
  ggplot(aes(x = Week, y = Score, group = ID)) +
  geom_line(aes(y        = pred_wrand,
                color    = "BLUP",
                size     = "BLUP",
                linetype = "BLUP")) +
  geom_line(aes(y        = pred_fixed,
                color    = "Marginal",
                size     = "Marginal",
                linetype = "Marginal")) +
  scale_color_manual(name   = "Type of Prediction",
                     values = c("BLUP"     = "gray50",
                                "Marginal" = "blue"))  +
  scale_size_manual(name    = "Type of Prediction",
                    values = c("BLUP"      = .8,
                               "Marginal" = .5))  +
  scale_linetype_manual(name   = "Type of Prediction",
                        values = c("BLUP"     = "dotted",
                                   "Marginal" = "solid")) +
  labs(x = "Weeks since baseline",
       y = "Hamilton Depression Scores")+
  theme(legend.position = c(0, 0),
        legend.justification = c(-0.1, -0.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(1.5, "cm"))

obs <- dtaL %>% 
  dplyr::group_by(Week) %>% 
  dplyr::summarise(observed = mean(Score, na.rm = TRUE))

effects::Effect(focal.predictors = "Week",
                mod = m0,
                xlevels = list(Week = 0:5)) %>% 
  data.frame() %>% 
  dplyr::rename(estimated = fit) %>% 
  dplyr::left_join(obs, by = "Week") %>% 
  dplyr::select(Week, observed, estimated) %>% 
  dplyr::mutate(diff = observed - estimated) %>% 
  pander::pander(caption = "Observed and Estimated Means")

Observed and Estimated Means
Week	observed	estimated	diff
0	23.44	23.55	-0.109
1	21.84	21.18	0.6652
2	18.31	18.8	-0.4928
3	16.42	16.42	-0.009553
4	13.62	14.05	-0.4303
5	11.95	11.67	0.2745

-在不同測量時間，估計平均HD量表分數，與測量到的HD量表分數，整體來說隨時間增加，分數會逐漸下降。

performance::icc(m0)

# Intraclass Correlation Coefficient

    Adjusted ICC: 0.463
  Unadjusted ICC: 0.319

-調整後的ICC為0.463，表示在相同時間不同病人的HD量表分數的相似性有46.3%。

# fit a linear model ignoring patient clusters
m00 <- lm(Score ~ Week, data = dtaL)

texreg::knitreg(list(m00, m0),
                custom.model.names = c("Linear model", "Multilevel"),
                single.row = TRUE,
                stars = numeric(0),
                caption = "Random Intercepts Models",
                caption.above = TRUE,
                custom.note = "")

Random Intercepts Models
	Linear model	Multilevel
(Intercept)	23.60 (0.55)	23.55 (0.64)
Week	-2.41 (0.18)	-2.38 (0.14)
R²	0.32
Adj. R²	0.32
Num. obs.	375	375
AIC		2294.73
BIC		2310.43
Log Likelihood		-1143.36
Num. groups: ID		66
Var: ID (Intercept)		16.45
Var: Residual		19.10

-以線性迴歸及多層次統計分析的比較，雖然兩種分析方式斜率及截距相似，但是線性迴歸不考慮個體差異，即多層次統計分析的隨機截距。

logLik(m00)

'log Lik.' -1199.86 (df=3)

sigma(m00)^2

[1] 35.3997

VarCorr(m0) %>% print(., comp=c("Variance","Std.Dev."), digits=4)

 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 16.45    4.055   
 Residual             19.10    4.370

lme4::VarCorr(m0) %>% 
  print(comp = c("Variance", "Std.Dev"))

 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 16.45    4.055   
 Residual             19.10    4.370

5.2 Random slopes and intercepts model

m1 <- lme4::lmer(Score ~ 1 + Week + (1 + Week | ID), data = dtaL)

texreg::knitreg(list(m0, m1),
                single.row = TRUE,
                stars = numeric(0),
                custom.model.names = c("Random Intercepts",
                                       " And Random Slopes"),
                caption = "MLM: Models fit w/REML",
                caption.above = TRUE,
                custom.note = "")

MLM: Models fit w/REML
	Random Intercepts	And Random Slopes
(Intercept)	23.55 (0.64)	23.58 (0.55)
Week	-2.38 (0.14)	-2.38 (0.21)
AIC	2294.73	2231.92
BIC	2310.43	2255.48
Log Likelihood	-1143.36	-1109.96
Num. obs.	375	375
Num. groups: ID	66	66
Var: ID (Intercept)	16.45	12.94
Var: Residual	19.10	12.21
Var: ID Week		2.13
Cov: ID (Intercept) Week		-1.48

-模型中包含隨機截距的情形下，在治療開始前，全體平均HD量表分數為23.55。時間每增加1單位，每位病人的HD量表分數會下降2.38。ICC為0.462(=16.45/(16.45+19.1))，表示在相同時間的兩個分數之間的相似性有46.2%。

-模型中有隨機截距跟隨機變項的情形下，在治療開始前，全體平均HD量表分數為23.58。時間每增加1單位，每位病人的HD量表分數會下降2.38。ICC為0.474(=12.94/(12.94+12.21+2.13))，表示在相同時間的兩個分數之間的相似性有47.4%。AIC為2231.92，比起m0，m1的模型適配度較佳。

anova(m0, m1, 
      model.names = c("Random Intercepts", "Random Intercepts/Slopes"),
      refit = TRUE) %>% 
  pander::pander(caption = "LRT: Significance of Random Slopes")

refitting model(s) with ML (instead of REML)

LRT: Significance of Random Slopes (continued below)
	npar	AIC	BIC	logLik	deviance	Chisq
Random Intercepts	4	2293	2309	-1143	2285	NA
Random Intercepts/Slopes	6	2231	2255	-1110	2219	66.15

	Df	Pr(>Chisq)
Random Intercepts	NA	NA
Random Intercepts/Slopes	2	4.319e-15

dtaL %>% 
  dplyr::mutate(pred_fixed = predict(m1, re.form = NA)) %>% # fixed effects only
  dplyr::mutate(pred_wrand = predict(m1)) %>%               # fixed and random effects together
  ggplot(aes(x = Week,
             y = Score,
             group = ID)) +
  geom_line(aes(y        = pred_wrand,
                color    = "BLUP",
                size     = "BLUP",
                linetype = "BLUP")) +
  geom_line(aes(y        = pred_fixed,
                color    = "Marginal",
                size     = "Marginal",
                linetype = "Marginal")) +
  scale_color_manual(name   = "Type of Prediction",
                     values = c("BLUP"     = "gray50",
                                "Marginal" = "blue"))  +
  scale_size_manual(name    = "Type of Prediction",
                    values = c("BLUP"      = .7,
                                "Marginal" = .5))  +
  scale_linetype_manual(name   = "Type of Prediction",
                        values = c("BLUP"     = "longdash",
                                   "Marginal" = "solid"))+
  labs(x = "Weeks Since Baseline",
       y = "Hamilton Depression Scores")+
  theme(legend.position = c(0, 0),
        legend.justification = c(-0.1, -0.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(1.5, "cm"))

fixef(m1)

(Intercept)        Week 
   23.57704    -2.37705

ranef(m1)$ID %>% head()

    (Intercept)      Week
101    1.057202 -2.115138
103    3.670790 -0.483248
104    2.672755 -1.500882
105   -3.041339  0.226450
106    0.315424  1.025475
107   -0.614899 -0.429738

# only the first 6 participants

coef(m1)$ID %>% head()

    (Intercept)     Week
101     24.6342 -4.49219
103     27.2478 -2.86030
104     26.2498 -3.87793
105     20.5357 -2.15060
106     23.8925 -1.35157
107     22.9621 -2.80679

coef(m1)$ID %>% 
  ggplot(aes(x = Week,
             y = `(Intercept)`)) +
  geom_point(pch=1, size=rel(.5)) +
  geom_hline(yintercept = fixef(m1)["(Intercept)"],
             linetype = "dotted") +
  geom_vline(xintercept = fixef(m1)["Week"],
             linetype = "dotted") +
  labs(subtitle = "Estimated random effects",
       x = "Weekly Change in Depression (Slopes)",
       y = "Baseline Depression Level (Intercepts)")

5.2.1 Problem

Fit a model with random intercepts and slopes but the slopes and intercepts are uncorrelated. Compare models.

m1_2 <- lme4::lmer(Score ~ 1 + Week  + (Week -1|ID) + (1| ID), data = dtaL)

texreg::knitreg(list(m1, m1_2),
                single.row = TRUE,
                stars = numeric(0),
                custom.model.names = c("Random Intercepts And Random Slopes",
                                       "Slopes and Intercepts Are Uncorrelated"),
                caption = "MLM: Models fit w/REML",
                caption.above = TRUE,
                custom.note = "")

MLM: Models fit w/REML
	Random Intercepts And Random Slopes	Slopes and Intercepts Are Uncorrelated
(Intercept)	23.58 (0.55)	23.57 (0.52)
Week	-2.38 (0.21)	-2.38 (0.20)
AIC	2231.92	2232.34
BIC	2255.48	2251.98
Log Likelihood	-1109.96	-1111.17
Num. obs.	375	375
Num. groups: ID	66	66
Var: ID (Intercept)	12.94
Var: ID Week	2.13	1.80
Cov: ID (Intercept) Week	-1.48
Var: Residual	12.21	12.62
Var: ID.1 (Intercept)		10.63

-m1:在治療開始前，全體平均HD量表分數為23.58。時間每增加1單位，每位病人的HD量表分數會下降2.38。ICC為0.474(=12.94/(12.94+12.21+2.13))，表示在相同時間的兩個分數之間的相似性有47.4%。AIC為2231.92。

-m1_2:在治療開始前，全體平均HD量表分數為23.57。時間每增加1單位，每位病人的HD量表分數會下降2.38。ICC為0.504(=10.63/(10.63+12.62+1.8))，表示在相同時間的兩個分數之間的相似性有50.3%。AIC為2232.34，m1與m1_2的模型適配度相近。

anova(m1, m1_2, 
      model.names = c("Random Intercepts/Slopes", "Slopes and Intercepts Are Uncorrelated"),
      refit = TRUE) %>% 
  pander::pander(caption = "LRT: Significance of Random Slopes")

refitting model(s) with ML (instead of REML)

LRT: Significance of Random Slopes (continued below)
	npar	AIC	BIC	logLik	deviance	Chisq
Slopes and Intercepts Are Uncorrelated	5	2231	2251	-1111	2221	NA
Random Intercepts/Slopes	6	2231	2255	-1110	2219	2.336

	Df	Pr(>Chisq)
Slopes and Intercepts Are Uncorrelated	NA	NA
Random Intercepts/Slopes	1	0.1264

m2 <- lme4::lmer(Score ~ 1 + Week * Endogenous + 
                (1 + Week | ID), 
                data = dtaL)

interactions::interact_plot(m2,
                            pred = Week,
                            modx = Endogenous,
                            interval = TRUE,
                   main.title = "Time by Diagnosis Effect")+
  labs(x="Time (weeks since baseline)",
       y="Mean depression score")+
  theme(legend.position=c(.2,.2))

anova(m1, 
      m2, 
      model.names = c("Time only", 
                      "Time X Dx")) %>% 
  pander::pander(caption = "LRT: Significance of Diagnosis by Time")

refitting model(s) with ML (instead of REML)

LRT: Significance of Diagnosis by Time
	npar	AIC	BIC	logLik	deviance	Chisq	Df	Pr(>Chisq)
Time only	6	2231	2255	-1110	2219	NA	NA	NA
Time X Dx	8	2231	2262	-1107	2215	4.108	2	0.1282

5.3 Models with a quadratic time trend

It might be better to augment the data with a column variable coding the quadratic trend in number of weeks since baseline.

m3 <- lme4::lmer(Score ~ 1 + Week + I(Week^2) + 
                 (1 + Week + I(Week^2) | ID), 
                 data = dtaL,
                 control = lmerControl(optimizer = "optimx",    # get it to converge
                 calc.derivs = FALSE,
                 optCtrl = list(method = "nlminb",
                                starttests = FALSE,
                                kkt = FALSE)))

texreg::knitreg(list(m1,
                     m3),
                custom.model.names = c("Linear Trend",
                                       "QUadratic Trend"),
                single.row = TRUE,
                stars = numeric(0),
                caption = "Linear versus quadratic trend",
                caption.above = TRUE, 
                custom.note = "")

Linear versus quadratic trend
	Linear Trend	QUadratic Trend
(Intercept)	23.58 (0.55)	23.76 (0.56)
Week	-2.38 (0.21)	-2.63 (0.48)
Week^2		0.05 (0.09)
AIC	2231.92	2231.62
BIC	2255.48	2270.89
Log Likelihood	-1109.96	-1105.81
Num. obs.	375	375
Num. groups: ID	66	66
Var: ID (Intercept)	12.94	10.75
Var: ID Week	2.13	6.86
Cov: ID (Intercept) Week	-1.48	-1.03
Var: Residual	12.21	10.51
Var: ID I(Week^2)		0.20
Cov: ID (Intercept) I(Week^2)		-0.10
Cov: ID Week I(Week^2)		-0.97

-QUadratic Trend:在治療開始前，全體平均HD量表分數為23.76。時間每增加1單位，每位病人的HD量表分數會下降2.63；平方項增加1單位，每位病人的HD量表分數會增加0.05。ICC為0.371(=10.51/(10.51+0.2+6.86+10.75))，表示在相同時間的兩個分數之間的相似性有37.1%。Linear Trend與QUadratic Trend的AIC相近，暫時判斷兩者模型適配度沒有優劣。

anova(m1, m3)

Data: dtaL
Models:
m1: Score ~ 1 + Week + (1 + Week | ID)
m3: Score ~ 1 + Week + I(Week^2) + (1 + Week + I(Week^2) | ID)
   npar  AIC  BIC logLik deviance Chisq Df Pr(>Chisq)
m1    6 2231 2255  -1110     2219                    
m3   10 2228 2267  -1104     2208 11.39  4     0.0225

-m3的AIC為2228，比起m1，m3的模型適配度較佳。

fixef(m3)

(Intercept)        Week   I(Week^2) 
 23.7607640  -2.6332347   0.0516462

coef(m3)$ID |> head()

    (Intercept)      Week I(Week^2)
101     25.1997 -5.306768  0.153506
103     27.5675 -3.485669  0.128337
104     26.0138 -3.088015 -0.179926
105     20.9825 -2.936235  0.162536
106     23.6232 -0.721621 -0.145821
107     22.6693 -1.943044 -0.187271

# subject 115
fun_115 <- function(Week){
  coef(m3)$ID["115", "(Intercept)"] +
  coef(m3)$ID["115", "Week"] * Week +
  coef(m3)$ID["115", "I(Week^2)"] * Week^2
}

# subject 610
fun_610 <- function(Week){
  coef(m3)$ID["610", "(Intercept)"] +
  coef(m3)$ID["610", "Week"] * Week +
  coef(m3)$ID["610", "I(Week^2)"] * Week^2
}

dtaL %>% 
  dplyr::mutate(pred_fixed = predict(m3, re.form = NA)) %>% # fixed effects only
  dplyr::mutate(pred_wrand = predict(m3)) %>%               # fixed and random effects together
  ggplot(aes(x = Week,
             y = Score,
             group = ID)) +
  stat_function(fun = fun_115) +          # add cure for ID = 115
  stat_function(fun = fun_610) +          # add cure for ID = 610
  geom_line(aes(y        = pred_fixed),
                color    = "blue",
                size     = 1.25)  +
  labs(x = "Weeks Since Baseline",
       y = "Hamilton Depression Scores",
       subtitle = "Marginal Mean show in thicker blue\nBLUPs for two of the participant in thinner black")+
  theme_minimal()+
  theme(legend.position = c(0, 0),
        legend.justification = c(-0.1, -0.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(1.5, "cm"))

-115病人與610病人的HD分數變化曲線，藍色線為全體變化趨勢。

dtaL %>% 
  dplyr::mutate(pred_fixed = predict(m3, re.form = NA)) %>% # fixed effects only
  dplyr::mutate(pred_wrand = predict(m3)) %>%               # fixed and random effects together
  ggplot(aes(x = Week,
             y = Score,
             group = ID)) +
  geom_line(aes(y        = pred_wrand,
                color    = "BLUP",
                size     = "BLUP",
                linetype = "BLUP")) +
  geom_line(aes(y        = pred_fixed,
                color    = "Marginal",
                size     = "Marginal",
                linetype = "Marginal")) +
  scale_color_manual(name   = "Type of Prediction",
                     values = c("BLUP"     = "gray50",
                                "Marginal" = "blue"))  +
  scale_size_manual(name    = "Type of Prediction",
                    values = c("BLUP"      = .5,
                                "Marginal" = 1.25))  +
  scale_linetype_manual(name   = "Type of Prediction",
                        values = c("BLUP"     = "longdash",
                                   "Marginal" = "solid"))  +
  labs(x = "Weeks Since Baseline",
       y = "Hamilton Depression Scores")+
  theme_minimal()+
  theme(legend.position = c(0, 0),
        legend.justification = c(-0.1, -0.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(1.5, "cm"))

6 References

Hedeker, D. (2004). An introduction to growth modeling. In D. Kaplan (Ed.), Quantitative Methodology for the Social Sciences. Thousand Oaks CA: Sage Publications. pdf

Reisby, N., Gram, L.F., Bech, P. et al. (1977). Imipramine: Clinical effects and pharmacokinetic variability. Psychopharmacology, 54, 263-272.

Riesby

Huang En-Li

2022-11-30

1 Introduction

2 Data

3 Summary statistics

3.1 Missing data

3.1.1 Problem

3.2 Covariances and correlations

3.2.1 Problem

3.2.2 Problem

3.2.3 Problem

4 Visualization

5 Models

5.1 Random intercepts only

5.1.1 Fitted values

5.2 Random slopes and intercepts model

5.2.1 Problem

5.3 Models with a quadratic time trend

6 References