12-example_hedeker4_reisby.Rmd

# MLM, Longitudinal: Hedeker & Gibbons - Depression

```{r, echo=FALSE}
library(pander)

panderOptions('digits', 2)
panderOptions('round', 2)
panderOptions('keep.trailing.zeros', TRUE)
panderOptions('keep.line.breaks', TRUE)
```

```{r, warning=FALSE, error=FALSE, message=FALSE, results='hide'}
library(tidyverse)    # all things tidy
library(pander)       # nice looking general tabulations
library(furniture)    # nice table1() descriptions
library(texreg)       # Convert Regression Output to LaTeX or HTML Tables

library(psych)        # contains some useful functions, like headTail()
library(sjstats)      # ICC calculations
library(effects)      # Estimated Marginal Means
library(performance)
library(VIM)          # Visualization and Imputation of Missing Values 
library(naniar)       # Summaries and Visualizations for Missing Data
library(corrplot)     # Visualize correlation matrix

library(lme4)         # Linear, generalized linear, & nonlinear mixed models
library(lmerTest)
library(optimx)
```

## Background

Starting in chapter 4, [@hedeker2006] details analysis of a psychiatric study described by [@reisby1977]. This study focuses on the relationship between Imipramine (IMI) and Desipramine (DMI) plasma levels and clinical response in 66 depressed inpatients (37 endogenous and 29 non-endogenous).

> Note: The IMI and DMI measures were only taken in the later weeks, but are not used here.

* `hamd` Hamilton Depression Scores (HD)


Independent or predictor variables:  

* `endog` Depression Diagnosis 
+ endogenous  
+ non-endogenous/reactive 

* IMI (imipramine) drug-plasma levels (µg/l)     
+ antidepressant given 225 mg/day, weeks 3-6   

* DMI (desipramine) drug-plasma levels (µg/l)
+ metabolite of imipramine

```{r}
data_raw <- read.table("data/riesby.dat.txt") %>% 
  dplyr::rename(id    = "V1",
                hamd  = "V2",
                endog = "V5",
                week  = "V4") %>%
  dplyr::select(-V3, -V6)
```



```{r}
data_raw %>% 
  psych::headTail(top = 11, bottom = 8)
```



### Long Format

```{r}
data_long <- data_raw %>%                        
  dplyr::filter(!is.na(hamd)) %>%                                # remove NA or missing scores
  dplyr::mutate(id = factor(id)) %>%                             # declare grouping var a factor
  dplyr::mutate(endog = factor(endog,                            # attach labels to a grouping variable
                               levels = c(0, 1),                 # order of the levels should match levels
                               labels = c("Reactive",            # order matters!
                                          "Endogenous"))) %>% 
  dplyr::mutate(hamd = as.numeric(hamd)) %>% 
  dplyr::select(id, week, endog, hamd) %>%                       # select the order of variables to include
  dplyr::arrange(id, week)                                       # sort observations  
```



```{r}
data_long %>% 
  psych::headTail(top = 11, bottom = 8)
```



### Wide Format

```{r}
data_wide <- data_long %>%     # save the dataset as a new name
  tidyr::pivot_wider(names_from = week,
                     names_prefix = "hamd_",
                     values_from = hamd)
```

Notice the missing values, displayed as `NA`.

```{r}
data_wide %>% 
  psych::headTail()
```

## Exploratory Data Analysis

### Diagnosis Group

```{r, results="asis"}
data_wide %>% 
  furniture::table1("Depression Type" = endog,
                    output = "markdown")
```

### Missing Data & Attrition

Plot the amount of missing values and the amount of each patter of missing values.

#### `VIM` package

```{r}
data_wide %>% 
  VIM::aggr(numbers = TRUE,   # shows the number to the far right
            prop = FALSE)     # shows counts instead of proportions
```

#### `naniar` package

```{r}
data_wide %>% 
  naniar::gg_miss_upset(sets = c("hamd_5_NA",
                                 "hamd_4_NA",
                                 "hamd_3_NA",
                                 "hamd_2_NA",
                                 "hamd_1_NA",
                                 "hamd_0_NA"),
                        keep.order = TRUE) 
```

### Depression Over Time, by Group

#### Table of Means

> Default = COMPLETE CASES ONLY!!!

