```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Overview

We will be using data from the same study (PANDAS) as from last week's challenge to practice making longitudinal multilevel models. As a reminder, we are not focusing on running models yet. These data include a variable for ID and gender, and then assess self-worth (SPPA_SWORTH), communication quality with mothers (MPACS), and generalized anxiety (SCARED_GAD). Self-worth, communication and GAD were assessed every 6-months at 6 different timepoints. The csv file is the same as the lesson mlm_datacleaning.csv and can be found here: https://columbia-sipps.github.io/workshops/coding_advanced/

Our question for today will be *how did self worth change over time?* We'll use longitudinal multilevel models in `rstanarm` to approach this

# Load in packages and data

```{r,warning=FALSE, results='hide', message=FALSE}
library(tidyverse)
library(rstanarm)
library(bayesplot)
library(tidybayes)

d = read_csv('https://columbia-sipps.github.io/workshop_files/coding-adv-04/mlm-data-cleaning.csv')
```

# Setup 

## Filter data for only IDs up to 100 (just for demo purposes of a smaller dataset)

```{r}
d <- d %>% 
  filter(ID <= 100)
```

## Select only the variables we will be plotting

Being: id, gender, and all of the self worth scores
```{r}
d <- d %>% 
  select(ID, contains("SPPA"))
```


## Rename your variables 

Rename them so that they are easier to understand 
```{r}
d <- d %>% 
# rename vectors
  rename(
    id = 1,                   
    sworth_1 = 2,   
    sworth_2 = 3,   
    sworth_3 = 4,   
    sworth_4 = 5,   
    sworth_5 = 6,   
    sworth_6 = 7
    )
```



## Transform our data from wide format to long format

Hint: use the pivot longer function. Note, we want there to be one column for time, and another column for self worth scores. There should be one row of data for each individual's self worth score at each timepoint (e.g., 6 rows per person if no missing data)
```{r}
d_long <- d %>% 
  tidyr :: pivot_longer(., cols = contains("sworth"), names_to = "time", values_to = "self_worth")
```

## Clean your data in the long format 

Drop the rows with NA values (the observations where there is no self worth score). Then recode time into a factor and rename the variables so that they are only numeric values that range from time 0 to time 5 instead of time 1 to time 6. 
```{r}
d_long <- d_long %>%
  dplyr::filter(!is.na(self_worth)) %>%
  dplyr::mutate(time = as.numeric(gsub('sworth_', '', time)))
```

## Plot the data 

```{r}
ggplot(data = d_long, aes(x = time, y = self_worth)) +
  geom_line(aes(group = id), color = "black") + geom_point(aes(group = id), size = .4) +
  facet_wrap("id") + # Group variable gives text labels to panels
  labs(x = "Time",
       y = "Self-worth (self-reported)",
       title = "")
```

# 1. Model changes in self worth quality over time with random intercepts only

## Make a random intercept model with rstanarm


`(1|id)` gives us random intercepts - intercepts for each person

```{r}
m1 = rstanarm::stan_glmer(data = d_long, self_worth ~ time + (1|id), cores = 4)
```

## Interpret the 'fixed' effects: what is the average change in self worth over time?

```{r}
# point estimates
fixef(m1)

# summary of full model
m1_summary = summary(m1, probs = c(.025, .975))
m1_summary[1:10, ]
```

> -  Woah! Why so much output here? Well, all of the `b[(Intercept) id:1]` parameters in the model are the intercepts for each different person (each id), so there are a lot of parameterss in this model!


## Plot the model predictions for the group average self worth quality (fixed effect) over time

Note that we have to add the `re_formula = NA` argument in the `add_fitted_draws()` function here to tell the model we want just the fixed effects (not predictions for specific participants)
```{r}
# define a data frame of timepoints to make predictions for for the group average
pred_frame = data.frame(time = 1:6)

# get predictions using tidybayes! 
# again this dataframe looks a little confuisng at first, but the outcome (self worth) is always in the .value column
fixed_preds_m1 = tidybayes::add_fitted_draws(m1, newdata = pred_frame, re_formula = NA)

```

## Make a plot of predicted self worth over time
```{r}
ggplot(fixed_preds_m1, aes(x = time, y= .value)) +
  tidybayes::stat_lineribbon(alpha = .4) +
  labs(x = 'Time', y = 'Self-reported self worth')
