```{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. 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 *among youth who went through gender transitions during the course of the study, how did communication quality with mothers change over time?* We'll use longitudinal multilevel models in `rstanarm` to approach this

# Load in packages and data

```{r}
library(tidyverse)
library(rstanarm)
library(bayesplot)
library(tidybayes)

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

# Setup (you don't need to modify this code)

## Filter data based on gender and include those with values of 4 or 5. 

Teens who identified as boys at T1 but then identified as girls or another gender at later timepoints were coded as 4, and teens who identified as girls at T1 but then identified as boys or another gender at a later timepoints were coded as 5. There will be a lot less adolescents we will be plotting here. 

```{r}
d <- d %>% 
  filter(GENDER_ACROSSTIME_T1T2T3T4 == 4 | GENDER_ACROSSTIME_T1T2T3T4 == 5)
```

## Select only the variables we will be plotting

Being: id, gender, and all of the communication with mothers scores
```{r}
d <- d %>% 
  select(ID, contains("MPACS"))
```


## Rename your variables 

Rename them so that they are easier to understand 
```{r}
d <- d %>% 
# rename vectors
  rename(
    id = 1,                   
    mom_comm_1 = 2,   
    mom_comm_2 = 3,   
    mom_comm_3 = 4,   
    mom_comm_4 = 5,   
    mom_comm_5 = 6,   
    mom_comm_6 = 7
    )
```

## Rescale communication

Right now it is on a scale of 20 to 100. Put it on a scale of 0 to 10
```{r}
d$mom_comm_1 <- 1/8 * (d$mom_comm_1 - 100) + 10
d$mom_comm_2 <- 1/8 * (d$mom_comm_2 - 100) + 10
d$mom_comm_3 <- 1/8 * (d$mom_comm_3 - 100) + 10
d$mom_comm_4 <- 1/8 * (d$mom_comm_4 - 100) + 10
d$mom_comm_5 <- 1/8 * (d$mom_comm_5 - 100) + 10
d$mom_comm_6 <- 1/8 * (d$mom_comm_6 - 100) + 10
```

## 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 communication scores. There should be one row of data for each individual's communication score at each timepoint (e.g., 6 rows per person if no missing data)
```{r}
d_long <- d %>% 
  tidyr :: pivot_longer(., cols = contains("mom_comm"), names_to = "time", values_to = "communication")
```

## Clean your data in the long format 

Drop the rows with NA values (the observations where there is no communication 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(communication)) %>%
  dplyr::mutate(time = as.numeric(gsub('mom_comm_', '', time)))
```

## Plot the data 

Note: some participants only have 1-2 timepoints of data here, but we can still include them in the model and not throw their data away! Multilevel models allow for different participants to have different amounts of data. Going forward, we might want to think more carefuly about *what it means* that these datapoints are missing though, and come up with a more nuanced strategy for dealing with  missing data. 

```{r}
ggplot(data = d_long, aes(x = time, y = communication)) +
  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 = "Communication Quality with Mothers",
       title = "")
```
# 1. Model changes in communication quality over time with random intercepts only

## Make a random intercept model with rstanarm

```{r}

```


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

Write a sentence or two describing what you found about changes in communication over time


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

```{r}

```

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

```{r}

```


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

## Make a random slope + intercept model with rstanarm

```{r}

```


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

Write a sentence or two describing what you found about changes in communication over time


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

```{r}

```


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

```{r}

```

## How does this spaghetti plot look different from the one with your model with only random intercepts? 


# 3. BONUS! Plot each person's raw data over time as a separate panel PLUS add the model-predicted changes over time for each participant on top. Try this for the random slopes model

Use the code below to get started with the plot of the raw data, but you'll have to add in the model-predicted lines:

```{r}
ggplot(data = d_long, aes(x = time, y = communication)) +
  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 = "Communication Quality with Mothers",
       title = "")
```