Note - the sample size at each time point is the same, but subjects are only included if they have data at every time point

```{r}
data_wide %>%          
  dplyr::group_by(endog) %>% 
  furniture::table1("Baseline" = hamd_0,
                    "Week 1" = hamd_1,
                    "Week 2" = hamd_2,
                    "Week 3" = hamd_3,
                    "Week 4" = hamd_4,
                    "Week 5" = hamd_5,
                    total = TRUE,
                    test = TRUE,
                    na.rm = TRUE,              # default: COMPLETE CASES ONLY!!!!!
                    digits = 2,
                    output = "markdown",
                    caption = "Hamilton Depression Scores Across Time, by Depression Type for Participants with 6 Complete Weeks")           
```

> Specify All data:

note - that the smaple sizes will be different for each time point

```{r, results='asis'}
data_wide %>%                
  dplyr::group_by(endog) %>% 
  furniture::table1("Baseline" = hamd_0,
                    "Week 1" = hamd_1,
                    "Week 2" = hamd_2,
                    "Week 3" = hamd_3,
                    "Week 4" = hamd_4,
                    "Week 5" = hamd_5,
                    total = TRUE,
                    test = TRUE,
                    na.rm = FALSE,              # default: COMPLETE CASES ONLY!!!!!
                    digits = 2,
                    output = "markdown",
                    caption = "Hamilton Depression Scores Across Time, by Depression Type for All Participants")  
```

#### Using the LONG formatted dataset

Each person's data is stored on multiple lines, one for each time point.

> FOR ALL DATA!

```{r}
data_sum_all <- data_long %>% 
  dplyr::group_by(endog, week) %>%                # specify the groups
  dplyr::filter(!is.na(hamd)) %>% 
  dplyr::summarise(hamd_n    = n(),                      # count of valid scores
                   hamd_mean = mean(hamd),               # mean score
                   hamd_sd   = sd(hamd),                 # standard deviation of scores
                   hamd_sem  = hamd_sd / sqrt(hamd_n))   # stadard error for the mean of scores

data_sum_all %>% 
  pander::pander()
```

> FOR COMPLETE CASES ONLY!!!

```{r}
data_sum_cc <- data_long %>% 
  dplyr::group_by(id) %>%               # group-by participant
  dplyr::filter(!is.na(hamd)) %>% 
  dplyr::mutate(person_vsae_n = n()) %>%    # count the number of valid VSAE scores
  dplyr::filter(person_vsae_n == 6) %>%     # restrict to only thoes children with 5 valid scores
  dplyr::group_by(endog, week) %>%                # specify the groups
  dplyr::summarise(hamd_n    = n(),                      # count of valid scores
                   hamd_mean = mean(hamd),               # mean score
                   hamd_sd   = sd(hamd),                 # standard deviation of scores
                   hamd_sem  = hamd_sd / sqrt(hamd_n))   # stadard error for the mean of scores

data_sum_cc %>% 
  pander::pander()
```

#### Person-Profile Plot or Spaghetti Plot

Use the data in LONG format.

A scatterplot of all observations of depression scores over time, joining the dots of each individual's data.

> NOTE: Not all lines have a point for every week!

```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd)) +
  geom_point() +
  geom_line(aes(group = id)) +   # join points that belong to the same "id"
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)") 
```

```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd,
             color = endog)) +    # color points and lines by the "endog" variable
  geom_line(aes(group = id)) +
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)") 
```

```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd)) +
  geom_line(aes(group = id)) +
  facet_grid( ~ endog) +           # side-by-side pandels by the "endog" variable
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)") 
```

```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week %>% factor(),
             y = hamd)) +
  geom_boxplot() +                   # compare center and spread
  facet_grid( ~ endog) +  
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)")         
```

```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week %>% factor(),
             y = hamd)) +
  geom_violin(fill = "gray75") +                   # similar to boxplots to show distribution
  stat_summary() +  
  stat_summary(aes(group = "endog"),
               fun = mean,
               geom = "line") +
  facet_grid( ~ endog)   +  
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)")          
```


```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd,
             color = endog)) +
  stat_summary() +
  stat_summary(aes(group = endog,
                   linetype = endog),
               fun = mean,
               geom = "line") +
  scale_linetype_manual(values = c("solid", "longdash")) +
  theme(legend.position = "bottom",
        legend.key.width = unit(2, "cm")) +  
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)") 
```



