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

## Overview

The questions and exercises below all include examples using a dataset put together in the Tottenham Lab of many Billboard Hot 100 songs from 1946-1983. Today in this challenge you'll want to work to answer several questions:

1. Did pop and rock songs on the Billboard charts tend to get *louder* on averaage from 1946-1983?
2. From the 1946-1983 Billboard charts, *which were louder on average -- pop or rock songs?*
3. Did changes in average loudness from 1946-1983 differ for pop and rock songs?


* To answer these questions, please use a combination of Bayesian regression models with `rstanarm::stan_glm()`, and plotting with `ggplot()`
* If anything is confusing or you get stuck, please use the previous lesson plans and challenge answer keys as a reference, or reach out to an instructor!
* Remember, **we do not expect you to memorize syntax, nor is this a crucial part of coding or stats**. Use google, stack overflow, documentation, previous lessons and challenges, colleagues, and instructors as a reference -- this is what 'real-life' coding is like!

## 0. Setup 

**You don't have to write code in this section, but it will be helpful to read through to understand the data a little better.** 

Load the needed packages (install if necessary and *ASK AN INSTRUCTOR* if you're getting stuck here

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

Now, we can load in the data from Github using `read_csv()`

```{r}
music = read_csv('../../../DANLAB/Projects/AMFM/music_databases/familiar_music_database.csv')
names(music)
```
We have a lot of columns here we could look at, but for today's challenge just focus on genre, loudnesss, and year. We'll also only look at the `Pop` and `Rock` genres for today

```{r}
# select just the relevant columns
music = dplyr::select(music, title, artist, year, genre, loudness) %>%
  # filter the data to only include Pop/Rock genre
  dplyr::filter(genre %in% c('Pop', 'Rock')) %>%
  # scale/center loudness
  dplyr::mutate(loudness = as.vector(scale(loudness, center = TRUE, scale = TRUE)))
```

```{r}
head(music)
```


## 1. Did pop and rock songs on the Billboard charts tend to get *louder* on averaage from 1946-1983?

* Build a Bayesian regression model using `rstanarm::stan_glm()` with 1 continuous predictor to understand relationships between`year` and `loudness` in the music dataset!
* Interpret your model parameters: *what can you conclude about average changes in loudness between 1946-1983? how strong is your evidence?*
* Graph your model predictions using `tidybayes` and `ggplot()`
  * Bonus: graph model predictions and the raw data on the same plot!


Fit the model
```{r, warning=FALSE, results='hide', message=FALSE}
year_loud_model = rstanarm::stan_glm(data = music, loudness ~ year)
```


```{r}
# use probs = c(.025, .975) to get the middle 95% of the posterior
# use digits=3 to see more digits after the decimal for each term
year_loud_summary = summary(year_loud_model, probs = c(.025, .975), digits = 3)
year_loud_summary
```

>- On average, the music was get louder in later years. On average, each year was associated with an average  loudness increase of `r year_loud_summary[2,1]`, 95% PI [`r year_loud_summary[2,4]`, `r year_loud_summary[2,5]`] (in standardized units)


```{r}
# Make a dataframe of years, then get model predictions for average loudness each year in the dataset
year_df = data.frame(year = unique(music$year))

# get model predictions for each year
year_loudness_predictions = tidybayes::add_fitted_draws(model = year_loud_model, newdata = year_df)


head(year_loudness_predictions)
```

What's going on in this dataframe of predictions?

* Tidybayes has pulled out 4000 *estimates of predicted loudness* for each year -- these estimates are in the `.value` column. In general with this `add_fitted_draws()` function, the `.value` column is the outcome variable in your regression model
* We don't have to pay much attention to the `.row`, `.chain`, and `.iteration` columns. The `.draw` column is helpful though, because it helps us see that there are exactly 4000 'draws' (estimate 'drawn' from the posterior distribution) for each value of `year`. This will be true in general - tidybayes will generate 4000 draws for every combination of predictors in the dataframe given as the `newdata` argument


**Now, let's plot the predictions, using `tidybayes::stat_lineribbon()`, which gives us a nice 'ribbon band' that summarizes the 4000 predictions at different value for year and makes a smooth 'ribbon' for different levels of posterior uncertainty. Here, the darker inner band shows the 80% posterior interval and the outer lighter band shows the 95% posterior interval

```{r}
ggplot(data = year_loudness_predictions, aes(x = year, y = .value)) +
  tidybayes::stat_lineribbon(.width = c(.8, .95)) +
  scale_fill_brewer() +
  labs(fill = 'Posterior Interval Level', x = 'Year', y = 'Song Loudness (standardized)') +
  theme_bw()
```

Now, with points for raw data too. This helps us see how large any changes in average volume are on the scale of the distribution of all songs

```{r}
ggplot(data = year_loudness_predictions, aes(x = year, y = .value)) +
  geom_point(data = music, aes(x = year, y = loudness), alpha = .4) + 
  tidybayes::stat_lineribbon(.width = c(.8, .95)) +
  scale_fill_brewer() +
  labs(fill = 'Posterior Interval Level', x = 'Year', y = 'Song Loudness (standardized)') +
  theme_bw()

```



## 2. From the 1946-1983 Billboard charts, *which were louder on average -- pop or rock songs?*

* Build a Bayesian regression model using `rstanarm::stan_glm()` with 1 categorical predictor to understand relationships between `genre` and `loudness` in the music dataset!
* Interpret your model parameters: *what can you conclude about average difference in loudness between rock versus pop songs? how strong is your evidence?*
* Graph your model predictions using `tidybayes` and `ggplot()`
  * Bonus: graph model predictions and the raw data on the same plot!
  

```{r, warning=FALSE, results='hide', message=FALSE}
genre_loud_model = rstanarm::stan_glm(data = music, loudness ~ genre)
```
  
  
```{r}
# use probs = c(.025, .975) to get the middle 95% of the posterior
# use digits=3 to see more digits after the decimal for each term
genre_loud_summary = summary(genre_loud_model, probs = c(.025, .975), digits = 3)
genre_loud_summary
```
  
  
>- On average, the rock songs were louder than pop. This is indicated by the `genreRock` term which describes differencess in loudness for rock songs relative to pop songs. On average, each rock was associated with an average loudness increase of `r genre_loud_summary[2,1]`, 95% PI [`r genre_loud_summary[2,4]`, `r genre_loud_summary[2,5]`] (in standardized units) compared to pop music


```{r}
# Make a dataframe of years, then get model predictions for average loudness each genre in the dataset
# We use stringsAsFactors=FALSE because the genre variable was treated as a string, not as a factor, in the raw data fed into the model
genre_df = data.frame(genre = c('Rock', 'Pop'), stringsAsFactors = FALSE)

# get model predictions for each genre
genre_loudness_predictions = tidybayes::add_fitted_draws(model = genre_loud_model, newdata = genre_df)


head(genre_loudness_predictions)
```

Like with the first model, we get 4000 posterior predictions for each level of the predictors in `newdata`. Here, the dataframe is 8000 rows -- 4000 predictions in the `.value` column for the averaage loudness of *Rock* and *Pop* songs, respectively. 
  
  
  
One way of plotting this -- we use `tidybayes::stat_halfeye()` to get plots of the full posterior distributions plus points (posterior mean) and eror bars (we use the `.width` parameter to set here that we want an 80% (thicker bar) and 95% (thinner bar) posterior interval)

```{r} 
ggplot(data = genre_loudness_predictions, aes(x = .value, fill = genre)) +
  tidybayes::stat_halfeye(.width = c(.8, .95)) +
  labs(y = 'Relative Posterior Density', 
       x = 'Estimated Loudness (standardized)', 
       fill = 'Song Genre') +
  theme_bw()
```

Or we could flip the axes the other way

```{r} 
ggplot(data = genre_loudness_predictions, aes(x = genre, y = .value, fill = genre)) +
  tidybayes::stat_halfeye(.width = c(.8, .95)) +
  labs(y = 'Estimated Loudness (standardized)', 
       x = 'Song Genre', 
       fill = 'Song Genre') +
  theme_bw() +
  theme(legend.position = 'none')
```

We can also add the raw data in as jittered points, or a boxplot.

```{r} 
ggplot(data = genre_loudness_predictions, aes(x = genre, y = .value, fill = genre)) +
  geom_jitter(data = music, aes(x = genre, y = loudness), width = .4, alpha = .2) +
  tidybayes::stat_halfeye(.width = c(.8, .95)) +
  labs(y = 'Estimated Loudness (standardized)', 
       x = 'Song Genre', 
       fill = 'Song Genre') +
  theme_bw() +
  theme(legend.position = 'none')
```

Hmm, looks kinda weird - why?

  * The average differences between pop/rock song volume are pretty small, and we're also pretty sure about them, so the posterior distributions now look very narrow and spiky once we add in all the raw datapoints (check out how the y axis has changed from the previous plot)
  * The 80% and 95% posterior intervals for the estimates of average loudness for pop/rock songs are small enough that we can't even really see them on this plot -- they blend in with the point estimates. 
  
With a boxplot instead
```{r} 
ggplot(data = genre_loudness_predictions, aes(x = genre, y = .value, fill = genre)) +
  geom_boxplot(data = music, aes(x = genre, y = loudness), width = .3, alpha = .2, position = position_nudge(-.3)) +
  tidybayes::stat_halfeye(.width = c(.8, .95)) +
  labs(y = 'Estimated Loudness (standardized)', 
       x = 'Song Genre', 
       fill = 'Song Genre') +
  theme_bw() +
  theme(legend.position = 'none')
```
  
## 3. Did changes in average loudness from 1946-1983 differ for pop and rock songs?

* Build a Bayesian regression model using `rstanarm::stan_glm()` with 1 categorical predictor, 1 continuous predictor, and an interaction term between them, to understand `loudness` as a function of `genre`, `year`, and their interaction.
* Interpret your model parameters: *how do changes in loudness over time from 1946-1983 differ for pop versus rock songs? how strong is your evidence?*
  * And -- how do the results you are seeing here add more nuance to your interpretations from parts 1 & 2. 
* Graph your model predictions using `tidybayes` and `ggplot()`
  * Bonus: graph model predictions and the raw data on the same plot!
  
  
```{r}
# here we center year first to make the interaction term more interpretable
music = dplyr::mutate(music, year_center = as.vector(scale(year, center = TRUE, scale = FALSE)))

# the mean year in the dataset, so year_scale=0 represents this year
mean(music$year)
```


```{r, warning=FALSE, results='hide', message=FALSE}
interaction_model = rstanarm::stan_glm(data = music, loudness ~ year_center*genre)
```


```{r}
# use probs = c(.025, .975) to get the middle 95% of the posterior
# use digits=3 to see more digits after the decimal for each term
interaction_summary = summary(interaction_model, probs = c(.025, .975), digits = 3)
interaction_summary
```

From this model, the `year_center` term now describes the estimated difference in loudness associated with 1 year specifically for pop songs. The `genreRock` term describes the estimated difference in rock > pop music loudness at year = 0 (1967, since we centered year), and the `year_center:genreRock` term describes how *different* the associations between year and loudness are for rock songs, compared to pop songs. In other words, this interaction term describes how  the *slope* of year (where x is year, y = loudness) differs for rock > pop music. For each of these terms, the 95% posterior interval excludes 0, so we might interpret these as consistent associations in the dataset. 

Overall, the estimates indicate:
* that pop songs tended to be louder in later years as indicated by the `year_center` parameter
* At `year_center=0`, rock songs tended to be louder than pop songs
* rock songs had a *more negative* association between year and loudness compared to pop songs

However, this is a little hard to parse from the model output alone, so let's plot!

First, make a dataframe of years AND genre using expand.grid() -- now our dataframe is ALL combinations of genre and year present in the dataset
```{r}
interaction_df = expand.grid(genre = c('Rock', 'Pop'), 
                       year_center = unique(music$year_center),
                       stringsAsFactors = FALSE) 

# for plotting purposes we'll want the actual year back in
interaction_df$year = interaction_df$year_center + mean(music$year)

head(interaction_df)
```


```{r}
# get model predictions for each combindation of genre and year
interaction_predictions = tidybayes::add_fitted_draws(model = interaction_model, newdata = interaction_df)
head(interaction_predictions)
```

Great, now we have 4000 posterior predictions for loudness for every combination of `genre` and `year_center`, with the predictions stored in the `.value` column. We can make a plot now very similar to the first one we made with `tidybayes::stat_lineribbon()` only this time we'll set `fill=genre`

```{r}
ggplot(data = interaction_predictions, aes(x = year, y = .value, fill = genre)) +
  tidybayes::stat_lineribbon(.width = c(.95), alpha = .8) +
  labs(fill = 'Song Genre', x = 'Year', y = 'Song Loudness (standardized)') +
  theme_bw()
```

This gives us a nice clear visualization showing that while pop songs tended to be louder in later years, this wasn't true for rock songs. However, rock songs in the earlier years were louder than pop, with pop songs eventually catching up in loudness around 1975. Now let's add the raw data in too


```{r}
ggplot(data = interaction_predictions, aes(x = year, y = .value, fill = genre)) +
  geom_point(data = music, aes(x = year, y = loudness, color = genre)) +
  tidybayes::stat_lineribbon(.width = c(.95), alpha = .8) +
  labs(fill = 'Song Genre', x = 'Year', y = 'Song Loudness (standardized)') +
  theme_bw() +
  guides(color = FALSE)
```

Apart from contextualizing these trends on the scaale of the data, this graph also reveals something else important when we see the raw data -- there really weren't many rock songs in the dataset at all before the the late 1950s! So, how might this impact our conclusions if we hoped to say that pop music increased in loudness over time while rock did not? We might have to scale back our inference and dig a little deeper into what kinds of songs were released over time, and what processes lead to different songs being labeled as 'Rock' versus 'Pop' in the first place (i.e. record label and radio promotion strategies, geographical locationsand identities of the musical artists, who was doing the genre 'labeling'). 


## 4. Wrap-up

Congrats! You've finished the challenge! 

* Before you move on, we recommend that you take 5 minutes to go back and comment your code -- as a gift to yourself when you come back later to run similar analyses! Commenting and documenting right away is almost always going to ultimately save you time and energy.
* For any concepts or pieces of code that don't make sense yet, write down your questions! Be sure to either ask a friend or an instructor right away, or ask in slack so you can learn.


  