```

a little nicer looking...
```{r}
ggplot(fixed_preds_m1, aes(x = time, y= .value)) +
  tidybayes::stat_lineribbon(alpha = .4) +
  labs(x = 'Time', y = 'Self-reported self worth') +
  scale_fill_brewer()
```


## Now, add 'spaghetti plot' lines for each participant's predicted communicaaion quality over time to your plot


We can get predictions for each point in the original dataset  by specificying the original dataset `d_long` as `newdata`

We also use the `median_qi` function to summarize and get the median and a posterior interval for each prediction. Now the `.value` column will be the median model prediction for each datapoint. We won't use `.upper` and `.lower` yet but we could later
```{r}
participant_preds_m1 = tidybayes::add_fitted_draws(model = m1, newdata = d_long) %>%
  median_qi(.value, .width = .95)
```

Now, add participant 'spaghetti' to our plot. We can see all the lines are parallel -- they have different intercepts (random intercepts!) but the same exact slope as the fixed effect. This is because there are no random slopes -- slopes cannot vary for different participants under this model.

```{r}
ggplot(fixed_preds_m1, aes(x = time, y= .value)) +
  geom_line(data = participant_preds_m1, aes(group = id), alpha = .5) +
  tidybayes::stat_lineribbon(alpha = .4) +
  labs(x = 'Time', y = 'Self-reported self worth') +
  scale_fill_brewer()
```

# 2. Model changes in self worth quality over time with random intercepts AND random slopes

## Make a random slope + intercept model with rstanarm

```{r}
m2 = rstanarm::stan_glmer(data = d_long, self_worth ~ time + (time|id), cores = 4)
```

## Interpret the 'fixed' effects: what is the average change in self worth over time?

```{r}
# point estimates
fixef(m1)

# summary of full model (a good way to shorten output)
g = summary(m2, probs = c(.025, .975))
g[1:5,]
```

Now, we have parameters like `b[time id:2]` for each participant! These are the random slopes!


## Plot the model predictions for the group average self worth quality (fixed effect) over time

Just like with the previous model

```{r}
fixed_preds_m2 = tidybayes::add_fitted_draws(m2, newdata = pred_frame, re_formula = NA)
```

```{r}
ggplot(fixed_preds_m2, aes(x = time, y= .value)) +
  tidybayes::stat_lineribbon(alpha = .4) +
  labs(x = 'Time', y = 'Self-reported self worth') +
  scale_fill_brewer()
```

## Now, add 'spaghetti plot' lines for each participant's predicted self worth quality over time to your plot

```{r}
participant_preds_m2 = tidybayes::add_fitted_draws(model = m2, newdata = d_long) %>%
  median_qi(.value, .width = .95)

ggplot(fixed_preds_m2, aes(x = time, y= .value)) +
  geom_line(data = participant_preds_m2, aes(group = id), alpha = .5) +
  tidybayes::stat_lineribbon(alpha = .4) +
  labs(x = 'Time', y = 'Self-reported self worth') +
  scale_fill_brewer()
```

> - Now, we can see the lines for the model predictions for different participants are not parallel! Random intercepts AND random slopes


## Create indiviudal panel plots to check our model fits

> - Now, we might want to make a panel plot to see how well our model fits each participant, comparing the model predicitons with the raw data

> - First we plot the raw data
> - Then we add in our predictions from m2 & m1

```{r,warning=FALSE, results='hide', message=FALSE}
panel_plt = ggplot(data = d_long, aes(x = time, y = self_worth)) +
  geom_line(aes(group = id, color = "Raw Daata")) + geom_point(aes(group = id), size = .4) + 
  geom_line(data = participant_preds_m1 , aes(x = time, y = .value, color = 'Random Intercept Model Prediction')) + 
  geom_line(data = participant_preds_m2, aes(x = time, y = .value, color = 'Random Slope Model Prediction')) + 
  facet_wrap("id") + # Group variable gives text labels to panels
  labs(x = "Time",
       y = "Self-reported self worth",
       title = "") 


# save to a file since it's a little hard to see inline
ggsave(panel_plt, file = 'panel_plt.png', height = 8, width = 15)
```

It's a little hard to see because the panels are small, but we can look to visually inspect which model seems to fit better!