```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd)) +
  geom_line(aes(group = id)) +
  facet_grid( ~ endog) +
  geom_smooth() +                     # DEFAULTS: method = "loess", se = TRUE, color = "blue"
  geom_smooth(method = "lm",
              se = FALSE,
              color = "hot pink")+  
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)") 
```

```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd)) +
  facet_grid( ~ endog) +
  geom_smooth(aes(color = "Flexible"),
              method = "loess",
              se = FALSE,) +
  geom_smooth(aes(color = "Linear"),
              method = "lm",
              se = FALSE) +
  scale_color_manual(name = "Smoother Type: ",
                     values = c("Flexible" = "blue", 
                                "Linear"   = "red")) +
  theme_bw() + 
  theme(legend.position = "bottom") +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)") 
```

```{r}
data_long %>% 
  ggplot(aes(x = week,
             y = hamd,
             group = endog,
             linetype = endog)) +
  geom_smooth(method = "loess",
              color = "black",
              alpha = .25) +
  theme_bw() +
  scale_linetype_manual(values = c("solid", "longdash")) +
  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")) +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Scale",
       linetype = "Type of Depression")
```

```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd,
             group = endog,
             linetype = endog)) +
  geom_smooth(method = "loess",
              color = "black",
              alpha = .25) +
  theme_bw() +
  scale_linetype_manual(values = c("solid", "longdash")) +
  theme(legend.position = c(1, 1),
        legend.justification = c(1.1, 1.1),
        legend.background = element_rect(color = "black"),
        legend.key.width = unit(2, "cm")) +
  labs(x = "Weeks Since Baseline",
       y = "EsStimated Hamilton Depression Scale",
       linetype = "Type of Depression")
```


```{r}
data_long %>% 
  dplyr::filter(!is.na(hamd)) %>% 
  ggplot(aes(x = week,
             y = hamd,
             group = endog,
             linetype = endog)) +
  geom_smooth(method = "lm",
              color = "black",
              alpha = .25) +
  theme_bw() +
  scale_linetype_manual(values = c("solid", "longdash")) +
  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")) +
  labs(x = "Weeks Since Baseline",
       y = "EStimated Hamilton Depression Scale",
       linetype = "Type of Depression")
```


## Patterns in the Outcome Over Time

### Variance-Covariance

#### Full Matrix

* Variances are down the diagonal    
+ Increasing variance over time violates the ANOVA assumption of *homogeity of variance*

```{r}
data_wide %>% 
  dplyr::select(starts_with("hamd_")) %>%  # just the outcome(s)
  cov(use = "complete.obs")  %>%           # covariance matrix, LIST-wise deletion
  round(3)
```

```{r}
data_wide %>% 
  dplyr::select(starts_with("hamd_")) %>%  # just the outcome(s)
  cov(use = "pairwise.complete.obs")  %>%  # covariance matrix, PAIR-wise deletion
  round(3)
```

#### Just Variances

Notice the variance in scores increases over time, which is seen in the side-by-side boxplots.

```{r}
data_wide %>% 
  dplyr::select(starts_with("hamd_")) %>%  # just the outcome(s)
  cov(use = "pairwise.complete.obs")  %>%  # covariance matrix, PAIR-wise deletion
  diag()                                   #  extracts just the variances
```

### Correlation

#### Full Matrix

Pairwise relationships are easier to eye-ball magnitude when presented as correlations, rather than covariances, due to the relative scale.

```{r}
data_wide %>% 
  dplyr::select(starts_with("hamd_")) %>% # just the outcome(s)
  cor(use = "complete.obs") %>%           # correlation matrix - LIST-wise deletion
  round(2)                                
```

```{r}
data_wide %>% 
  dplyr::select(starts_with("hamd_")) %>% # just the outcome(s)
  cor(use = "pairwise.complete.obs") %>%   # correlation matrix - PAIR-wise deletion
  round(2)                                
```

#### Visualization

Looking for patterns is always easier with a plot. All RM or mixed ANOVA assume sphericity or compound symmetry, meaning that all the correlations in the matrix would be the same. This is not the case for these data. Instead we see a classic pattern of corralary decay. Measures taken close in time, say 1 week apart, exhibit the highest degree of correlation. The farther apart in time that two measures are taken, the less correlated they are. Note that the adjacent measures become more highly correlated, too. This can be due to attrition; later time points having a smaller sample size.

```{r}
data_wide %>% 
  dplyr::select(starts_with("hamd_")) %>% # just the outcome(s)
  cor(use = "pairwise.complete.obs") %>%   # correlation matrix
  corrplot::corrplot.mixed(upper = "ellipse")
```

### For Each Group

It can be a good ideal to investigate if the groups exhibit a similar pattern in correlation.

> Reactive Depression

```{r}
data_wide %>% 
  dplyr::filter(endog == "Reactive") %>%    # filter observations for the REACTIVE group
  dplyr::select(starts_with("hamd_")) %>% 
  cor(use = "pairwise.complete.obs") %>%   
  corrplot::corrplot.mixed(upper = "ellipse")
```

> Endogenous Depression

```{r}
data_wide %>% 
  dplyr::filter(endog == "Endogenous") %>%    # filter observations for the Endogenous group
  dplyr::select(starts_with("hamd_")) %>%
  cor(use = "pairwise.complete.obs") %>%   
  corrplot::corrplot.mixed(upper = "ellipse")
```

## MLM - Null or Emptly Models

### Fit the model

> Random Intercepts, with Fixed Intercept and Time Slope (i.e. Trend)....\@hedeker2006 section 4.3.5, starting on page 55

Since this situation deals with longitudinal data, it is more appropriate to start off including the time variable in the null model as a fixed effect only.

```{r}
fit_lmer_week_RI_reml <- lmerTest::lmer(hamd ~ week + (1|id), 
                                        data = data_long,
                                        REML = TRUE)
```

### Table of Parameter Estimates




```{r mlmRI, results='asis'}
texreg::screenreg(fit_lmer_week_RI_reml, 
                single.row = TRUE,
                caption = "MLM: Random Intercepts Null Model fit w/REML",
                caption.above = TRUE,
                custom.note = "Reproduction of Hedeker's table 4.3 on page 55, except using REML here instead of ML")
```

On average, patients start off with HDRS scores of `r fixef(fit_lmer_week_RI_reml)["(Intercept)"] %>% round(2)` and then change by `r fixef(fit_lmer_week_RI_reml)["week"] %>% round(2)` points each week. This weekly improvement of about 2 points a week is statistically significant via the Wald test.

### Estimated Marginal Means Plot

> Multilevel model on page 55 [@hedeker2006]

$$
\hat{y} = `r fixef(fit_lmer_week_RI_reml)["(Intercept)"] %>% round(3)` + 
`r fixef(fit_lmer_week_RI_reml)["week"] %>% round(3)` week
$$

The fastest way to plot a model is to use the `sjPlot::plot_model()` function.

> Note: you can't use `interactions::interact_plot()` since there is only one predictor (i.e. independent variable) in this model.

```{r}
sjPlot::plot_model(fit_lmer_week_RI_reml,
                   type = "pred",
                   terms = c("week"))
```

### Estimated Marginal Means and Emperical Bayes Plots

With a bit more code we can plot not only the **marginal model** (fixed effects only), but add the Best Linear Unbiased Predictions (**BLUP**s) or person-specific specific models (both fixed and random effects).

```{r}
data_long %>% 
  dplyr::mutate(pred_fixed = predict(fit_lmer_week_RI_reml, 
                                     re.form = NA,           # fixed effects only
                                     newdata = .)) %>% 
  dplyr::mutate(pred_wrand = predict(fit_lmer_week_RI_reml,  # fixed and random effects both
                                     newdata = .)) %>%              
  ggplot(aes(x = week,
             y = hamd,
             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")) +
  theme_bw() +
  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"))  +
  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")) +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)")
```

Notice that in this model, all the BLUPs are parallel. That is because we are only letting the intercept vary from person-to-person while keeping the effect of time (slope) constant.

> Reproduce Table 4.4 on page 55 [@hedeker2006]

One way to judge a model is to compare the estimated means to the observed means to see how accuratedly they are represented by the model. This excellent fit of the estimated marginal means to the observed data supports the hypothesis that the change in depression across time is LINEAR.

```{r}
obs <- data_long %>% 
  dplyr::group_by(week) %>% 
  dplyr::summarise(observed = mean(hamd, na.rm = TRUE)) 

effects::Effect(focal.predictors = "week",
                mod = fit_lmer_week_RI_reml,
                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")
```

### Intra-individual Correlation (ICC)

```{r}
performance::icc(fit_lmer_week_RI_reml)
```

**Interpretation** Just less than a third ($\rho_c = .32$) in baseline depression is explained by person-to-person differences. Thus, subjects display considerable heterogeneity in depression levels.

This value of .46is an oversimplification of the correlation matrix above and may be thought of as the expected correlation between two randomly drawn weeks for any given person.

```{r}
performance::r2(fit_lmer_week_RI_reml)
```

**Interpretation**  Linear growth accounts for 31% of the variance in Hamilton Depression Scores across the six weeks.  Linear growth **AND** person-to-person differences account for a total 63% of this variance.

This value of 46% is an oversimplification of the correlation matrix above and may be thought of as the expected correlation between two randomly drawn weeks for any given person.

```{r}
performance::r2(fit_lmer_week_RI_reml)
```

> Note: The marginal $R^2$ considers only the variance of the fixed effects, while the conditional $R^2$ takes both the fixed and random effects into account. The random effect variances are actually the mean random effect variances, thus the $R^2$ value is also appropriate for mixed models with random slopes or nested random effects (see Johnson 2014).

**Interpretation** Linear growth accounts for 31% of the variance in Hamilton Depression Scores across the six weeks. Linear growth AND person-to-person differences account for a total 63% of this variance.

> Note: The marginal R^2 considers only the variance of the fixed effects, while the conditional R^2 takes both the fixed and random effects into account. The random effect variances are actually the mean random effect variances, thus the R^2 value is also appropriate for mixed models with random slopes or nested random effects (see Johnson 2014).



### Compare to the Single-Level Null: No Random Effects

> Simple Linear Regression, [@hedeker2006]

To compare, fit the single level regression model

```{r}
fit_lm_week_ml <- lm(hamd ~ week,
                     data = data_long)
```

#### Table of Parameter Estimates



```{r olsvsmlm, results='asis'}
texreg::screenreg(list(fit_lm_week_ml, fit_lmer_week_RI_reml),
                custom.model.names = c("Single-Level", "Multilevel"),
                single.row = TRUE,
                caption = "MLM: Longitudinal Null Models",
                caption.above = TRUE,
                custom.note = "The singel-level model treats are observations as being independent and unrelated to each other, even if they were made on the same person.")
```

For the multilevel model, the Wald tests indicated the fixed intercept is significant *(no surprised that the depressions scores are not zero at baseline)*. More of note is the significance of the fixed effect of time. This signifies that depression scores are declining over time. On average, patients are improving (Hamilton Depression Scores get smaller) across time, by an average of 2.4'ish points a week.

For the multilevel model, the Wald tests indicated the fixed intercept is significant (no surprised that the depressions scores are not zero at baseline). More of note is the significance of the fixed effect of time. This signifies that depression scores are declining over time. On average, patients are improving (Hamilton Depression Scores get smaller) across time, by an average of 2.4’ish points a week.



#### Residual Variance

Note, the fixed estimates are very similar for the two models, but the standard errors are different. Additionally, whereas the single-level regression lumps all remaining variance together ($\sigma^2$), the multilevel model seperates it into within-subjects ($\sigma^2_{u0}$ or $\tau_{00}$) and between-subjects variance ($\sigma^2_{e}$ or $\sigma^2$).

```{r}
sigma(fit_lm_week_ml)^2
```

```{r}
lme4::VarCorr(fit_lmer_week_RI_reml) %>% # in longitudinal data, a group of observations = a participant or person
  print(comp = c("Variance", "Std.Dev"))
```



> "One statistician's error term is another's career!"
>
> [@hedeker2006], page 56

## MLM: Add Random Slope for Time (i.e. Trend)

### Fit the Model

```{r}
fit_lmer_week_RIAS_reml <- lmerTest::lmer(hamd ~ week + (week|id), #     MLM-RIAS
                                          data = data_long,
                                          REML = TRUE)
```




```{r mlmRS, results='asis'}
texreg::screenreg(list(fit_lmer_week_RI_reml, 
                     fit_lmer_week_RIAS_reml),
                custom.model.names = c("Random Intercepts",
                                       " And Random Slopes"),
                single.row = TRUE,
                caption = "MLM: Null models fit w/REML",
                caption.above = TRUE,
                custom.note = "Hedeker table 4.4 on page 55 and table 4.5 on page 58, except using REML here instead of ML")
```



Visually, we can see that the unexplained or residual variance is less (12.21 vs 19.10) for the model that includes person-specific slopes (trajectories over time).

Note: the negative covariance between random intercepts and random slopes ($\sigma_{u01} = \tau_{01} = -1.48$):


> "This suggests that patients who are initially more depressed (i.e. greater intercepts) improve at a greater rate (i.e. more pronounced negative slopes). An alternative explainatio, though,is that of a floor effect due to the HDRS rating scale. Simply put, patients with less depressed intitial scores have a more limited range of lower scores than those with higher initial scores."
>
> [@hedeker2006], page 58

### Likelihood Ratio Test

```{r}
anova(fit_lmer_week_RI_reml, fit_lmer_week_RIAS_reml, 
      model.names = c("RI", "RIAS"),
      refit = FALSE) %>% 
  pander::pander(caption = "LRT: Assess Significance of Random Slopes")
```

Including the random slope for time significantly improved the model fit via the formal **Likelihood Ratio Test**. This rejects the assumption of compound symmetry.


```{r}
performance::compare_performance(fit_lmer_week_RI_reml, 
                                 fit_lmer_week_RIAS_reml,
                                 rank = TRUE)
```




### Estimated Marginal Means Plot

```{r}
sjPlot::plot_model(fit_lmer_week_RIAS_reml,
                   type = "pred",
                   terms = c("week"))
```

Adding the random slopes didn't change the estimates for the fixed effects much.

### Estimated Marginal Means and Emperical Bayes Plots

```{r}
data_long %>% 
  dplyr::mutate(pred_fixed = predict(fit_lmer_week_RIAS_reml, 
                                     re.form = NA,               # fixed effects only
                                     newdata = .)) %>% 
  dplyr::mutate(pred_wrand = predict(fit_lmer_week_RIAS_reml,    # fixed and random effects together
                                     newdata = .)) %>%               
  ggplot(aes(x = week,
             y = hamd,
             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")) +
  theme_bw() +
  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"))  +
  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")) +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score")
```

BLUPs are also refered to as **Empirical Bayes Estimates** and may be extracted from a model fit. In this cases there will be a specific intercept (`(Intercept)`) and time slope (`week`) for each individual or person (`id`).

#### Fixed Effects

Marginal Model = within-subject effects

```{r}
fixef(fit_lmer_week_RIAS_reml)
```

#### Random Effects

between-subjects effects

```{r}
ranef(fit_lmer_week_RIAS_reml)$id %>% 
  head()                                 # only the first 6 participants
```

#### BLUPs or Empirical Bayes Estimates

fixed effects + random effects

```{r}
coef(fit_lmer_week_RIAS_reml)$id %>% 
  head()                                 # only the first 6 participants
```

We can create a scatterplot of these to see the correlation between them.

```{r}
coef(fit_lmer_week_RIAS_reml)$id %>% 
  ggplot(aes(x = week,
             y = `(Intercept)`)) +
  geom_point() +
  geom_hline(yintercept = fixef(fit_lmer_week_RIAS_reml)["(Intercept)"],
             linetype = "dashed") +
  geom_vline(xintercept = fixef(fit_lmer_week_RIAS_reml)["week"],
             linetype = "dashed") +
  geom_smooth(method = "lm") +
  labs(title = "Hedeker's Figure 4.4 on page 59",
       subtitle = "Reisby data: Estimated random effects",
       x = "Week Change in Depression",
       y = "Baseline Depression Level") +
  theme_bw()

```

## MLM: Coding of Time

So far we have used the variable `week` to denote time as weeks since baseline = `week` $\in 0, 1, 2, 3, 4, 5$.

But We could CENTER week at the middle of the study (`week` = 2.5).



### Fit the Model

```{r}
fit_lmer_week_RIAS_reml_wc <- lme4::lmer(hamd ~ I(week-2.5) + (I(week-2.5)|id), #     MLM-RIAS
                                         data = data_long,
                                         REML = TRUE)
```

### Table of Parameter Estimates


```{r mlmwc, results='asis'}
texreg::screenreg(list(fit_lmer_week_RIAS_reml, 
                     fit_lmer_week_RIAS_reml_wc),
                custom.model.names = c("Random Intercepts",
                                       " And Random Slopes"),
                single.row = TRUE,
                caption = "MLM: Null models fit w/REML",
                caption.above = TRUE,
                custom.note = "Hedeker table table 4.5 on page 58 and table 4.6 on page 61, except using REML here instead of ML")
```

* Unchanged    
+ model fit: AIC, BIC, -2LL, residual variance    
+ fixed effect of week    
+ variance for random intercepts     

* Changed    
+ fixed intercept    
+ variance for random slopes   
+ covariance between random intercepts and random slopes

    -   model fit: AIC, BIC, -2LL, residual variance\
    -   fixed effect of week\
    -   variance for random intercepts

-   Changed

    -   fixed intercept\
    -   variance for random slopes\
    -   covariance between random intercepts and random slopes

### Estimated Marginal Means and Emperical Bayes Plots

```{r}
data_long %>% 
  dplyr::mutate(pred_fixed = predict(fit_lmer_week_RIAS_reml_wc, 
                                     re.form = NA, 
                                     newdata = .)) %>%               # fixed effects only
  dplyr::mutate(pred_wrand = predict(fit_lmer_week_RIAS_reml_wc, 
                                     newdata = .)) %>%               # fixed and random effects together
  ggplot(aes(x = week,
             y = hamd,
             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")) +
  theme_bw() +
  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"))  +
  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")) +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Score (HD)")
```

**Again, centering time doesn't change the interpretation at all, since there are no interactions.**

## MLM: Effect of DIagnosis on Time Trends (Fixed Interaction)

The researcher specifically wants to know if the trajectory over time differs for the two types of depression. This translates into a fixed effects interaction between time and group.

> Start by comapring random intercepts only (RI) to a random intercetps and slopes (RIAS) model.

### Fit the Models

```{r}
fit_lmer_week_RIAS_ml <- lmerTest::lmer(hamd ~ week + (week|id), 
                                        data = data_long,
                                        REML = FALSE)

fit_lmer_wkdx_RIAS_ml <- lmerTest::lmer(hamd ~ week*endog + (week|id), 
                                        data = data_long,
                                        REML = FALSE)
```

### Estimated Marginal Meanse Plot

```{r}
interactions::interact_plot(fit_lmer_wkdx_RIAS_ml,
                            pred = week,
                            modx = endog,
                            legend.main = "Type of Depression",
                            interval = TRUE,
                            main.title = "Hedeker's Table 4.7 on page 64") +
  theme_bw() +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Marginal Mean\nHamilton Depression Scores (HD)") +
  theme(legend.background = element_rect(color = "black"),
        legend.position = c(0, 0),
        legend.justification = c(-0.1, -0.1),
        legend.key.width = unit(1.75, "cm"))
```





### Likelihood Ratio Test

```{r}
anova(fit_lmer_week_RIAS_ml, 
      fit_lmer_wkdx_RIAS_ml, 
      model.names = c("Just Time", 
                      "Time X Dx")) %>% 
  pander::pander(caption = "LRT: Assess Significance of Diagnosis Moderation of Trend")
```


```{r}
performance::compare_performance(fit_lmer_week_RIAS_ml, 
                                 fit_lmer_wkdx_RIAS_ml,
                                 rank = TRUE)
```



The more complicated model (including the moderating effect of diagnosis) is NOT supported.


The more complicated model (including the moderating effect of diagnosis) is **NOT** supported.

## MLM: Quadratic Trend

### Fit the Model

```{r}
fit_lmer_quad_RIAS_ml <- lme4::lmer(hamd ~ week + I(week^2) + (week+ I(week^2)|id), 
                                    data = data_long,
                                    REML = FALSE,
                                    control = lmerControl(optimizer = "optimx",    # get it to converge
                                                          calc.derivs = FALSE,
                                                          optCtrl = list(method = "nlminb",
                                                                         starttests = FALSE,
                                                                         kkt = FALSE))) 
```

### Table of Parameter Estimates


```{r mlmquad, results='asis'}
texreg::screenreg(list(fit_lmer_week_RIAS_ml,
                     fit_lmer_quad_RIAS_ml),
                custom.model.names = c("Linear Trend",
                                       "QUadratic Trend"),
                single.row = TRUE,
                caption = "MLM: RIAS models fit w/ML",
                caption.above = TRUE,
                custom.note = "Hedeker table 4.5 on page 58 and table 5.1 on page 84")
```

### Likelihood Ratio Test

```{r}
anova(fit_lmer_week_RIAS_ml, fit_lmer_quad_RIAS_ml)
```

```{r}
performance::compare_performance(fit_lmer_week_RIAS_ml, 
                                 fit_lmer_quad_RIAS_ml,
                                 rank = TRUE)
```



Even though the Wald test did not find the quadratic fixed time trend to be significant at the population level (marginal), the LRT and Bayes Factor both find that including the quadratic terms improves the model’s fit.




### Estimated Marginal Means Plot

```{r}
fixef(fit_lmer_quad_RIAS_ml)
```

```{r}
sjPlot::plot_model(fit_lmer_quad_RIAS_ml,
                   type = "pred",
                   terms = "week")
```

At the population level, the curviture is very slight.

### BLUPs or Emperical Bayes Estimates

```{r}
coef(fit_lmer_quad_RIAS_ml)$id
```

For Illustration, two cases have been hand selected: `id` = 115 and 610.

```{r}
fun_115 <- function(week){
  coef(fit_lmer_quad_RIAS_ml)$id["115", "(Intercept)"] +
    coef(fit_lmer_quad_RIAS_ml)$id["115", "week"] * week +
    coef(fit_lmer_quad_RIAS_ml)$id["115", "I(week^2)"] * week^2
}


fun_610 <- function(week){
  coef(fit_lmer_quad_RIAS_ml)$id["610", "(Intercept)"] +
    coef(fit_lmer_quad_RIAS_ml)$id["610", "week"] * week +
    coef(fit_lmer_quad_RIAS_ml)$id["610", "I(week^2)"] * week^2
}
```

```{r, fig.cap="Two Example BLUPS for two different participants"}
data_long %>% 
  dplyr::mutate(pred_fixed = predict(fit_lmer_quad_RIAS_ml, 
                                     re.form = NA,
                                     newdata = .)) %>% # fixed effects only
  dplyr::mutate(pred_wrand = predict(fit_lmer_quad_RIAS_ml,
                                     newdata = .)) %>%               # fixed and random effects together
  ggplot(aes(x = week,
             y = hamd,
             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)  +
  theme_bw() +
  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")) +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Scores (HD)",
       title = "Similar to Hedeker's Figure 5.3 on page 84",
       subtitle = "Marginal Mean show in thicker blue\nBLUPs for two of the participant in thinner black")
```

These two individuals have quite different curvatures and illustrated how this type of curvatures in person-specific trajectories may end up cancelling each other out to arrive at a fairly linear marginal model.

### Estimated Marginal Means and Emperical Bayes Plots


Note: although the BLUPs are shown for all participants, the predictions are just connects and are therefore slightly jagged and now smoother like the lines on the plot above.

```{r, fig.cap = "EStimated curvilinear trends"}
data_long %>% 
  dplyr::mutate(pred_fixed = predict(fit_lmer_quad_RIAS_ml, 
                                     re.form = NA,
                                     newdata = .)) %>% # fixed effects only
  dplyr::mutate(pred_wrand = predict(fit_lmer_quad_RIAS_ml,
                                     newdata = .)) %>%               # fixed and random effects together
  ggplot(aes(x = week,
             y = hamd,
             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")) +
  theme_bw() +
  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"))  +
  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")) +
  labs(x = "Weeks Since Baseline",
       y = "Estimated Hamilton Depression Scores",
       title = "Hedeker's Figure 5.4 on page 85")
```

At the person-level, the curvature is very diverse (heterogeneous). Some individuals have accelerating downward tend while other have accelerating upward trends.

The improvement that the curvi-linear model provides in describing change across time is perhaps modest.