---
title: "Top Bops and Pop Flops: What Makes a Hit Song?"
author: "Eric Lee"
date: 'Jun 2, 2024'
output:
  tufte::tufte_handout: default
  toc: true
  toc_depth: 2
  tufte::tufte_html: 
    citation_package: natbib
    latex_engine: xelatex
  tufte::tufte_book:
    citation_package: natbib
    latex_engine: xelatex
---

```{r setup, include=FALSE,echo=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(dplyr)
library(readr)
library(tidyr)
library(gridExtra)
library(tufte)
library(reshape2)
library(kableExtra)
# invalidate cache when the tufte version changes
knitr::opts_chunk$set(cache.extra = packageVersion('tufte'))
```

```{=tex}
\tableofcontents
\newpage
\renewcommand{\textsc}[1]{\smallcapsspacing{\textsmallcaps{#1}}}
```
# Abstract

Noting the evolving nature of the music industry in the modern era and the common desire among artists to create widely popular songs, this report investigates the factors that make a viral hit and how those factors have changed over time. To conduct this exploration, this report jointly uses data from the Billboard Hot 100 and audio features from the Spotify API to gather diverse features for all charting songs from 2000 to 2019, with two response variables: a song's peak chart ranking and total time spent on the charts. This data is preprocessed and cleaned to adjust for errors and outliers, and preliminary analysis is conducted on each feature with visualizations. This report uses forward selection to determine the significant features for the two response variables, and uses linear models and F-tests to verify the significance of the selected features. Further visualizations are plotted and analyzed both in aggregate and year-by-year for each such feature. Such analysis shows that the most important factors in determining chart performance pertain to the genre, energetic nature, production, and sensual impact of each song. The impact of these features shift over time with relatively steady trends, and impact peak performance and chart longevity in different ways. These results hold greater significance in enabling budding artists and the music industry to create targeted music more likely to become successful in the context of recent trends.

# Research Question and Significance

The music industry is a dynamic and ever-evolving landscape, characterized by the constant pursuit of creating hit songs that capture the attention of audiences worldwide. With the rise of digital streaming platforms and social media, the process of achieving viral success has become increasingly complex. Although countless songs are written and released by aspiring pop stars every day, only a select few are successful in ascending to the top of the charts.

In particular, we see a wide range of interesting observations. Although the top charts mainly feature well-established pop stars, there are scattered instances of indie artists making their first big break. Some artists enjoy remarkable consistency and longevity on the charts, while others only feature as one-hit wonders before fading into obscurity. While songs of different genres, backgrounds, and styles achieve virality, many of the most popular songs adhere to familiar structures and conventions. Moreover, while English songs dominate the charts, there are also appearances from songs in international languages like Korean pop and Latin bachata. These phenomena all lead to one central question: **exactly which factors make a song a viral hit, especially in the modern era?**

The answer to this question holds significance in many regards. First, the ability to release a chart-topping hit is critically important to the livelihoods of budding artists and established stars alike. Producing a song requires time and fixed costs, so maximizing the popularity of each song is necessary for artists to sustain themselves. Specifically, a viral hit can bring an artist not only high streaming numbers that translate directly to earnings, but also a level of exposure that brings brand deals and fame. Second, these insights are crucial in elevating the music of artists who do not fit the mold of the status quo, such as indie artists, international singers, and stylistic pioneers. By discovering the factors behind hit songs and applying that knowledge to bring exposure to unique songs and artists, we can bring innovation to the music industry while elevating underrepresented voices in a saturated field. Third, viral hits often transcend mere commercial success and shape popular culture, so understanding the recipe behind viral hits can yield insights about cultural dynamics, social discourse, and popular ideas. Especially since chart-topping hits are often continuously played in phones, cars, and stores all over the world, their influence is extremely pervasive. Moreover, on social media platforms, these hits govern the top trends with regard to language, dance, and humor.

# Context and Novelty

Given the significance of this topic, there have been many prior efforts in this field to predict and analyze the popularity of songs based on various audio features. These articles are valuable in informing a strong foundation and logical approach to exploring this question, but they leave an evident gap in the literature investigating the post-2000 era of streaming and social media with a comprehensive dataset of both the top charts and audio features.

An early effort in 2017 was Nijkamp's comprehensive analysis[^1] of the relationship between the audio features and stream counts of 1000 songs sourced from Spotify. In this analysis, Nijkamp focused on the impact of audio features in isolation, discovering that certain features were positively or negatively correlated to higher stream counts. While this feature-wise approach offer insights with regard to the relative effects of each audio feature, the analysis is limited with regard to evaluating these features in tandem. Furthermore, a further limitation lies in the observation that the feature of interest—stream counts—was taken from Spotify and thus excludes other important factors in popularity such as stream counts on other platforms, plays on social media, and plays on radio. In my project, I aim to expand on this article's feature-wise analysis with additional evaluation of multi-feature models and use data from the Billboard charts[^2] whose metric for song popularity factors in plays and streams from other platforms.

[^1]: [Sean Miller](https://www.kaggle.com/datasets/thedevastator/billboard-hot-100-audio-features?select=Hot+Stuff.csv) (2024) *Billboard Hot 100 Weekly Charts with Spotify Audio Features*, Kaggle

[^2]: [Eric Lee](https://drive.google.com/file/d/1GFjUR-kvfFS_BnVX0GUzVwL5YF-JrBNF/view?usp=sharing) *Data Processing.ipynb*

A related 2020 analysis[^3] expands on this research by investigating trends in the featured songs on the Billboard charts. In particular, the article finds that in recent decades, there is more homogeneity with regard to the number of unique artists represented on the top charts. Further, this analysis finds that recent songs are less likely to stay on the chart for a long time after debut. These insights are notable and point to the need for more focused explorafigtion into the popular hits of the last two decades to explain these changes. Further, while this analysis solely focused on insights to be gleaned from Billboard data alone, I look to delve deeper by combining this dataset with Spotify data in order to analyze features at the level of individual songs.

[^3]: [Eric Lee](https://drive.google.com/file/d/109z3NgYKWcM2CdBDgkIHL2F0XhPkvTPd/view?usp=sharing) *Processed Data.csv*

Many other research efforts investigating this question used the Million Song Database provided by Columbia University[^4], such as articles by Nasreldin[^5] and Pham et al[^6]. This dataset is a comprehensive dataset of one million songs published until 2011, containing information about audio features and metadata for each song. These articles corroborate in identifying genre and artist familiarity as prominent features in determining popularity, informing my hypotheses for my analysis. However, this dataset is inherently limited due to its scope in only including songs up to 2011 (missing more than a decade of music up to 2024) and its lack of information regarding performance on the charts. Therefore, I aim to fill this gap to glean insights about the last decade of music that these articles are unable to address with a more updated and comprehensive dataset which includes information such as weeks on the top charts.

[^4]: [Eric Lee](https://docs.google.com/spreadsheets/d/1bgyyLFO86t_X-8LI7TTKgepAEeYuyizLdKRhps-vtVs/edit?usp=sharing) *Data Dictionary*

[^5]: [Mohamed Nasreldin](https://towardsdatascience.com/song-popularity-predictor-1ef69735e380) (2018) *Song Popularity Predictor*, Medium

[^6]: [Pham et al.](https://cs229.stanford.edu/proj2015/140_report.pdf) (2015) *Predicting Song Popularity*, Stanford Department of Computer Science

# Hypotheses

Considering the many attributes of each viral song, there are many potential hypotheses to investigate.

First, it makes sense to explore whether certain audio features translate to more success and popularity on the charts. For example, it seems likely that songs with characteristics related to high energy and catchiness are more likely to be top hits due to their affinity to being played on the radio and in stores. Specifically, I hypothesize that songs with higher values for `danceability`, `energy`, `valence`, and `tempo` as defined by Spotify will be more likely to reach higher peaks on the Billboard charts and stay on the Billboard charts for a longer duration.

Second, I hope to investigate how the genre impacts popularity and longevity on the charts. While many viral songs fall into the general category of "pop", there are exceptions. Further, there are many subgenres within pop, and it is relevant to see whether certain styles make a song more conducive to become a hit (ex. Latin pop, reggaeton, etc.). In particular, I hypothesize that songs that are in the genre of `pop` have higher peaks and longer stays on the charts. Further, I hypothesize that songs in other recognized "non-`pop`" genres like rap and country have perform better on the charts than songs who are not in any commonly recognized genre (but not better than `pop`).

Third, it is also prudent to consider other factors such as the effects of collaborations with popular artists and the recognition gained from appearing on the charts to begin with. On this axis, comparisons will be made between songs with and without famous collaborators. Here, I hypothesize that songs including collaborations are more likely to reach success on the charts in comparison to songs by solo artists.

# Data Processing

My original dataset[^7] included two tables: one containing weekly data regarding song rankings on the Billboard Hot 100 charts, and the other containing data for audio features sourced from Spotify. To investigate my hypotheses regarding the influence of audio features, song genre, and collaboration on chart performance and longevity, I combined these two tables into a merged table and created additional variables. For purposes of reproducibility, see the references in the margin to access the Python code[^8] used for the data processing and the resulting dataset[^9] (CSV file).

[^7]: [Sean Miller](https://www.kaggle.com/datasets/thedevastator/billboard-hot-100-audio-features?select=Hot+Stuff.csv) (2024) *Billboard Hot 100 Weekly Charts with Spotify Audio Features*, Kaggle

[^8]: [Eric Lee](https://drive.google.com/file/d/1GFjUR-kvfFS_BnVX0GUzVwL5YF-JrBNF/view?usp=sharing) *Data Processing.ipynb*

[^9]: [Eric Lee](https://drive.google.com/file/d/109z3NgYKWcM2CdBDgkIHL2F0XhPkvTPd/view?usp=sharing) *Processed Data.csv*

## Merging Tables

In particular, I started by creating new columns `Year`, `Month`, and `Day` from the original column `WeekID` to interpret date data as integers and at a greater level of granularity. Then, I used the new `Year` field to filter out all rows not between 2000 and 2019 (the desired time period). Afterward, I combined duplicate rows in the chart dataset by keeping the entry with the highest `PeakPosition` and counting the total number of `WeeksOnChart`, with the correct date values to match.

Then, I combined the updated `PeakPosition` and `WeeksOnChart` tables, only keeping rows where the unique key `SongID` is present in both tables. After renaming columns to be more readable, the resulting dataset contains each unique song that appeared at least once on the Hot 100 charts between 2000 and 2019, with information about its audio features, peak chart performance, the date of that peak, and the total number of weeks on the chart (in the 2000-2019 timeframe).

## Creating New Variables

Next, I extracted genre data from the `GenreList` variable by creating new columns for the most popular genres to define broader genre categories: `Pop`, `Country`, `Rock`, `R&B`, `RapHipHop`, and `Indie`. For each of these genres, I looked for the presence of certain target substrings within the `GenreList`, adding value 1 if the substring is found and 0 otherwise. For example, if a song's `GenreList` had the substring `pop` or `dance`, the new `Pop` column would have value 1. Additionally, I added two more columns related to genre: `OtherGenre` as a binary indicator for when a song fails to fit in any of the above genres, and `NoGenre` as a binary indicator for when a song had an empty `GenreList` from Spotify.

Finally, I created new columns to investigate variables related to artist collaboration from the `Performer` column. In the original dataset, all featured collaborators in the `Performer` variable are listed after the string 'Featuring', with additional '&' and ',' added when there are three or more featured performers. Using splitting functions for these delimiter strings, I created a list of featured performers for each song not including the original performer(s) in a new column `CollaboratorsList`. An additional column that was created is `Collaborators`, counting the number of collaborators in `CollaboratorsList`. Further, I created the additional column `FamousCollaborator` as an indicator to whether any of the collaborators in the song's `CollaboratorsList` had a prior hit on the charts, which was processed in Python via for loops based on prior `Performers` sorted by `Month`, `Year`, and `Day`.

# Final Dataset and Data Dictionary

As a result of the data processing steps described above, I arrived at the final dataset which can be loaded into R with the code below (with some renaming of columns for readability). Note that the relevant CSVs are attached in the margins above for reproducibility. This final dataset has 8,664 rows (songs) and 34 columns.

```{r, echo=FALSE}
path <- "/Users/ericlee/Desktop/Spring 2024/DATASCI 120/"
data <- read.csv(paste0(path, "DATASCI 120 - Processed Data.csv"))
names(data)[names(data) == 'Peak.Position'] <- 'PeakPosition'
names(data)[names(data) == 'Weeks.on.Chart'] <- 'WeeksOnChart'
names(data)[names(data) == 'Genre.List'] <- 'GenreList'
names(data)[names(data) == 'Duration..ms.'] <- 'Duration'
names(data)[names(data) == 'Time.Signature'] <- 'TimeSignature'
names(data)[names(data) == 'R.B'] <- 'R&B'
names(data)[names(data) == 'Rap.Hip.Hop'] <- 'RapHipHop'
names(data)[names(data) == 'Other.Genre'] <- 'OtherGenre'
names(data)[names(data) == 'No.Genre'] <- 'NoGenre'
names(data)[names(data) == 'Collaborators.List'] <- 'CollaboratorsList'
names(data)[names(data) == 'Famous.Collaborator'] <- 'FamousCollaborator'
data <- data %>%
  filter(Year != 2020)
```

For more detailed information regarding the variables included in the dataset and the corresponding descriptions, types, ranges, and reproducibility notes, please refer to the data dictionary[^10] I created. After running `summary(data)`, I see that there are no obvious mistakes or nonsense values in the data, as all the ranges make sense. More detailed analysis of the distribution of the data with regard to each feature follows below. Further, the information about variable ranges gleaned from this command have been updated into the data dictionary.

[^10]: [Eric Lee](https://docs.google.com/spreadsheets/d/1bgyyLFO86t_X-8LI7TTKgepAEeYuyizLdKRhps-vtVs/edit?usp=sharing) *Data Dictionary*

## Quality Control

With regard to quality control of this dataset, an inspection of the table shows no obvious nonsense values. Additionally, the reputability of the sources (Billboard and Spotify) gives confidence that the chance for errors in the data is relatively low. Even in the case where there were inherent errors in the original data tables, the data processing described above addresses many of these concerns. Namely, since the merge between tables was an inner join on the unique key `SongID`, songs with incorrectly spelled IDs or mismatches between tables were excluded from the final dataset. Further, given that all rows representing songs outside of the years 2000-2020 was sliced out of the final dataset, there is a reduced chance of error caused by limitations in tracking or technology.

Although there is low likelihood for errors in the dataset, we note that some rows in the final table have `NA` values for certain audio features due to limitations in accessing the Spotify API. We observe that there are 7,849 rows of complete data (without `NA` values for any variable) out of the 8,664 rows in the final dataset, meaning there are 815 rows with at least one missing value.

Given that a very high percentage (90.5%) of the rows in the final dataset have values for each variable, all models in the analysis conducted below simply exclude the rows with missing values. Additionally, we see in Figure 1 that the 815 songs with at least one missing value are spread relatively evenly across each year without a notable trend, indicating that these missing values can attributed to randomness as opposed to a meaningful confounding variable. Testing label as in \ref{testing figure}.

```{r, echo=FALSE, cache=TRUE, fig.margin=TRUE, fig.cap="Distribution of songs with one or more N/A values per year.\\label{testing figure}"}
data_with_na <- data[!complete.cases(data), ]
ggplot(data_with_na, aes(x = Year)) + 
  geom_bar(fill = 'blue') + 
  labs(x = "Year", y = "Frequency of Songs with 1+ NA Value", title = "Number of Songs with 1+ NA Value by Year") +
  scale_x_continuous(limits = c(1999, 2020), breaks = seq(2000, 2019, 1)) +
 theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1))
```

Thus, the full linear model conducted on all variables utilizes 90.5% of the original dataset, with all other models focusing on a subset of the variables utilizing at least 7,849 rows. Since analyses will be conducted with a shifting focus on each variable, we leave `NA` values in this main `data` table while removing the relevant rows in individual code snippets by variable below.

# Summary of Variables

Next, we take a deeper dive into the data with exploratory data analysis. In this analysis, chart performance (`PeakPerformance`) and chart longevity (`WeeksOnChart`) are considered the outcome variables while the other numeric/boolean variables are considered the explanatory variables.

## Outcome Variables

We use two variables to quantify virality and success on the charts, namely peak position and longevity. Specifically, peak position refers to the highest rank a given song achieves throughout its stay on the Hot 100 charts. First, we investigate the distribution of the peak positions of songs across this time period in Figure 2. As expected, we see that there are many more songs that peak near the bottom of the charts, and increasingly fewer that reach the very top of the charts near the top position. However, it is interesting to note that there is a slight spike around the peak position of 75, indicating that it may be easier for songs to reach that area of the charts before plateauing.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin = TRUE, fig.cap="Distribution of songs by peak position."}
ggplot(data, aes(x = PeakPosition)) +
  geom_histogram(binwidth = 5, fill = "skyblue", color = "black") +
  labs(title = "Distribution of Songs by Peak Position", x = "Peak Position", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_x_reverse()
```

We also investigate the distribution of songs by longevity, namely in weeks remaining on the charts. The longevity of a song is defined by the cumulative sum of the weeks it spends on the Hot 100 before being removed. The general trend of the histogram in Figure 3 is unsurprising, that the majority of songs fall off the charts very quickly after debut, with fewer and fewer remaining for a very long time (50+ weeks) as seen by the heavy right tail. Although, there is a very noticeable spike at 20 weeks, indicating that many songs in the database fell off of the charts after exactly 20 weeks. This unusual behavior can be attributed to Billboard's chart policy that songs which are "descending are removed from the Billboard Hot 100...after 20 weeks and if ranking below No. 50" as a way to combat chart inertia and feature new songs gaining popularity and momentum[^11]. This policy does not necessitate any major changes in our analysis, as it represents the phenomenon of "slipping off the charts" due to obscurity. However, we should remain cognizant that the 20 week mark is likely to appear very frequently, and represents a song that could have potentially lasted slightly longer than 20 weeks (before officially falling off).

[^11]: [Billboard](https://www.billboard.com/billboard-charts-legend/) (2024) *Billboard Charts Legend: Recurrent Rules*, Billboard

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin = TRUE, fig.cap="Distribution of songs by total number of weeks on chart."}
ggplot(data, aes(x = WeeksOnChart)) +
  geom_histogram(binwidth = 3, fill = "tan", color = "black") +
  labs(title = "Distribution of Songs by Longevity", x = "Weeks on Chart", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 
```

## Audio Features

Next, with a focus on gaining insight into our first hypothesis regarding the importance of different audio features, we explore the distributions of each audio feature alongside its relationship with the outcome variables.

Before delving into individual audio features, as an exploratory overview we examine the boxplots for each audio feature with relation to peak chart performance and chart longevity. First, we consider the audio features generated by Spotify as scores between 0 and 1. This display gives a preliminary look at the general center and spread of each of these variables. Notably, `instrumentalness` seems to be at 0 for the vast majority of songs, while `mode` seems to be at 1 for most songs. Additionally, the centers for `acousticness`, `liveness`, and `speechiness` are relatively low, indicating that these features have a lower general score for recent songs that reach the charts. No outliers are noticed here, and all data falls within 0 and 1 (inclusive), as expected.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width = 7, fig.height = 5, fig.cap="Boxplots for audio features with scores between 0 and 1."}
columns_of_interest <- c('Danceability', 'Energy', 'Mode', 'Speechiness', 'Acousticness', 'Instrumentalness', 'Liveness', 'Valence')
data_long <- pivot_longer(na.omit(data[, columns_of_interest]), all_of(columns_of_interest), names_to = "Variable", values_to = "Value")

audio_feat_boxplot <- ggplot(data_long, aes(x = Variable, y = Value)) +
  geom_boxplot(fill = "grey") +
  labs(title = "Boxplots for Audio Features Between 0 and 1", x = "Audio Feature", y = "Value") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
    theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

audio_feat_boxplot
```

## Outliers

Similarly, we can see a summary of the center and spread of the other audio features (not between 0 and 1) by examining the boxplots below. Apart from the observations of the different centers of each features, there are a few notable outliers that need to be addressed.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.cap="Boxplots for audio features duration, key, and loudness."}
library(patchwork)

boxplot_duration <- ggplot(data[!is.na(data$Duration), ], aes(x = "", y = Duration)) +
  geom_boxplot(fill = "grey") +  
  labs(title = "Duration", x = "Duration", y = "Value") +  
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

boxplot_key <- ggplot(data[!is.na(data$Key), ], aes(x = "", y = Key)) +
  geom_boxplot(fill = "grey") +  
  labs(title = "Key", x = "Key", y = "Value") +  
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

boxplot_loudness <- ggplot(data[!is.na(data$Loudness), ], aes(x = "", y = Loudness)) +
  geom_boxplot(fill = "grey") +  
  labs(title = "Loudness", x = "Loudness", y = "Value") +  
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

combined_boxplots_1 <- boxplot_duration + boxplot_key + boxplot_loudness 

print(combined_boxplots_1)

```

```{r, message = FALSE, warning=FALSE, echo=FALSE, fig.cap="Boxplots for audio features tempo and time signature."}
library(patchwork)

boxplot_tempo <- ggplot(data[!is.na(data$Tempo), ], aes(x = "", y = Tempo)) +
  geom_boxplot(fill = "grey") +  
  labs(title = "Tempo", x = "Tempo", y = "Value") +  
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

boxplot_ts <- ggplot(data[!is.na(data$TimeSignature), ], aes(x = "", y = TimeSignature)) +
  geom_boxplot(fill = "grey") +  
  labs(title = "Time Signature", x = "Time Signature", y = "Value") +  
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

combined_boxplots_2 <- boxplot_tempo + boxplot_ts

print(combined_boxplots_2)
```

First, as seen in the Duration boxplot, there is a song that is extremely long. Using code, this outlier is R. Kelly's song, "Trapped in the Closet". Checking with the original Spotify track[^12] shows that this song is indeed 16 minutes (992160 ms), so this outlier is not born from error.

[^12]: [R. Kelly](https://open.spotify.com/track/45LrQ3tg1z8plpuQRCuSwE?si=ff2b73f4e7c548ba) (2005) *Trapped In The Closet*, Spotify

```{r, echo=FALSE}
# Investigate outlier of R. Kelly's Trapped in the Closet with 992160ms (16 min)
d1 <- data[!is.na(data$Duration), ]
duration_outlier <- d1[d1$Duration > 900000, ]
# print(duration_outlier)
```

Second, there is a very notable outlier near the bottom of the Loudness boxplot. This outlier is Billie Eilish's song, "Listen Before I Go" with -23.023 dB. This value is unverifiable with other sources, so the value is reverted to `NA` out of an abundance of caution.

```{r, echo=FALSE}
# Investigate outlier of Billie Eilish's Listen Before I Go with -23.023 loudness
d2 <- data[!is.na(data$Loudness), ]
loudness_outlier <- d2[d2$Loudness < -20, ]
# print(loudness_outlier)
data$Loudness[data$SongID == "Listen Before I GoBillie Eilish"] <- NA
```

Third, there is a seemingly nonsense value in the time signature boxplot showing a song with time signature 0 (which is not possible). This point turns out to be "Imma Be" by "The Black Eyed Peas", which indeed has 0 as the value for `TimeSignature` in the data table. This is clearly an error, as the song has a time signature of 4 beats per bar [^13]. Therefore, we change the value for this song from 0 to 4 for `TimeSignature`.

[^13]: [SongBPM](https://songbpm.com/the-black-eyed-peas/imma-be) *Song Metrics*, Imma Be

```{r, echo=FALSE}
# Investigate outlier of Black Eyed Peas' Imma Be with 0 time signature
d3 <- data[!is.na(data$TimeSignature), ]
ts_outlier <- d3[d3$TimeSignature == 0, ]
# print(ts_outlier)
data$TimeSignature[data$SongID == "Imma BeThe Black Eyed Peas"] <- 4
```

In addition to considering all the outliers, we notice that the distribution for all variables is approximately normal, with the exception of Instrumentalness, Acousticness, Liveness, and Speechiness (which are skewed with a heavy tail). Thus, we apply a log transformation to these four variables such that the resulting distribution is much more normal. With these transformed features being approximately normal, we gain confidence in their application within a linear model. The figures showing the distributions for all of these audio features can be found in Figures 38-39 in the Appendix.

## Genres

In the final dataset, a percentage of songs were not labeled with any data, and these songs are classified as `NoGenre`. For the other songs which have at least one labeled genre, they are labeled as `Known Genre` if they fit into `Pop`, `Country`, `Rock`, `Blues`, `R&B`, `RapHipHop`, or `Indie`, and `Other Genre` if they fail to fit into any of the above categories. The bulk of the songs fit into at least one of the aforementioned known genres, and that hence these 7 genres encompass the vast majority of the songs that end up on the Hot 100 charts. The following table displays the distribution of songs across genres (noting that each song can belong to multiple genres). More detailed boxplots for each of these genre variables can be found in Figures 40-41 in the Appendix.

```{r, echo=FALSE, warning=FALSE, message=FALSE}
genre_df <- data.frame(
  Genre = c("Pop", "Country", "Rock", "R&B", "Indie", "Rap/Hip-Hop", "Blues", "OtherGenre", "NoGenre"),
  PercentageInGenre = c("65.6%", "18.1%", "20.5%", "13.0%", "2.2%", "38.9%", "0.3%", "7.7%", "4.6%"),
  PercentageNotInGenre = c("34.4%", "81.9%", "79.5%", "87.0%", "97.8%", "61.1%", "99.7%", "92.3%", "95.4%"))

colnames(genre_df) <- c("Genre", "% of Songs in Genre", "% of Songs Not in Genre")

knitr::kable(genre_df, align = 'c') %>%
  column_spec(1, width = "110px") %>%
  column_spec(2, width = "90px") %>%
  column_spec(3, width = "90px") 
```

```{r, echo=FALSE}
Pop_count <- sum(data$Pop == 1)
# print(Pop_count / nrow(data))
# 0.6558572
```

```{r, echo=FALSE}
Country_count <- sum(data$Country == 1)
# print(Country_count / nrow(data))
# 0.1813725
```

```{r, echo=FALSE}
Rock_count <- sum(data$Rock == 1)
# print(Rock_count / nrow(data))
# 0.2051282
```

```{r, echo=FALSE}
Blues_count <- sum(data$Blues == 1)
# print(Blues_count / nrow(data))
# 0.003393665
```

```{r, echo=FALSE}
RnB_count <- sum(data$`R&B` == 1)
# print(RnB_count / nrow(data))
# 0.1303419
```

```{r, echo=FALSE}
RHH_count <- sum(data$RapHipHop == 1)
# print(RHH_count / nrow(data))
# 0.388888
```

```{r, echo=FALSE}
Indie_count <- sum(data$Indie == 1)
# print(Indie_count / nrow(data))
# 0.0224987
```

```{r, echo=FALSE}
other_count <- sum(data$OtherGenre == 1)
# print(other_count / nrow(data))
# 0.07717446
```

```{r, echo=FALSE}
nogenre_count <- sum(data$NoGenre == 1)
# print(nogenre_count / nrow(data))
# 0.04587733
```

## Collaboration

Note that although a minority of songs (19.9%) feature famous collaborators, the songs that do have the ability to reach the top of the charts. It is relevant to note that the average peak position of all songs without famous collaborators is 75.3505, which is only slightly worse than the average peak position of all songs with famous collborators, which is 74.6282.

```{r, echo=FALSE}
zero <- mean(data$PeakPosition[data$FamousCollaborator == 0], na.rm = TRUE)
one <- mean(data$PeakPosition[data$FamousCollaborator == 1], na.rm = TRUE)

FamousCollab_count <- sum(data$FamousCollaborator == 1)
# print(FamousCollab_count / nrow(data))
#0.0.1998492
```

In addition to considering the effects of the presence of famous collaborators, it is also important to consider the total number of collaborators on a given song (who may or may not be famous or previously featured on the charts). We see that the vast majority of songs do not have collaborators, and for the ones who do, the number of collaborators is usually one and sometimes two. Instances of three or more collaborators are a rare occurrence on the charts. The distribution of these collaboration variables can be found in the Appendix in Figure 42.

```{r, echo=FALSE}
path <- "/Users/ericlee/Desktop/Spring 2024/DATASCI 120/"
data <- read.csv(paste0(path, "DATASCI 120 - Processed Data.csv"))
names(data)[names(data) == 'Peak.Position'] <- 'PeakPosition'
names(data)[names(data) == 'Weeks.on.Chart'] <- 'WeeksOnChart'
names(data)[names(data) == 'Genre.List'] <- 'GenreList'
names(data)[names(data) == 'Duration..ms.'] <- 'Duration'
names(data)[names(data) == 'Time.Signature'] <- 'TimeSignature'
names(data)[names(data) == 'R.B'] <- 'R&B'
names(data)[names(data) == 'Rap.Hip.Hop'] <- 'RapHipHop'
names(data)[names(data) == 'Other.Genre'] <- 'OtherGenre'
names(data)[names(data) == 'No.Genre'] <- 'NoGenre'
names(data)[names(data) == 'Collaborators.List'] <- 'CollaboratorsList'
names(data)[names(data) == 'Famous.Collaborator'] <- 'FamousCollaborator'
data <- data %>%
  filter(Year != 2020)
data$Speechiness <- log(data$Speechiness + 0.01)
data$Acousticness <- log(data$Acousticness + 0.01)
data$Instrumentalness <- log(data$Instrumentalness + 0.01)
data$Liveness <- log(data$Liveness + 0.01)
```

# Determining Significant Features

After seeing the initial distributions and relationships between variables with evidence supporting the current three hypotheses, we use linear models to further identify the most positively impactful variables on peak chart performance and chart longevity. We start by running a baseline linear model on all features in the dataset (only including rows without `NA` values) to extract insights about the full model. Then, we use forward selection for feature selection and conduct F-tests with reduced linear models only incorporating the selected features to analyze the significance of the results. Afterward, we conduct this same analysis for each subset of feature sorted by hypothesis (audio features, genre, and collaboration), by training a baseline linear model, conducting feature selection, and assessing reduced models with F-tests.

## Models and Feature Selection

We begin by training two baseline linear model on all features in the dataset, with one predicting `PeakPosition` and another predicting `WeeksOnChart` (only on rows with no `NA` values).

The full baseline linear model for `PeakPosition` is `full_base_peak_model`, whose summary is below. The and p-values for the coefficients of this naive baseline provide that variables `Energy`, `Liveness`, `OtherGenre`, `Pop`, `Country`, `Rock`, `R&B`, `RapHipHop`, `Indie`, and `Collaborators` are significant in predicting `PeakPosition`. Notably, many of these variables pertain to genre rather than audio features, giving evidence that genre is a dominant factor in peak performance.

Next, we train the full baseline linear model titled `full_base_weeks_model` for `WeeksOnChart`, whose summary is below. As seen by the p-values for the coefficients in the summary, the most impactful variables with regard to chart longevity seem to be `Duration`, `Danceability`, `Energy`, `Loudness`, `Speechiness`, `Acousticness`, `Liveness`, `Valence`, `Pop`, `Country`, `Rock`, and `R&B`. In contrast to the model trained to predict peak performance on the charts, this subset of variables features many more audio features and far less genre features. This observation supports the idea that while genre is important in determining peak performance on the charts, the audio features of the song have a notable impact on chart longevity.

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.cap="Coefficients for the full models for peak position and weeks on charts."}
full_model_df <- data.frame(
  PeakCoefficients = c("Intercept", "Duration", "Danceability", "Energy", "Key", "Loudness", "Mode", "Speechiness", "Acousticness", "Instrumentalness", "Liveness", "Valence", "Tempo", "TimeSignature", "OtherGenre", "NoGenre", "Pop", "Country", "Rock", "Blues", "R&B", "RapHipHop", "Indie", "FamousCollaborator", "Collaborators"),
  PeakEstimates = c("7.324e+001", "-8.319e+006", "-1.159e+000", "8.815e+000", "-1.027e-02", "9.040e-02", "6.193e-03", "2.201e-01", "3.362e-01", "8.457e-01", "-1.254e+00", "1.993e+00", "-2.982e-03", "-1.554e+00", "7.458e+00", "4.744e-01", "1.872e+00", "1.400e+01", "6.407e+00", "1.970e+00", "6.653e+00", "3.437e+00", "1.002e+01", "-3.109e-01", "2.094e+00"),
  PeakSE = c("6.311e+00", "6.134e-06", "2.345e+00", "2.765e+00", "7.693e-02", "1.833e-01", "6.124e-01", "4.842e-01", "2.490e-01", "5.276e-01", "4.615e-01", "1.538e+00", "9.628e-03", "1.030e+00", "1.773e+00", "2.490e+00", "7.677e-01", "9.677e-01", "7.928e-01", "4.777e+00", "8.616e-01", "7.920e-01", "1.870e+00", "1.177e+00", "7.489e-01"),
  WeeksEstimates = c("-8.383e+00", "2.110e-05", "8.046e+00", "-1.295e+00", "1.366e-02", "2.061e-01", "1.483e-01", "-1.260e+00", "-2.850e-01", "-5.013e-01", "-8.759e-01", "3.728e+00", "-5.402e-03", "8.280e-01", "-1.924e+00", "1.575e+00", "7.247e-01", "2.133e+00", "3.534e+00", "6.709e-01", "3.256e+00", "-2.800e-01", "-1.709e+00", "1.445e-02", "5.380e-02"),
  WeeksSE = c("2.906e+00", "2.825e-06", "1.080e+00", "1.273e+00", "3.542e-02", "8.443e-02", "2.820e-01", "2.230e-01", "1.146e-01", "2.429e-01", "2.125e-01", "7.084e-01", "4.433e-03", "4.743e-01", "8.165e-01", "1.147e+00", "3.535e-01", "4.456e-01", "3.651e-01", "2.200e+00", "3.968e-01", "3.647e-01", "8.610e-01", "5.422e-01", "3.448e-01")
)
colnames(full_model_df) <- c("Variables", "PeakPosition Coefficient", "PeakPosition Std. Error", "WeeksOnChart Coefficient", "WeeksOnChart Std. Error")

knitr::kable(full_model_df, align = 'c') %>%
  column_spec(1, width = "100px") %>%
  column_spec(2, width = "70px") %>%
  column_spec(3, width = "70px") %>%
  column_spec(4, width = "70px") %>%
  column_spec(5, width = "70px")
```

```{r, echo=FALSE}
# Remove any rows with NA in any column 
full_data <- data[complete.cases(data), ]
full_base_peak_model <- lm(PeakPosition ~ Duration + Danceability + Energy + Key + Loudness + Mode + Speechiness + Acousticness + Instrumentalness + Liveness + Valence + Tempo + TimeSignature + OtherGenre + NoGenre + Pop + Country + Rock + Blues + `R&B` + RapHipHop + Indie + FamousCollaborator + Collaborators, data = full_data)
# summary(full_base_peak_model)
```

```{r, echo=FALSE}
# Remove any rows with NA in any column 
full_data <- data[complete.cases(data), ]
full_base_weeks_model <- lm(WeeksOnChart ~ Duration + Danceability + Energy + Key + Loudness + Mode + Speechiness + Acousticness + Instrumentalness + Liveness + Valence + Tempo + TimeSignature + OtherGenre + NoGenre + Pop + Country + Rock + Blues + `R&B` + RapHipHop + Indie + FamousCollaborator + Collaborators, data = full_data)
# summary(full_base_weeks_model)
```

Next, we conduct feature selection using forward selection (with default metric AIC) for the `PeakPosition` model and `WeeksOnChart` model.

```{r, echo=FALSE, warning=FALSE, message=FALSE}
selected_model_df <- data.frame(
  PeakCoefficients = c("Intercept", "Duration", "Danceability", "Energy", "Key", "Loudness", "Mode", "Speechiness", "Acousticness", "Instrumentalness", "Liveness", "Valence", "Tempo", "TimeSignature", "OtherGenre", "NoGenre", "Pop", "Country", "Rock", "Blues", "R&B", "RapHipHop", "Indie", "FamousCollaborator", "Collaborators"),
  PeakSelection = c("65.299", "—", "—", "8.397", "—", "—", "—", "—", "—", "—", "-1.213", "2.124", "—", "-1.640", "7.635", "—", "1.760", "13.981", "6.330", "—", "6.492", "3.369", "10.130", "—", "1.870"),
  WeeksSelection = c("-0.109", "0.002", "8.332", "—", "—", "0.143", "—", "-1.380", "-0.231", "-0.545", "-0.091", "3.544", "—", "—", "-1.166", "—", " 0.764", "2.185", "3.575", "—", "3.279", "—", "-1.662", "—", "—"))

colnames(selected_model_df) <- c("Variables", "PeakPosition Coefficient", "WeeksOnChart Coefficient")

knitr::kable(selected_model_df, align = 'c') %>%
  column_spec(1, width = "110px") %>%
  column_spec(2, width = "90px") %>%
  column_spec(3, width = "90px") 
```

```{r, echo=FALSE}
# Intercept-only model 
m0 <- lm(PeakPosition ~ 1, data = full_data)

# AIC used as score by default
forward <- step(m0, 
                scope=list(lower= ~ 1,
                           upper = ~ Duration + Danceability + Energy + Key + Loudness + Mode + Speechiness + Acousticness + Instrumentalness + Liveness + Valence + Tempo + TimeSignature + OtherGenre + NoGenre + Pop + Country + Rock + Blues + `R&B` + RapHipHop + Indie + FamousCollaborator + Collaborators),
                direction="forward", data=full_data, trace=0)
# forward
```

```{r, echo=FALSE}
# Intercept-only model 
m0 <- lm(WeeksOnChart ~ 1, data = full_data)

# AIC used as score by default
forward <- step(m0, 
                scope=list(lower= ~ 1,
                           upper = ~ Duration + Danceability + Energy + Key + Loudness + Mode + Speechiness + Acousticness + Instrumentalness + Liveness + Valence + Tempo + TimeSignature + OtherGenre + NoGenre + Pop + Country + Rock + Blues + `R&B` + RapHipHop + Indie + FamousCollaborator + Collaborators),
                direction="forward", data=full_data, trace=0)
# forward
```

We see that the subset of features chosen for `PeakPosition` is `Country`, `Rock`, `R&B`, `Collaborators`, `Energy`, `Indie`, `RapHipHop`, `OtherGenre`, `Pop`, `Liveness`, `Valence`, `Instrumentalness`, and `TimeSignature`. The subset of features selected for `WeeksOnChart` is `Rock`, `Valence`, `R&B`, `Duration`, `Speechiness`, `Danceability`, `Country`, `Acousticness`, `Pop`, `Liveness`, `Tempo`, `Loudness`, `Energy`, `TimeSignature`, and `Indie`.

## F-Tests

Using these selected features, we train new reduced models using the complements of these feature subsets for `PeakPosition` and `WeeksOnChart`, titled `reduced_base_peak_model` and `reduced_base_weeks_model` respectively.

```{r, echo=FALSE, warning=FALSE, message=FALSE}
full_data <- data[complete.cases(data), ]
reduced_base_peak_model <- lm(PeakPosition ~ Duration + Danceability + Key + Loudness + Mode + Speechiness + Acousticness + Tempo + NoGenre + Blues + FamousCollaborator, data = full_data)
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
full_data <- data[complete.cases(data), ]
reduced_base_weeks_model <- lm(WeeksOnChart ~ Key + Mode + Instrumentalness + OtherGenre + NoGenre + Blues + RapHipHop + FamousCollaborator + Collaborators, data = full_data)
```

Then, we conduct an F-test of the full models against the reduced models for `PeakPosition` and `WeeksOnChart` respectively.

```{r, echo=FALSE, warning=FALSE, message=FALSE}
selected_model_df <- data.frame(
  Variable = c("PeakPosition", "WeeksOnChart"),
  Degrees = c("13", "15"),
  F_Val = c("34.037", "32.33"), 
  P_Val = c("< 2.2e-16", "< 2.2e-16")
)

colnames(selected_model_df) <- c("Outcome Variable", "Degrees of Freedom", "F Value", "Probability > F")

knitr::kable(selected_model_df, align = 'c') %>%
  column_spec(1, width = "110px") %>%
  column_spec(2, width = "90px") %>%
  column_spec(3, width = "90px") %>%
  column_spec(4, width = "90px")
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
# anova(reduced_base_peak_model, full_base_peak_model)
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
# anova(reduced_base_weeks_model, full_base_weeks_model)
```

We see that for both `PeakPosition` and `WeeksOnChart`, the p-value is near 0 so we reject the null hypothesis. In other words, this result suggests that there is a statistically significant difference between the full model and reduced model and that the missing variables (the selected features) contribute significantly to explaining the variance in `PeakPosition` and `WeeksOnChart` beyond what is explained in the reduced models.

Thus, when considering all features, there is evidence to see that the most important features impacting `PeakPosition` are `Duration`, `Liveness`, `Valence`, `Instrumentalness`, and `TimeSignature` for audio features, `Country`, `Rock`, `R&B`, `RapHipHop`, `OtherGenre`, and `Pop` for genre, and `Collaborators` for collaboration. Notably, genre seems to be an important factor in reaching the top of the charts. Likewise, this test gives support to the conclusion that the most important features in determining chart longevity as `WeeksOnChart` are `Valence`, `Duration`, `Speechiness`, `Danceability`, `Acousticness`, `Liveness`, `Tempo`, `Loudness`, `Energy`, and `TimeSignature` for audio features, `Rock`, `R&B`, `Country`, `Pop`, and `Indie` for genre, and no features for collaboration. It is interesting to observe that the addition of a famous collaborator or the number of collaborators seems to play little role in impacting chart longevity.

## Analysis of Correlation

Before continuing to delve deeper into analysis for each of these selected features, we first check the correlation between each variable in this set to minimize redundancy. In order to do this, we generate the correlation matrix and plot the pair-wise correlations between selected features in a heatmap.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.cap="Heatmap for correlation between all variables."}
correlation_matrix <- cor(full_data[, c('Energy', 'Liveness', 'Valence', 'Instrumentalness', 'TimeSignature', 'Duration', 'Speechiness', 'Danceability', 'Acousticness', 'Tempo', 'Loudness', 'Country', 'Rock', 'R&B','Indie', 'RapHipHop', 'OtherGenre', 'Pop', 'Collaborators')])
heatmap(correlation_matrix, Colv = NA, Rowv = NA, cexRow = 0.4, cexCol = 0.4)

# correlation_matrix['Energy', 'Loudness']
# correlation_matrix['RapHipHop', 'Speechiness']
```

In observing the heatmap, we observe that the majority of the variables have low correlation, as seen by the light yellow and light golden shading thoughout the heatmap. There a couple of squares in the heatmap that stand out as being particularly dark, indicating high correlation.

The most obvious pair of highly correlated features is `Energy` and `Loudness`, which have a correlation of 0.713. Additionally, two squares that are slightly darker represent `RapHipHop` and `Speechiness` which have a correlation of 0.513. All other features have low to very moderate correlation with value 0.4 or lower.

Noting that energy overlaps with volume, we consider the impact of removing Loudness from the linear model (given that `Energy` encapsulates volume, but also more about the nature of the music). Similarly, given that `Speechiness` is a measure of the presence of human speech in a song, we are curious about the redundancy of this feature with regard to the genre feature of `RapHipHop`. These features only overlap in the `WeeksOnChart` model, so we conduct an F-test to examine if the model suffers as a result of removing the features `Speechiness` and `Loudness` from that model.

First, we see that when we conduct an F-test with a reduced model predicting WeeksOnChart removing the predictor `Loudness`, the resulting p-value is 0.0043. This is statistically significant, so we reject the null hypothesis and note that removing `Loudness` from the original model significantly reduces the quality of the model fit. Similarly, we repeat this F-test procedure with a reduced model predicting `WeeksOnChart` with the predictor `Speechiness` removed. The resulting p-value is 1.502e-09, which is extremely small and indicates that removing `Speechiness` likewise reduces the model fit significantly.

```{r, echo = FALSE, warning=FALSE, message=FALSE}
full_data <- data[complete.cases(data), ]
selected_weeks_model <- lm(WeeksOnChart ~ Rock + Valence + `R&B` + Duration + Speechiness + Danceability + Country + Acousticness + Pop + Liveness + Tempo + Loudness + Energy + TimeSignature + Indie, data = full_data)
# In reduced model, remove Loudness
a_reduced_selected_weeks_model <- lm(WeeksOnChart ~ Rock + Valence + `R&B` + Duration + Speechiness + Danceability + Country + Acousticness + Pop + Liveness + Tempo + Energy + TimeSignature + Indie, data = full_data)
# Remove Speechiness
b_reduced_selected_weeks_model <- lm(WeeksOnChart ~ Rock + Valence + `R&B` + Duration + Danceability + Country + Acousticness + Pop + Liveness + Tempo + Loudness + Energy + TimeSignature + Indie, data = full_data)
# anova(a_reduced_selected_weeks_model, selected_weeks_model)
```

```{r, echo = FALSE, warning=FALSE, message=FALSE}
# anova(b_reduced_selected_weeks_model, selected_weeks_model)
```

```{r, echo=FALSE, warning=FALSE, message=FALSE}
selected_model_df <- data.frame(
  Variable = c("PeakPosition (no Loudness)", "WeeksOnChart (no Speechiness)"),
  Degrees = c("1", "1"),
  F_Val = c("8.453", "38.977"), 
  P_Val = c("3.654e-03", "4.526e-10")
)

colnames(selected_model_df) <- c("Model", "Degrees of Freedom", "F Value", "Probability > F")

knitr::kable(selected_model_df, align = 'c') %>%
  column_spec(1, width = "120px") %>%
  column_spec(2, width = "90px") %>%
  column_spec(3, width = "90px") %>%
  column_spec(4, width = "90px")
```

Therefore, we see by the results of these F-tests that it makes sense to include `Speechiness` and `Loudness` in the model, and that these prediction variables contribute meaningfully in predicting `WeeksOnChart`. With this insight and noting the low correlation values for all other pairwise combinations of the selected features, we have confidence that this set of features is indeed significant and proceed with analyzing the individual features more deeply.

## Principal Component Analysis

To further validate the selected set of features, we run a principal component analysis for each linear model (with the selected features for `PeakPosition` and `WeeksOnChart` respectively, and plot the points along the first two principal components.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap = 'Principal component analysis for the PeakPosition variables.'}
library(ggfortify)
full_data <- data[complete.cases(data), ]
imp_data <- full_data[, c('Energy', 'Liveness', 'Valence', 'Instrumentalness', 'TimeSignature', 'Duration', 'Speechiness', 'Danceability', 'Acousticness', 'Tempo', 'Loudness', 'Country', 'Rock', 'R&B','Indie', 'RapHipHop', 'OtherGenre', 'Pop', 'Collaborators', 'PeakPosition', 'WeeksOnChart')] 

pca_res <- prcomp(scale(imp_data[, c('PeakPosition', 'Danceability', 'Country', 'Rock', 'R&B', 'Energy', 'Indie', 'Collaborators', 'OtherGenre', 'RapHipHop', 'Liveness', 'Pop', 'TimeSignature', 'Valence')]), scale. = FALSE)

# Display by color, hard to see
custom_breaks_1 <- c(0, 50, 100)
imp_data$PeakPosition <- cut(imp_data$PeakPosition, breaks = custom_breaks_1, labels = FALSE)
imp_data$PeakPosition <- factor(imp_data$PeakPosition, levels = c(1, 2), labels = c("Rank #1 - Rank #50",  'Rank #51 - Rank #100'))

autoplot(pca_res, data = imp_data, colour = 'PeakPosition', alpha = 0.5) +
  theme_minimal() + 
  labs(x = "Principal Component 1 (18.18%)", y = "Principal Component 2 (11.59%)") +  
  theme(
    panel.background = element_rect(fill = "white"),
  )
```

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.height=4, fig.width=6, fig.cap = 'Principal component analysis for the WeeksOnChart variables.'}
library(ggfortify)
full_data <- data[complete.cases(data), ]
imp_data <- full_data[, c('Energy', 'Liveness', 'Valence', 'Instrumentalness', 'TimeSignature', 'Duration', 'Speechiness', 'Danceability', 'Acousticness', 'Tempo', 'Loudness', 'Country', 'Rock', 'R&B','Indie', 'RapHipHop', 'OtherGenre', 'Pop', 'Collaborators', 'PeakPosition', 'WeeksOnChart')] 

pca_res <- prcomp(scale(imp_data[, c('WeeksOnChart', 'Danceability', 'Rock', 'R&B', 'Speechiness', 'Valence', 'Duration', 'Country', 'Liveness', 'Loudness', 'Pop', 'Instrumentalness', 'Acousticness', 'Indie', 'OtherGenre', 'TimeSignature')]), scale. = FALSE)

# Display by color, hard to see
custom_breaks_2 <- c(0, 26, 52, 100)
imp_data$WeeksOnChart <- cut(imp_data$WeeksOnChart, breaks = custom_breaks_2, labels = FALSE)
imp_data$WeeksOnChart <- factor(imp_data$WeeksOnChart, levels = c(1, 2, 3), labels = c('0-26 Weeks', '27-52 Weeks', '53+ Weeks'))

autoplot(pca_res, data = imp_data, colour = 'WeeksOnChart', alpha = 0.5) + 
  theme_minimal() + 
  labs(x = "Principal Component 1 (14.36%)", y = "Principal Component 2 (10.70%)") +  
  theme(
    panel.background = element_rect(fill = "white"),
  )
```

First, we examine the principal component analysis for `PeakPosition`. We see that while the points are relatively clustered together, we see that the higher ranked songs (from 1-50) cluster toward the left hand side of the plot while the lower ranked songs (from 51-100) cluster toward the right hand side of the plot. While there is not a clear separation between the two groups, we note a limitation in that the first two principal components only capture 18.18% and 11.59% of the variance respectively. This indicates high dimensionality in the data which may not be able to be easily expressed on such a two-dimensional plot. Thus, while we cannot actively conclude that this principal component analysis validates this set of variables as a significant one for prediction of peak position, the clustering we see (albeit with overlap) lends support to this set of variables.

Similarly, with the set of significant variables in determining WeeksOnChart, we see similar overlapping points but slight clustering behavior. In particular, in this PCA, the songs that stayed on the charts for over half a year (green and blue) cluster toward the bottom right side of the plot while the songs with shorter stays are more toward the center and top of the plot. Here, we run into a similar issue as the first two principal components only capture 14.36% and 10.7% of the variance respectively. Similarly, while acknowledging this limitation, we consider the visual presence of some grouping in classifying a song chart longevity to be support toward the selected set of features.

# Analysis of Significant Features

With our analysis above determining minimal redundancy and validating significance of the features in this set, we investigate each of these features with regard to the response variables. For each such variable, we consider **six** pieces of evidence to help us draw a corroborated conclusion about the relation between that song and virality as a whole:

**(1) Coefficient of variable in linear model predicting chart peak on the whole dataset.** We consider the coefficients of the linear model predicting `PeakPosition` from the set of selected variables specified above.

**(2) Coefficient of variable in linear model predicting chart longevity on the whole dataset.** We consider the coefficients of the linear model predicting `PeakPosition` from the set of selected variables specified above.

**(3) Coefficient of variable in linear models predicting chart peak, separated by year.** We partition the dataset by year, and then train 20 linear models (one for each year between 2000 and 2019) based on the selected features to predict `PeakPosition`. Then, we graph the change in coefficients with regard to the particular variable by year as a means of gauging how the magnitude and sign of the coefficients change year over year.

**(4) Coefficient of variable in linear models predicting chart longevity, separated by year.** From the same procedure above, we do the same for models predicting `WeeksOnChart` year over year.

**(5) Distribution of songs in top quartile of the variable, by year.** We consider the distribution of songs across years in the top quartile (taken across the whole dataset) of the given variable. The resulting plot informs whether more or fewer songs are trending toward having that variable be high as time passes.

**(6) Average variable value in top quartile, by year.** We consider the average variable value of the songs in the top quartile of that variable by year to observe changes in absolute value over time. For a given varible (ex. energy), this is to account for the scenarios where it could be possible that less and less songs have high energy, but the average value of the high energy songs is actually increasing.

```{r, echo=FALSE}
peak_00 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2000, ])
peak_01 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2001, ])
peak_02 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2002, ])
peak_03 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2003, ])
peak_04 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2004, ])
peak_05 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2005, ])
peak_06 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2006, ])
peak_07 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2007, ])
peak_08 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2008, ])
peak_09 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2009, ])
peak_10 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2010, ])
peak_11 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2011, ])
peak_12 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2012, ])
peak_13 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2013, ])
peak_14 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2014, ])
peak_15 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2015, ])
peak_16 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2016, ])
peak_17 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2017, ])
peak_18 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2018, ])
peak_19 <- lm(formula = PeakPosition ~ Country + Rock + `R&B` + Energy + 
    Indie + Collaborators + OtherGenre + RapHipHop + Pop + Liveness + 
    Acousticness + Valence + Instrumentalness, data = full_data[full_data$Year == 2019, ])

peak_coef_data <- data.frame(
  coef(peak_00),
  coef(peak_01), 
  coef(peak_02),
  coef(peak_03),
  coef(peak_04),
  coef(peak_05),
  coef(peak_06),
  coef(peak_07),
  coef(peak_08),
  coef(peak_09),
  coef(peak_10),
  coef(peak_11),
  coef(peak_12), 
  coef(peak_13),
  coef(peak_14),
  coef(peak_15),
  coef(peak_16),
  coef(peak_17),
  coef(peak_18),
  coef(peak_19))
peak_coef_data <- t(peak_coef_data)
```

```{r, echo = FALSE}
weeks_00 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2000, ])
weeks_01 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2001, ])
weeks_02 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2002, ])
weeks_03 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2003, ])
weeks_04 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2004, ])
weeks_05 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2005, ])
weeks_06 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2006, ])
weeks_07 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2007, ])
weeks_08 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2008, ])
weeks_09 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2009, ])
weeks_10 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2010, ])
weeks_11 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2011, ])
weeks_12 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2012, ])
weeks_13 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2013, ])
weeks_14 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2014, ])
weeks_15 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2015, ])
weeks_16 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2016, ])
weeks_17 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2017, ])
weeks_18 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2018, ])
weeks_19 <- lm(WeeksOnChart ~ Rock + `R&B` + Danceability + Speechiness + 
    Duration + Valence + Country + Acousticness + Liveness + 
    Pop + Indie + TimeSignature + Tempo + OtherGenre + Energy + 
    Loudness + RapHipHop, data = full_data[full_data$Year == 2019, ])

weeks_coef_data <- data.frame(
  coef(weeks_00),
  coef(weeks_01), 
  coef(weeks_02),
  coef(weeks_03),
  coef(weeks_04),
  coef(weeks_05),
  coef(weeks_06),
  coef(weeks_07),
  coef(weeks_08),
  coef(weeks_09),
  coef(weeks_10),
  coef(weeks_11),
  coef(weeks_12), 
  coef(weeks_13),
  coef(weeks_14),
  coef(weeks_15),
  coef(weeks_16),
  coef(weeks_17),
  coef(weeks_18),
  coef(weeks_19))
weeks_coef_data <- t(weeks_coef_data)
```

As a general note for this section, keep in mind that `PeakPosition` has the highest ranking at 1 and the lowest ranking at 100, such that lower/negative coefficients are predictive of higher rankings (and not lower ones). Also, note that some features will not have all of the six pieces of evidence described above as certain models do not have certain features and for the binary genre features analysis of quartiles is less relevant (as all values are either 0 or 1).

# Insights About Song Characteristic

## Energy

For this variable, we note that the average energy level has dropped as time has gone on, with lower levels in more recent years.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average energy levels for songs by year."}
energy_data <- data %>% 
  select(SongID, Month, Year, Energy) %>% 
  filter(!is.na(Energy) & !is.nan(Energy))

counts_by_year <- energy_data %>%
  count(Year, name = "Total_Count")

average_energy_per_year <- energy_data %>%
  group_by(Year) %>%
  summarise(Average_Energy = mean(Energy, na.rm = TRUE))

# Plot the histogram
ggplot(average_energy_per_year, aes(x = Year, y = Average_Energy)) +
  geom_col(fill='blue') +
  labs(title = "Average Energy per Year",
       x = "Year",
       y = "Average Energy") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

For `Energy`, the coefficient in predicting `PeakPosition` is 9.800 and the coefficient in predicting `WeeksOnChart` is -3.780. This indicates that for the model across the 20 year period, higher energy values were predictive of lower chart rankings and shorter total stays on the chart.

Examining the development of coefficients over time, we see that in the early 2000s, higher energy values indicated higher chart performance but shorter stays on the charts. Then, in the late 2000s to 2015, there was a reversal as higher energy values predict lower chart peaks, with still shorter stays on the charts. From 2015 onwards, there is one final reversal as the trend reverts back to higher energy values being tied to higher peak rankings as well as an upward trend with regard to chart longevity.

In examining the distribution of songs in the top energy quartile over time, we clearly see a trend in an increasing percentage of high energy songs from 2000 to 2012 followed by a drop from 2013 to 2017, and a gradual upward curve from 2018 to 2019. The plot on the right side indicates that over time, the level of energy in "high energy" songs on the charts has been steadily decreasing.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Distribution and average energy per year only for songs in the top energy quartile across the two-decade dataset."}
peak_energy_coef <- peak_coef_data[, "Energy"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Energy_Coef = peak_energy_coef)

energy_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Energy_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Energy Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "Energy Coefficient") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

weeks_energy_coef <- weeks_coef_data[, "Energy"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Energy_Coef = weeks_energy_coef)

energy_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Energy_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Energy Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "Energy Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

energy_data <- data %>% 
  select(SongID, Month, Year, Energy) %>% 
  filter(!is.na(Energy) & !is.nan(Energy))

quartiles <- quantile(energy_data$Energy, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
energy_data$Quartile <- cut(energy_data$Energy, breaks = quartiles, labels = FALSE)

counts_by_year <- energy_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- energy_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_energy_data <- energy_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_energy_data %>%
  group_by(Year) %>%
  summarise(Average_Energy = mean(Energy, na.rm = TRUE))

energy_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in\nTop Energy Quartile", title = "Percentage of Songs per Year in\nTop Energy Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

energy_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Energy)) +
  geom_point() +
  labs(x = "Year", y = "Average Energy in Annual\nTop Quartile", title = "Average Energy in Annual Top\nQuartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "gam", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

grid.arrange(energy_coef_1, energy_coef_2, energy_plot1, energy_plot2, ncol = 2)
```

## Valence

For the PeakPosition model, the coefficient for Valence is 2.1387, meaning that songs with higher musical positiveness generally indicated lower chart peaks across the time period. Meanwhile, the coefficient for Valence in the WeeksOnChart model is 4.161, indicating that higher valence songs predicted longer stays.

An analysis of the PeakPosition plot on the left over time shows that from 2000 to 2005, high valence songs were predictive of songs that had progressively lower chart peaks. Then, in 2006, the trend has since been steady in high valence being indicative of increasingly higher chart peaks. Meanwhile, the sinusoidal trend shape for the WeeksOnChart model shows that in recent years, high valence has been indicative of increasingly shorter stays on the charts.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average valence levels for songs by year."}
valence_data <- data %>% 
  select(SongID, Month, Year, Valence) %>% 
  filter(!is.na(Valence) & !is.nan(Valence))

counts_by_year <- valence_data %>%
  count(Year, name = "Total_Count")

average_valence_per_year <- valence_data %>%
  group_by(Year) %>%
  summarise(Average_Valence = mean(Valence, na.rm = TRUE))

# Plot the histogram
ggplot(average_valence_per_year, aes(x = Year, y = Average_Valence)) +
  geom_col(fill='blue') +
  labs(title = "Average Valence per Year",
       x = "Year",
       y = "Average Valence") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

The distribution of high valence songs over the years shows that over time, there has been a steady decrease in the percentage of songs each year that make the charts with valence in the top quartile. The plot on the right side adds to this observation, showing that within the high valence songs, the songs are becoming lower in valence with regard to value as time passes.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.fullwidth=FALSE, fig.cap="Distribution and average valence per year only for songs in the top valence quartile across the two-decade dataset."}

peak_valence_coef <- peak_coef_data[, "Valence"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Valence_Coef = peak_valence_coef)

valence_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Valence_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Valence Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "Valence Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

weeks_valence_coef <- weeks_coef_data[, "Valence"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Valence_Coef = weeks_valence_coef)

valence_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Valence_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Valence Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "Valence Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

valence_data <- data %>% 
  select(SongID, Month, Year, Valence) %>% 
  filter(!is.na(Valence) & !is.nan(Valence))

quartiles <- quantile(valence_data$Valence, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
valence_data$Quartile <- cut(valence_data$Valence, breaks = quartiles, labels = FALSE)

counts_by_year <- valence_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- valence_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_valence_data <- valence_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_valence_data %>%
  group_by(Year) %>%
  summarise(Average_Valence = mean(Valence, na.rm = TRUE))

valence_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in\nTop Valence Quartile", title = "Percentage of Songs per Year\nin Top Valence Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

valence_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Valence)) +
  geom_point() +
  labs(x = "Year", y = "Average Valence in Annual\nTop Quartile", title = "Average Valence in Annual Top\nQuartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "lm", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

grid.arrange(valence_coef_1, valence_coef_2, valence_plot1, valence_plot2, ncol = 2)
```

## Danceability

The Danceability variable, only found in the WeeksOnChart model across two decades, has a coefficient of 7.001. This indicates that throughout the entire dataset, higher danceability scores generally predicted a higher number of weeks a song stays on the Hot 100.

Taking a closer look at the change in coefficients over time shows a more nuanced story: namely, the danceability coefficient in predicting WeeksOnChart increased steadily from 2000 to 2012, and then has been on a steady decline since. This indicates that from 2000 to 2012, songs deemed to be very danceable were increasingly indicative of high success in staying on the chart for a long time. However, since 2013, danceable songs have become predictive of increasingly shorter chart stays.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average danceability levels for songs by year."}
danceability_data <- data %>% 
  select(SongID, Month, Year, Danceability) %>% 
  filter(!is.na(Danceability) & !is.nan(Danceability))

counts_by_year <- danceability_data %>%
  count(Year, name = "Total_Count")

average_danceability_per_year <- danceability_data %>%
  group_by(Year) %>%
  summarise(Average_Danceability = mean(Danceability, na.rm = TRUE))

# Plot the histogram
ggplot(average_danceability_per_year, aes(x = Year, y = Average_Danceability)) +
  geom_col(fill='blue') +
  labs(title = "Average Danceability per Year",
       x = "Year",
       y = "Average Danceability") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

In examining the distributions of songs over time with high danceability values, we can clearly observe that the concentration of very danceable songs declined from 2000 to 2011, at which point the trend reversed and began to steadily increase. As seen on the right hand side, the average values also mirrored this trend. This means that since 2011, danceable songs represent an increasing proportion of charting songs each year. Moreover, in this time period the danceable songs have only gotten more danceable with regard to value.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Distribution of average danceability per year only for songs in the top danceability quartile across the two-decade dataset."}
weeks_dance_coef <- weeks_coef_data[, "Danceability"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Dance_Coef = weeks_dance_coef)

dance_coef <- ggplot(weeks_plot_data, aes(x = Year, y = Dance_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Danceability Coefficient Over the Years (WeeksOnChart Model)", x = "Year", y = "Danceability Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

danceability_data <- data %>% 
  select(SongID, Month, Year, Danceability) %>% 
  filter(!is.na(Danceability) & !is.nan(Danceability))

quartiles <- quantile(danceability_data$Danceability, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE, duplicates.ok = TRUE)
danceability_data$Quartile <- ifelse(danceability_data$Danceability >= quartiles[4], 4, 0)

counts_by_year <- danceability_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- danceability_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_danceability_data <- danceability_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_danceability_data %>%
  group_by(Year) %>%
  summarise(Average_Danceability = mean(Danceability, na.rm = TRUE))

danceability_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in\nTop Danceability Quartile", title = "Percentage of Songs per Year in\nTop Danceability Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

danceability_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Danceability)) +
  geom_point() +
  labs(x = "Year", y = "Average Danceability in Annual\nTop Quartile", title = "Average Danceability in Annual\nTop Quartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "loess", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))
grid.arrange(dance_coef, danceability_plot1, danceability_plot2, layout_matrix = layout_matrix)
```

## Tempo

For Tempo, the coefficient with regard to this variable on the overall model for WeeksOnChart is -7.128. Given that this coefficient is negative, we see that throughout 2000-2019, a slower song (lower tempo value) is predictive of a greater number of weeks on the Hot 100.

This relationship is amplified over time, as the downward trend starting from 2007 shows steadily decreasing tempo coefficients that indicate that lower tempos increasingly predict longer reigns on the charts.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average tempo for songs by year."}
tempo_data <- data %>% 
  select(SongID, Month, Year, Tempo) %>% 
  filter(!is.na(Tempo) & !is.nan(Tempo))

counts_by_year <- tempo_data %>%
  count(Year, name = "Total_Count")

average_tempo_per_year <- tempo_data %>%
  group_by(Year) %>%
  summarise(Average_Tempo = mean(Tempo, na.rm = TRUE))

# Plot the histogram
ggplot(average_tempo_per_year, aes(x = Year, y = Average_Tempo)) +
  geom_col(fill='blue') +
  labs(title = "Average Tempo per Year",
       x = "Year",
       y = "Average Tempo") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

The figure on the left shows that, for two periods 2000-2009 and 2012-2019, very fast songs make up an increasing percentage of charting songs per year. In contrast, as time goes by, the fast songs that make the charts are becoming steadily slower.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Distribution and average tempo per year only for songs in the top energy quartile across the two-decade dataset."}
weeks_tempo_coef <- weeks_coef_data[, "Tempo"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Tempo_Coef = weeks_tempo_coef)

tempo_coef <- ggplot(weeks_plot_data, aes(x = Year, y = Tempo_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Tempo Coefficient Over the Years (WeeksOnChart Model)", x = "Year", y = "Tempo Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

tempo_data <- data %>% 
  select(SongID, Month, Year, Tempo) %>% 
  filter(!is.na(Tempo) & !is.nan(Tempo))

quartiles <- quantile(tempo_data$Tempo, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE, duplicates.ok = TRUE)
tempo_data$Quartile <- ifelse(tempo_data$Tempo >= quartiles[4], 4, 0)

counts_by_year <- tempo_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- tempo_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_tempo_data <- tempo_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_tempo_data %>%
  group_by(Year) %>%
  summarise(Average_Tempo = mean(Tempo, na.rm = TRUE))

tempo_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in Top\nTempo Quartile", title = "Percentage of Songs per Year in Top\nTempo Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

tempo_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Tempo)) +
  geom_point() +
  labs(x = "Year", y = "Average Tempo in Annual\nTop Quartile", title = "Average Tempo in Annual Top\nQuartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "lm", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))
grid.arrange(tempo_coef, tempo_plot1, tempo_plot2, layout_matrix = layout_matrix)
```

## Duration

The Duration variable, which is only in the WeeksOnChart model over all two decades, has a coefficient of 0.0002 (which is small in magnitude due to the large values of duration in milliseconds). This positive coefficient indicates that longer songs are indicative of better success in terms of longevity on the charts throughout the time period.

However, an examination of duration by year shows a steady increase from 2000 to 2008, and then a steady decrease from 2009 to 2019 in terms of Duration coefficients. This trend indicates that the duration of a song predicted longer stays on the charts to an increasing degree until 2008, at which point shorter songs indicated longevity to a greater extent as time went on.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average duration of songs by year."}
duration_data <- data %>% 
  select(SongID, Month, Year, Duration) %>% 
  filter(!is.na(Duration) & !is.nan(Duration))

counts_by_year <- duration_data %>%
  count(Year, name = "Total_Count")

average_duration_per_year <- duration_data %>%
  group_by(Year) %>%
  summarise(Average_Duration = mean(Duration, na.rm = TRUE))

# Plot the histogram
ggplot(average_duration_per_year, aes(x = Year, y = Average_Duration)) +
  geom_col(fill='blue') +
  labs(title = "Average Duration per Year",
       x = "Year",
       y = "Average Duration") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

Looking at the distribution of songs by year in the top quartile of longest songs, it is clear to see a decreasing trend over time. In other words, as we move from 2000 to 2019, fewer of the songs in each year that make the charts are very long. Additionally, the figure on the right shows us that even within the long songs that make it onto the charts in recent years, those long songs are becoming shorter.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Distribution and average duration per year only for songs in the top duration quartile across the two-decade dataset."}
weeks_duration_coef <- weeks_coef_data[, "Duration"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Duration_Coef = weeks_duration_coef)

duration_coef <- ggplot(weeks_plot_data, aes(x = Year, y = Duration_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Duration Coefficient Over the Years (WeeksOnChart Model)", x = "Year", y = "Duration Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

duration_data <- data %>% 
  select(SongID, Month, Year, Duration) %>% 
  filter(!is.na(Duration) & !is.nan(Duration))

quartiles <- quantile(duration_data$Duration, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE, duplicates.ok = TRUE)
duration_data$Quartile <- ifelse(duration_data$Duration >= quartiles[4], 4, 0)

counts_by_year <- duration_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- duration_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_duration_data <- duration_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_duration_data %>%
  group_by(Year) %>%
  summarise(Average_Duration = mean(Duration, na.rm = TRUE))

duration_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in Top\nDuration Quartile", title = "Percentage of Songs per Year in\nTop Duration Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

duration_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Duration)) +
  geom_point() +
  labs(x = "Year", y = "Average Duration in Annual\nTop Quartile", title = "Average Duration in Annual Top\nQuartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "loess", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(duration_coef, duration_plot1, duration_plot2, layout_matrix = layout_matrix)
```

# Insights About Song Production

## Instrumentalness

The Instrumentalness variable is not included in the WeeksOnChart model, and in the PeakPosition model over the whole time period, its coefficient is 5.812. This positive value indicates that throughout the twenty years, songs with high scores in Instrumentalness have been indicative of lower chart peaks.

Further, examining the change in this coefficient year over year, we see a slight trend upward with a steady increase in years 2017-2019. This trend indicates that as we move through time, high instrumentalness becomes predictive of increasingly lower peak rankings on the charts.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average instrumentalness levels for songs by year."}
instrumentalness_data <- data %>% 
  select(SongID, Month, Year, Instrumentalness) %>% 
  filter(!is.na(Instrumentalness) & !is.nan(Instrumentalness))

counts_by_year <- instrumentalness_data %>%
  count(Year, name = "Total_Count")

average_instrumentalness_per_year <- instrumentalness_data %>%
  group_by(Year) %>%
  summarise(Average_Instrumentalness = mean(Instrumentalness, na.rm = TRUE))

# Plot the histogram
ggplot(average_instrumentalness_per_year, aes(x = Year, y = Average_Instrumentalness)) +
  geom_col(fill='blue') +
  labs(title = "Average Instrumentalness per Year",
       x = "Year",
       y = "Average Instrumentalness") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

Examining the distribution of songs with high instrumentalness across the years, there is a clear downward trend in the percentage of songs released each year with high instrumentalness values. However, as seen in the plot on the right side, starting from 2015, the average value of the instrumentalness scores is increasing. This means that although it is lesson common to see a song with high instrumentalness reach the charts in the present-day, the highly instrumental songs that do make it on the charts are more instrumental than they were in previous years.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Distribution and average log instrumentalness per year only for songs in the top instrumentalness quartile across the two-decade dataset."}
peak_inst_coef <- peak_coef_data[, "Instrumentalness"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Inst_Coef = peak_inst_coef)

inst_coef <- ggplot(peak_plot_data, aes(x = Year, y = Inst_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "log(Instrumentalness) Coefficient Over the Years (PeakPosition Model)", x = "Year", y = "log(Instrumentalness) Coefficient") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

instrumentalness_data <- data %>% 
  select(SongID, Month, Year, Instrumentalness) %>% 
  filter(!is.na(Instrumentalness) & !is.nan(Instrumentalness))

quartiles <- quantile(instrumentalness_data$Instrumentalness, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE, duplicates.ok = TRUE)
instrumentalness_data$Quartile <- ifelse(instrumentalness_data$Instrumentalness >= quartiles[4], 4, 0)

counts_by_year <- instrumentalness_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- instrumentalness_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_instrumentalness_data <- instrumentalness_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_instrumentalness_data %>%
  group_by(Year) %>%
  summarise(Average_Instrumentalness = mean(Instrumentalness, na.rm = TRUE))

inst_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in Top\nInstrumentalness Quartile", title = "Percentage of Songs per Year in\nTop Instrumentalness Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

inst_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Instrumentalness)) +
  geom_point() +
  labs(x = "Year", y = "Average Instrumentalness in\nAnnual Top Quartile", title = "Average Instrumentalness in\nAnnual Top Quartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "loess", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(inst_coef, inst_plot1, inst_plot2, layout_matrix = layout_matrix)
```

## Acousticness

The variable measuring Acousticness is not included in the PeakPosition overall model, but has a coefficient of -3.810 on the overall WeeksOnChart model. This negative coefficient indicates that a higher acousticness score in this twenty year period is generally predictive of worse outcomes with regard to chart longevity.

Taking a look at the trend over time shows a steady decrease in coefficients from 2001 to 2013, at which point the trend becomes less clear. This indicates that in this 12 year period, acoustic songs become increasingly predictive of less longevity.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average acousticness levels for songs by year."}
acousticness_data <- data %>% 
  select(SongID, Month, Year, Acousticness) %>% 
  filter(!is.na(Acousticness) & !is.nan(Acousticness))

counts_by_year <- acousticness_data %>%
  count(Year, name = "Total_Count")

average_acousticness_per_year <- acousticness_data %>%
  group_by(Year) %>%
  summarise(Average_Acousticness = mean(Acousticness, na.rm = TRUE))

# Plot the histogram
ggplot(average_acousticness_per_year, aes(x = Year, y = Average_Acousticness)) +
  geom_col(fill='blue') +
  labs(title = "Average Acousticness per Year",
       x = "Year",
       y = "Average Acousticness") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

Notably, the percentage of songs per year that are classified as highly acoustic (in the top quartile) has been increasing from 2013 onwards. There is also a very clear positive and linear relationship in the plot for average values, showing that of the highly acoustic songs that chart each year, as time goes on these songs are becoming acoustic to a greater extent.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Distribution and average log acousticness per year only for songs in the top acousticness quartile across the two-decade dataset."}
weeks_ac_coef <- weeks_coef_data[, "Acousticness"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Ac_Coef = weeks_ac_coef)

ac_coef <- ggplot(weeks_plot_data, aes(x = Year, y = Ac_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "log(Acousticness) Coefficient Over the Years (WeeksOnChart Model)", x = "Year", y = "log(Acousticness) Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

acousticness_data <- data %>% 
  select(SongID, Month, Year, Acousticness) %>% 
  filter(!is.na(Acousticness) & !is.nan(Acousticness))

quartiles <- quantile(acousticness_data$Acousticness, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE, duplicates.ok = TRUE)
acousticness_data$Quartile <- ifelse(acousticness_data$Acousticness >= quartiles[4], 4, 0)

counts_by_year <- acousticness_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- acousticness_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_acousticness_data <- acousticness_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_acousticness_data %>%
  group_by(Year) %>%
  summarise(Average_Acousticness = mean(Acousticness, na.rm = TRUE))

ac_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in Top\nAcousticness Quartile", title = "Percentage of Songs per Year in Top\nAcousticness Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

ac_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Acousticness)) +
  geom_point() +
  labs(x = "Year", y = "Average Acousticness in Annual\nTop Quartile", title = "Average Acousticness in Annual Top\nQuartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "lm", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 


layout_matrix <- rbind(c(1, 1),
                       c(2, 3))
grid.arrange(ac_coef, ac_plot1, ac_plot2, layout_matrix = layout_matrix)
```

## Loudness

The coefficient for the Loudness variable int he comprehensive WeeksOnChart model is 2.375, which suggests that louder songs across all years are more predictive of better chart longevity.

Examining the coefficients over the years gives a marked upward trend over time. This trend shows that as time goes by, loud songs become predictive of successful chart longevity to an increasing extent.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE,  fig.cap="Average loudness levels for songs by year."}
loudness_data <- data %>% 
  select(SongID, Month, Year, Loudness) %>% 
  filter(!is.na(Loudness) & !is.nan(Loudness))

counts_by_year <- loudness_data %>%
  count(Year, name = "Total_Count")

average_loudness_per_year <- loudness_data %>%
  group_by(Year) %>%
  summarise(Average_Loudness = mean(Loudness, na.rm = TRUE))

# Plot the histogram
ggplot(average_loudness_per_year, aes(x = Year, y = Average_Loudness)) +
  geom_col(fill='blue') +
  labs(title = "Average Loudness per Year",
       x = "Year",
       y = "Average Loudness") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

Meanwhile, viewing the distribution of loud songs by year shows that the percentage of very loud charting songs per year steadily increased from 2000 to 2009, decreased from 2010 to 2017, and begins to increase again from 2018 onwards. The figure plotting the average loudness value for loud songs shows a mirrored trend, meaning that when the percentage of loud charting songs increases, the loud songs get increasingly louder (and vice versa).

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width=7, fig.height=5,fig.cap="Distribution and average loudness per year only for songs in the top loudness quartile across the two-decade dataset."}

weeks_loud_coef <- weeks_coef_data[, "Loudness"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Loud_Coef = weeks_loud_coef)

loud_coef <- ggplot(weeks_plot_data, aes(x = Year, y = Loud_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Loudness Coefficient Over the Years (WeeksOnChart Model)", x = "Year", y = "Loudness Coefficient") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

loudness_data <- data %>% 
  select(SongID, Month, Year, Loudness) %>% 
  filter(!is.na(Loudness) & !is.nan(Loudness))

quartiles <- quantile(loudness_data$Loudness, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE, duplicates.ok = TRUE)
loudness_data$Quartile <- ifelse(loudness_data$Loudness >= quartiles[4], 4, 0)

counts_by_year <- loudness_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- loudness_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_loudness_data <- loudness_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_loudness_data %>%
  group_by(Year) %>%
  summarise(Average_Loudness = mean(Loudness, na.rm = TRUE))

loudness_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in Top\nLoudness Quartile", title = "Percentage of Songs per Year in Top\nLoudness Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

loudness_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Loudness)) +
  geom_point() +
  labs(x = "Year", y = "Average Loudness in Annual\nTop Quartile", title = "Average Loudness in Annual Top\nQuartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "loess", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(loud_coef, loudness_plot1, loudness_plot2, layout_matrix = layout_matrix)
```

## Liveness

For `Liveness`, the coefficient for the overall peak position model is -4.503, indicating that across the 20 year period liveness is indicative of higher chart performance. Further, the coefficient for the overall `WeeksOnChart` model is -3.150, meaning that songs with higher liveness scores across the time period generally had shorter overall stays on the charts.

In looking at the plot on the left side below, we clearly see that the coefficients for liveness across models for PeakPosition increased steadily from 2000 to 2012, and then decreased steadily from 2013 to 2019. This trend indicates that as we approach the present day, high liveness scores in songs are increasingly predictive of high rankings. The WeeksOnChart model supports this conclusion, as alhthough there is less of a clear trend, the points are clearly moving upward starting from 2016 such that songs with high liveness scores are more indicative of long stays on the charts in recent years.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average liveness levels for songs by year."}
liveness_data <- data %>% 
  select(SongID, Month, Year, Liveness) %>% 
  filter(!is.na(Liveness) & !is.nan(Liveness))

counts_by_year <- liveness_data %>%
  count(Year, name = "Total_Count")

average_liveness_per_year <- liveness_data %>%
  group_by(Year) %>%
  summarise(Average_Liveness = mean(Liveness, na.rm = TRUE))

# Plot the histogram
ggplot(average_liveness_per_year, aes(x = Year, y = Average_Liveness)) +
  geom_col(fill='blue') +
  labs(title = "Average Liveness per Year",
       x = "Year",
       y = "Average Liveness") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

In examining the distribution of songs with high liveness over the years, it is clear to see th at there is a trend moving toward a higher concentration of songs with high liveness from 2000 to 2009, and then a relatively steady decrease from 2010 to 2019. The plot on the right adds to this data, showing that alongside this timeline the value of liveness in high liveness songs has followed a similar trend.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 7, fig.height = 5, fig.cap="Distribution and average log liveness per year only for songs in the top liveness quartile across the two-decade dataset."}

peak_liveness_coef <- peak_coef_data[, "Liveness"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Liveness_Coef = peak_liveness_coef)

liveness_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Liveness_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "log(Liveness) Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "log(Liveness) Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

weeks_liveness_coef <- weeks_coef_data[, "Liveness"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Liveness_Coef = weeks_liveness_coef)

liveness_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Liveness_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "log(Liveness) Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "log(Liveness) Coefficient") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

liveness_data <- data %>% 
  select(SongID, Month, Year, Liveness) %>% 
  filter(!is.na(Liveness) & !is.nan(Liveness))

quartiles <- quantile(liveness_data$Liveness, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
liveness_data$Quartile <- cut(liveness_data$Liveness, breaks = quartiles, labels = FALSE)

counts_by_year <- liveness_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- liveness_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

q4_liveness_data <- liveness_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_liveness_data %>%
  group_by(Year) %>%
  summarise(Average_Liveness = mean(Liveness, na.rm = TRUE))

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

liveness_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in Top\nLiveness Quartile", title = "Percentage of Songs per Year in\nTop Liveness Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

liveness_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Liveness)) +
  geom_point() +
  labs(x = "Year", y = "Average Liveness in Annual\nTop Quartile", title = "Average Liveness in Annual Top\nQuartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "loess", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 7), axis.text.y = element_text(size = 7)) 

grid.arrange(liveness_coef_1, liveness_coef_2, liveness_plot1, liveness_plot2, ncol = 2)
```

## Speechiness

The coefficient for Speechiness in the WeeksOnChart model in the 2000-2019 period (not included in the PeakPosition model) is -7.129. This negative coefficient indicates that within this time period, songs that were speechy were predicted to perform more poorly in terms of staying on the charts.

We gain additional insight from looking at the coefficients for Speechiness over time. In the plot below, we see that there is a steady upward trend starting from 2000 to 2005, and then less of a clear trend following. This observation shows that songs with lots of human speech became predictive of increasingly long stays on the charts until 2005, at which point the indication varies per year.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.margin=TRUE, fig.cap="Average speechiness for songs by year."}
speechiness_data <- data %>% 
  select(SongID, Month, Year, Speechiness) %>% 
  filter(!is.na(Speechiness) & !is.nan(Speechiness))
counts_by_year <- speechiness_data %>%
  count(Year, name = "Total_Count")

average_speechiness_per_year <- speechiness_data %>%
  group_by(Year) %>%
  summarise(Average_Speechiness = mean(Speechiness, na.rm = TRUE))

# Plot the histogram
ggplot(average_speechiness_per_year, aes(x = Year, y = Average_Speechiness)) +
  geom_col(fill='blue') +
  labs(title = "Average Speechiness per Year",
       x = "Year",
       y = "Average Speechiness") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme_minimal() +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
```

Furthermore, as we look at the distribution of songs with high speechiness values (such as rap songs) across the years, we clearly see an increase in the concentration of very speechy songs around 2003, a dip between 2006 and 2014, and then a steady increase up to 2019. The figure on the right matches this trend as well. This shows that leading up to 2003 and following 2014, not only were speechy song indicative of increasingly longer stays on the charts, but they were also becoming increasingly speechy on average.

```{r, echo=FALSE, message = FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Distribution and average speechiness per year only for songs in the top speechiness quartile across the two-decade dataset."}

weeks_speechiness_coef <- weeks_coef_data[, "Speechiness"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Speechiness_Coef = weeks_speechiness_coef)

speech_coef <- ggplot(weeks_plot_data, aes(x = Year, y = Speechiness_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "log(Speechiness) Coefficient Over the Years (WeeksOnChart Model)", x = "Year", y = "log(Speechiness) Coefficient") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

speechiness_data <- data %>% 
  select(SongID, Month, Year, Speechiness) %>% 
  filter(!is.na(Speechiness) & !is.nan(Speechiness))

quartiles <- quantile(speechiness_data$Speechiness, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE, duplicates.ok = TRUE)
speechiness_data$Quartile <- ifelse(speechiness_data$Speechiness >= quartiles[4], 4, 0)

counts_by_year <- speechiness_data %>%
  count(Year, name = "Total_Count")

counts_quartile_4 <- speechiness_data %>%
  filter(Quartile == 4) %>%
  count(Year, name = "Count_Quartile_4")

merged_data <- left_join(counts_by_year, counts_quartile_4, by = "Year")

merged_data <- merged_data %>%
  mutate(Percentage_Quartile_4 = (Count_Quartile_4 / Total_Count) * 100)

q4_speechiness_data <- speechiness_data %>%
  filter(Quartile == 4) 

q4_by_year <- q4_speechiness_data %>%
  group_by(Year) %>%
  summarise(Average_Speechiness = mean(Speechiness, na.rm = TRUE))

sp_plot1 <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Quartile_4)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs in Top\nSpeechiness Quartile", title = "Percentage of Songs per Year in\nTop Speechiness Quartile") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 


sp_plot2 <- ggplot(q4_by_year, aes(x = Year, y = Average_Speechiness)) +
  geom_point() +
  labs(x = "Year", y = "Average Speechiness in Annual\nTop Quartile", title = "Average Speechiness in Annual\nTop Quartile by Year") + 
  theme_minimal() + 
  geom_smooth(method = "loess", se = FALSE) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1)) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(speech_coef, sp_plot1, sp_plot2, layout_matrix = layout_matrix)
```

# Insights About Genre

## Pop

For the overall PeakPosition model, the coefficient is 1.798, and for the overall WeeksOnChart model, the coefficient for Pop is 0.692. This indicates that throughout the time period, a song being labeled as pop indicates lower peaks on the charts, but positively contributes to the prediction of the number of weeks on the charts.

Taking a look at the coefficients over time, we observe that there is a steady increase in the coefficients for predicting PeakPosition, and a steady decrease for the coefficient of Pop in predicting the chart longevity. This observation provides evidence that as time passes, a song labeled in the genre of pop is predictive of lower chart peak and a shorter duration of stay to a greater degree.

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=7, fig.height=5, fig.cap="Change in coefficient for pop variable in the linear models trained to predict peak position and total weeks on the charts."}

total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

pop_songs <- data %>%
  group_by(Year) %>%
  summarise(Pop_Songs = sum(Pop == 1))

merged_data <- merge(total_songs, pop_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_Pop = (Pop_Songs / Total_Songs) * 100)

pop_plot <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Pop)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Pop Songs", title = "Percentage of Pop Songs by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

peak_pop_coef <- peak_coef_data[, "Pop"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Pop_Coef = peak_pop_coef)

pop_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Pop_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Pop Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "Pop Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

weeks_pop_coef <- weeks_coef_data[, "Pop"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Pop_Coef = weeks_pop_coef)

pop_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Pop_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Pop Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "Pop Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(pop_plot, pop_coef_1, pop_coef_2, layout_matrix = layout_matrix)
```

It is also important to consider that the percentage of charting songs that are in the pop genre steadily increased from 2000 to 2012, but has been steadily decreasing since

## Country

Across the 2000-2019 time period, the coefficient for Country is 13.959 in predicting peak ranking and 2.290 in predicting weeks on the chart. We note that the coefficient 13.959 is relatively large (especially compared to Pop), and its positive nature indicates that a song labeled as country is predicted to be have much more mediocre peak rankings. However, the positive value for WeeksOnChart indicates positive prediction for longer chart stays across the entire period.

With further inspection year-by-year, we see that over time, there is a clear trend with increasing coefficients in the PeakPosition model, and increasing coefficitions in the WeeksOnChart model since 2015 (with a spike prior from 2005 - 2014). This combination of evidence suggests that over time, the label of being a country genre is predictive to an increasing degree of low chart performance, but since 2014, the label more so indicates long stays on the charts.

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=7, fig.height=5, fig.cap="Change in the coefficient for the country variable in the linear models predicting peak position and total weeks on the charts."}

total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

country_songs <- data %>%
  group_by(Year) %>%
  summarise(Country_Songs = sum(Country == 1))

merged_data <- merge(total_songs, country_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_Country = (Country_Songs / Total_Songs) * 100)

country_plot <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Country)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Country Songs", title = "Percentage of Country Songs by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

peak_country_coef <- peak_coef_data[, "Country"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Country_Coef = peak_country_coef)

country_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Country_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Country Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "Country Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

weeks_country_coef <- weeks_coef_data[, "Country"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Country_Coef = weeks_country_coef)

country_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Country_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Country Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "Country Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(country_plot, country_coef_1, country_coef_2, layout_matrix = layout_matrix)
```

Regarding the distribution of country songs over time, we see that there is a relatively high concentration of charting country songs in the early 2000s (roughly 1 in 5 songs on the charts were tagged country) which decreased until 2010, had a resurgence from 2011 to 2013, and then dipped again to the present day which sees an increase in 2019. Overall, the percentage of country songs per year for charting songs decreased through this two decade period, but seems to be increasing again due to the significantly higher percentage in 2019.

## Rock

The coefficient for Rock is 6.382 in the PeakPosition overall model, and 3.558 in the overall WeeksOnChart model. These coefficients indicate that across 2000-2019, songs in the genre Rock are predicted to reach lower peaks but have longer stays on the Billboard Hot 100.

Delving into the year-by-year coefficient charts shows a clear increase in coefficients for PeakPosition until 2013, followed by a steep decrease. The WeeksOnChart figure follows a similar pattern from 2013 onwards, with a steady decrease in coefficients. This indicates that from 2013, the label of Rock became less indicative of a low peak, but more indicative of a short stay on the charts.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width=7, fig.height=5, fig.cap="Change in coefficient for rock variable in the linear models trained to predict peak position and total weeks on chart."}

total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

rock_songs <- data %>%
  group_by(Year) %>%
  summarise(Rock_Songs = sum(Rock == 1))

merged_data <- merge(total_songs, rock_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_Rock = (Rock_Songs / Total_Songs) * 100)

rock_plot <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Rock)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Rock Songs", title = "Percentage of Rock Songs by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

peak_rock_coef <- peak_coef_data[, "Rock"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Rock_Coef = peak_rock_coef)

rock_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Rock_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Rock Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "Rock Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

weeks_rock_coef <- weeks_coef_data[, "Rock"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Rock_Coef = weeks_rock_coef)

rock_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Rock_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Rock Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "Rock Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

grid.arrange(rock_plot, rock_coef_1, rock_coef_2, layout_matrix = layout_matrix)
```

The distribution of the proportion of rock songs on the charts per year shows that since 2005, there has been a notable and steady decrease in the percentage of charting rock songs. Notably, in 2018 and 2019, less than 8% of the songs that reached the charts were tagged as any type of rock

## R&B

For the overall set of the data, the coefficient for R&B was 6.561 in predicting PeakPosition and 3.134 in predicting WeeksOnChart. These coefficients suggest that the R&B label is predictive of lower chart peaks and higher number of weeks on the chart in general.

To expand on this insight, we take a look at coefficients over time. We see that the coefficient increases over time in predicting the peak position on the charts, and also seems to increase (albeit less steadily) in predicting longevity. In turn, these observations indicate an increasing extent of the general trend observed in the overall model: over time, an R&B label indicates lower peaks and longer stays to a greater degree.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 7, fig.height = 5, fig.cap="Change in coefficient for R&B variable in the linear models trained to predict peak position and total weeks on charts."}
total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

rnb_songs <- data %>%
  group_by(Year) %>%
  summarise(RnB_Songs = sum(`R&B` == 1))

merged_data <- merge(total_songs, rnb_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_RnB = (RnB_Songs / Total_Songs) * 100)

rnb_plot <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_RnB)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of R&B Songs", title = "Percentage of R&B Songs by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

peak_rnb_coef <- peak_coef_data[, "`R&B`"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, RnB_Coef = peak_rnb_coef)

rnb_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = RnB_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "R&B Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "R&B Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

weeks_rnb_coef <- weeks_coef_data[, "`R&B`"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, RnB_Coef = weeks_rnb_coef)

rnb_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = RnB_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "R&B Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "R&B Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(rnb_plot, rnb_coef_1, rnb_coef_2, layout_matrix = layout_matrix)
```

In taking a look at the distribution of the percentage of R&B songs per year, we see a steady decline as time passes (with the exception of a short rise from 2013 to 2016).

## Indie

For the Indie label on the entire dataset, the coefficient in the PeakPosition model is 10.124 and the coefficient in the WeeksOnChart model is -1.768. These observations indicate that for all years, generally an Indie classification is predictive of lower peaks and shorter stays.

The decrease from year 2015 onward in the coefficients for the PeakPosition model shows that the impact on the peak being lower is reduced as time moves forward. Further, the rise in coefficients for the WeeksOnChart model from 2005 to 2011 hits a plateau until the end of the time period, showing that moving more toward recent times Indie records have not been as predictive of very short chart stays.

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width = 7, fig.height = 5, fig.cap="Change in coefficient for indie variable in the linear models trained to predict peak position and total weeks on chart."}

total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

indie_songs <- data %>%
  group_by(Year) %>%
  summarise(Indie_Songs = sum(Indie == 1))

merged_data <- merge(total_songs, indie_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_Indie = (Indie_Songs / Total_Songs) * 100)

indie_plot <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Indie)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Indie Songs", title = "Percentage of Indie Songs by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 
peak_indie_coef <- peak_coef_data[, "Indie"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Indie_Coef = peak_indie_coef)

indie_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Indie_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Indie Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "Indie Coefficient") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

weeks_indie_coef <- weeks_coef_data[, "Indie"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Indie_Coef = weeks_indie_coef)

indie_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Indie_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Indie Coefficient Over the Years\n(WeeksOnChart Model)", x = "Year", y = "Indie Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(indie_plot, indie_coef_1, indie_coef_2, layout_matrix = layout_matrix)
```

Note that the percentage of indie songs out of all chart songs is relatively low, with a maximum of \~5% in 2006. In observing the variation between years, we see that although the proportion of charting indie songs is low in the early 2000s and near 2020, there are spikes and falls in a relatively random manner between 2005 and 2017.

## RapHipHop

From 2000 to 2020, the coefficients for Rap in the models are 3.416 for PeakPosition and -5.662 for WeeksOnChart. Incorporating songs throughout the entire two decade period, in general a rap song is predicted to have a slightly lower peak and shorter stay on the charts.

However, analyzing the time series representation gives us further insights. We see a notable increase in the coefficients for determining chart peaks from 2004 to 2014, and then a rapid drop since then. These trends indicate that as time approaches 2020, rap songs are increasingly predicted to chart highly. The plot for the WeeksOnChart model also shows a steady rise in coefficients until 2013, and then a rapid descent since then. This indicates an increase in the predictive strength for Rap for longer chart stays as we approach the current day.

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width=7, fig.height=5, fig.cap="Change in coefficient for rap/hip-hop variable in the linear models trained to predict peak position and total weeks on chart."}

total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

raphiphop_songs <- data %>%
  group_by(Year) %>%
  summarise(RapHipHop_Songs = sum(RapHipHop == 1))

merged_data <- merge(total_songs, raphiphop_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_RapHipHop = (RapHipHop_Songs / Total_Songs) * 100)

rhh_plot <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_RapHipHop)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Rap and Hip Hop Songs", title = "Percentage of Rap and Hip Hop Songs by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

peak_rhh_coef <- peak_coef_data[, "RapHipHop"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, RHH_Coef = peak_rhh_coef)

rhh_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = RHH_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "RapHipHop Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "RapHipHop Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

weeks_rhh_coef <- weeks_coef_data[, "RapHipHop"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, RHH_Coef = weeks_rhh_coef)

rhh_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = RHH_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "RapHipHop Coefficient Over the\nYears (WeeksOnChart Model)", x = "Year", y = "RapHipHop Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(rhh_plot, rhh_coef_1, rhh_coef_2, layout_matrix = layout_matrix)
```

Taking a look at the distribution of rap and hip-hop songs as a proportion of the total number of charting songs each year, we clearly see an increase over time that accelerates in 2014 and peaks at 2018 (with a dip in 2019).

## OtherGenre

We consider the OtherGenre variable, which represents a song without any of the other established genre labels seen above (intuitively obscure songs). The coefficient for the overall model in predicting peak position is 7.550, and the coefficient for predicting chart longevity is -1.281. These coefficients support the notion that this label is indicative of lower peaks and shorter stays.

In looking at the OtherGenre coefficient values over time, we see a similar trend in both models: a dip starting in 2000 that increases to a peak in 2009, a subsequent dip that reaches a bottom in 2014, and an increase from that point forward until 2019. Focusing on the last increase in coefficients from 2014, this period of time has this label increasingly be predictive of worse peaks and less successful reigns.

```{r, echo=FALSE, warning=FALSE, message=FALSE, fig.width = 7, fig.height = 5, fig.cap="Change in coefficient for OtherGenre in the linear models trained to predict peak position and weeks on the chart."}

total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

othergenre_songs <- data %>%
  group_by(Year) %>%
  summarise(OtherGenre_Songs = sum(OtherGenre == 1))

merged_data <- merge(total_songs, othergenre_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_OtherGenre = (OtherGenre_Songs / Total_Songs) * 100)

other_plot <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_OtherGenre)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Other Genre\nSongs", title = "Percentage of Other Genre Songs by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

peak_other_coef <- peak_coef_data[, "OtherGenre"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Other_Coef = peak_other_coef)

other_coef_1 <- ggplot(peak_plot_data, aes(x = Year, y = Other_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "OtherGenre Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "OtherGenre Coefficient") +
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

weeks_other_coef <- weeks_coef_data[, "OtherGenre"]
years <- 2000:2019
weeks_plot_data <- data.frame(Year = years, Other_Coef = weeks_other_coef)

other_coef_2 <- ggplot(weeks_plot_data, aes(x = Year, y = Other_Coef)) +
  geom_point() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "OtherGenre Coefficient Over the\nYears (WeeksOnChart Model)", x = "Year", y = "OtherGenre Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(other_plot, other_coef_1, other_coef_2, layout_matrix = layout_matrix)
```

The percentage of other genre songs per year is highest in the early 2000s, and has relatively steadily decreased since then. In more recent years, roughly 5-10% of songs are not categorized by any of the other existing labels.

# Insights About Collaboration

In the two-decade model for determining peak position, the coefficient for the variable Collaborators counting the number of collaborators on each charting track is 1.888. This positive value suggests that an increased number of collaborators is predictive of decreased chart peaks.

Looking further at this variable across each year tells us that the coefficients increased in the early 2000s, and has been decreasing as a trend since then (particularly in the more recent years). This trend indicates that in the years that are closer to 2020, a higher number of collaborators becomes less extreme in indicating lower peak performance on the Hot 100.

```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.width = 7, fig.height = 5, fig.cap="Change in coefficient for number of collaborators in the linear model trained to predict peak position."}
collab_data <- data %>% 
  select(SongID, Month, Year, Collaborators) %>% 
  filter(!is.na(Collaborators) & !is.nan(Collaborators))

counts_by_year <- collab_data %>%
  count(Year, name = "Total_Count")

average_collab_per_year <- collab_data %>%
  group_by(Year) %>%
  summarise(Average_Collab = mean(Collaborators, na.rm = TRUE))

collab_plot <- ggplot(average_collab_per_year, aes(x = Year, y = Average_Collab)) +
  geom_col(fill='blue') +
  labs(title = "Average Collaborators per Year",
       x = "Year",
       y = "Average Collaborators") +
  theme_minimal() + 
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

total_songs <- data %>%
  group_by(Year) %>%
  summarise(Total_Songs = n())

collab_songs <- data %>%
  group_by(Year) %>%
  summarise(Collab_Songs = sum(Collaborators >= 1))

merged_data <- merge(total_songs, collab_songs, by = "Year", all.x = TRUE)

merged_data <- merged_data %>%
  mutate(Percentage_Collab = (Collab_Songs / Total_Songs) * 100)

collab_plot_more <- ggplot(merged_data, aes(x = factor(Year), y = Percentage_Collab)) +
  geom_bar(stat = "identity", fill = "blue") +
  labs(x = "Year", y = "Percentage of Songs with\n1+ Collaborators", title = "Percentage of Songs with 1+\nCollaborators by Year") + 
  theme_minimal() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

peak_collab_coef <- peak_coef_data[, "Collaborators"]
years <- 2000:2019
peak_plot_data <- data.frame(Year = years, Collab_Coef = peak_collab_coef)

collab_peak_plot <- ggplot(peak_plot_data, aes(x = Year, y = Collab_Coef)) +
  geom_point() +
  geom_smooth(method = "loess", se = FALSE) +
  labs(title = "Collaborators Coefficient Over the Years\n(PeakPosition Model)", x = "Year", y = "Collaborators Coefficient") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 12), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 9), axis.text.y = element_text(size = 9)) 

layout_matrix <- rbind(c(1, 1),
                       c(2, 3))

grid.arrange(collab_plot, collab_plot_more, collab_peak_plot, layout_matrix = layout_matrix)
```

Aside from songs with no collaborators, the most common number of collaborators to have on a track is 1. In examining this distribution of songs on the charts with exactly 1 collaborator, we note that this type of collaboration was far less frequent in the 2000s, but steadily grew in popularity over time until the present day.

# Results

After conducting this analysis of my data with regard to feature selection, statistical tests, and visual exploration, a number of insights have become clear.

1.  The features and patterns that impact longevity on the charts (WeeksOnChart) are totally different from those that determine peak ranking on the charts (PeakPosition). This can be seen by the different features selected for these different response variables, contrasting trends shown in the plots for the same variable, and contrasting trends over time for variables in the same set.
2.  The most significant features in creating chart hits that peak highly pertain to the genre, production, and energetic character of the song. This is seen by evidence from feature selection, results from F-tests, EDA plots, and later significant feature analysis plots.
3.  The most significant features in creating chart hits that last a long time pertain to the genre, sensual impact (ex. loudness, tempo), and production of the tracks. This insight is grounded by the selected features, the exploratory data analysis and distributions, the results from the F-tests, and the more detailed time analysis above.
4.  Collaboration seems to have a relatively minimal impact on overall chart performance, and is a fairly common practice. This can be seen by the exclusion of collaboration variables in feature selection, verification by F-tests, and the plots for the collaboration variables.
5.  The impact of different features change significantly over time, often in well-defined trends, in a way that impacts both chart peak and longevity (often together but not always). This is mostly seen int he later analysis with year-by-year data of each significant feature in the Analysis of Significant Features section.

# Corroborating Resources

These results are corroborated by Viner's 2020 analysis [^14] of the impact of audio features on Billboard chart data on a smaller timeframe (2010 to 2020) and with relatively elementary analysis (simple averaging of audio features). Although we pursue different methodologies, these respective analyses identify qualities relating to the energetic character and production of a song to be most critical with regard to its success. Furthermore, the aforementioned article by Pham et al.[^15] similarly supports the findings above, echoing the conclusion that the most important features in determining chart performance lie in audio features.

[^14]: [Josh Viner](https://joshdviner.medium.com/what-makes-a-hit-song-analyzing-data-from-the-billboard-hot-100-chart-74c1d5ad3fa3) (2020) *What Makes a Hit Song: Analyzing Data from the Billboard Hot 100 Chart*, Medium

[^15]: [Pham et al.](https://cs229.stanford.edu/proj2015/140_report.pdf) (2015) *Predicting Song Popularity*, Stanford Department of Computer Science

In addition, another relevant source to compare results regarding my hypothesis is from the Spotify Newsroom[^16] with regard to the impact of collaboration with famous artists on the success of songs. This article is published from data analysts at Spotify, and analyzes a number of top hits throughout the desired time period with further context about each individual song. As seen by the plots and trends demonstrated for the selected significant features, the insights found by Spotify closely mirror those in my analysis. These results add another dimension of nuance to the impact of collaboration on the success of top songs, and help to prove my hypothesis more convincingly.

[^16]: [Spotify For The Record](https://newsroom.spotify.com/2023-05-24/crossover-collaborations-genre-collabs-streaming-data-spotify/) (2023) *The Crossover Effect: Artist Collaborations Thrive on Spotify*, Spotify

# Conclusion

## Limitations

While the results found have been corroborated by multiple sources with in-depth analysis taking multiple variables into consideration, this analysis is inherently limited by three main factors.

1.  Although the source of this combined dataset is Spotify and Billboard, large and reputable companies, there are still many other avenues through which people listen to music. Relying on data from Spotify and Billboard means that music more streamed through other platforms (ex. Apple Music, Soundcloud) is not considered in the analysis. In particular, there may be a bias toward established artists, as niche/indie music websites (ex. Soundcloud) or music by smaller artists (ex. covers on YouTube) are likely to be undercounted by Spotify/Billboard.
2.  While the dataset I analyzed includes a more recent time period (2000-2019) than other comparable studies, my analysis is limited to music before 2019. Thus, though it is likely the trends found in my analysis are likely to continue in the present day, there is not a strict guarantee that this analysis can be effectively extrapolated to draw conclusions about the music after 2020 without further work on recent years.
3.  In using audio features that are unilaterally determined by Spotify, these scores are subjective in favor of Spotify's aims and at risk of being biased. This is particularly true for scores for subjective features like energy or danceability as compared to objectively determinable features such as tempo and song duration. In particular, this is a notable concern because the principal component analysis conducted shows that while there is a degree of clustering, there is some overlap that can point to redundancy or other flaws inherent in the features.

## Ethical Considerations

Noting the limitations of this analysis, we consider the critical ethical issues that may stem as a result.

1.  In particular, relying on Spotify and Billboard for analysis may cause representative harms in creating a recipe for popular music. Given that Spotify and Billboard are companies with operations centered in America and mostly working in English, it is likely that these platforms are more accessible in North America than in other places worldwide. Consequently, the scores and rankings found in the dataset are likely to be more representative of the tastes and preferences of urban Western listeners who have better access to Spotify. This means that the results of this analysis is likely biased toward English songs which are popular in America, suggesting that artists emulate this style of music to achieve success on the charts. This may be harmful with regard to failing to adequately represent the listening tastes and musical creations of artists around the world.
2.  With this analysis being focused on music only from 2000-2019, the data and resulting analysis may inherently reflect discriminatory tastes prevalent during this period of time. As we move closer to the present day, the music industry has become more diverse with the inclusion of global music, unique singers, and united efforts against problematic musicians and producers (ex. Ke\$ha's producer, Dr. Luke). Excluding the popular music from 2020 onwards means excluding the musical tastes of a more conscious population of listeners and may potentially amplify the work of troubled musicians and their beliefs.
3.  By relying on Spotify's subjective scores for audio features, this analysis is vulnerable to echoing and amplifying any existing biases within these features. For example, if the models used to determine the scores within Spotify are trained mostly on American music, then the lower scores assigned to other diverse music will influence the findings in this analysis to discourage emulating other styles of music. Clearly, this can lead to the minimization of music from other countries, cultures, and backgrounds which is a critical representative harm.

## Future Research and Broader Significance

To address the limitations and mitigate the ethical concerns discussed above, future work should be done in a similar methodology to this report that incorporates music from 2020 onwards from multiple sources of data beyond Spotify and Billboard. Ideally, a future analysis includes top songs aggregating the features from multiple music libraries across different countries up to 2024 (and beyond).

Despite these limitations, the many modes of statistical analysis with consideration of multiple features and support from corroborating research shows that the results from this analysis holds insights that can be applied as a broad recipe for making a hit song in the current era. These insights are extremely significant in allowing up-and-coming artists make their first hit to secure their livelihood, boost their career, and share the story behind their music. Additionally, by sharing this information about top hits widely, we even the playing field in the music industry to allow smaller artists to succeed while introducing listeners to more diverse music representative of different cultures and styles.\
\newpage

# References

[Sean Miller](https://www.kaggle.com/datasets/thedevastator/billboard-hot-100-audio-features?select=Hot+Stuff.csv) (2024) *Billboard Hot 100 Weekly Charts with Spotify Audio Features*, Kaggle

[Mohamed Nasreldin](https://towardsdatascience.com/song-popularity-predictor-1ef69735e380) (2018) *Song Popularity Predictor*, Medium

[Pham et al.](https://cs229.stanford.edu/proj2015/140_report.pdf) (2015) *Predicting Song Popularity*, Stanford Department of Computer Science

[Sean Miller](https://www.kaggle.com/datasets/thedevastator/billboard-hot-100-audio-features?select=Hot+Stuff.csv) (2024) *Billboard Hot 100 Weekly Charts with Spotify Audio Features*, Kaggle

[Billboard](https://www.billboard.com/billboard-charts-legend/) (2024) *Billboard Charts Legend: Recurrent Rules*, Billboard

[R. Kelly](https://open.spotify.com/track/45LrQ3tg1z8plpuQRCuSwE?si=ff2b73f4e7c548ba) (2005) *Trapped In The Closet*, Spotify

[SongBPM](https://songbpm.com/the-black-eyed-peas/imma-be) *Song Metrics*, Imma Be

[Josh Viner](https://joshdviner.medium.com/what-makes-a-hit-song-analyzing-data-from-the-billboard-hot-100-chart-74c1d5ad3fa3) (2020) *What Makes a Hit Song: Analyzing Data from the Billboard Hot 100 Chart*, Medium

[Pham et al.](https://cs229.stanford.edu/proj2015/140_report.pdf) (2015) *Predicting Song Popularity*, Stanford Department of Computer Science

[Spotify For The Record](https://newsroom.spotify.com/2023-05-24/crossover-collaborations-genre-collabs-streaming-data-spotify/) (2023) *The Crossover Effect: Artist Collaborations Thrive on Spotify*, Spotify

\newpage

# Appendix

In the appendix below, you can find the histograms and boxplots for the following variables:

-   Duration

-   Danceability

-   Energy

-   Key

-   Loudness

-   Valence

-   Tempo

-   Time Signature

-   Speechiness (with and without log transformation)

-   Acousticness (with and without log transformation)

-   Instrumentalness (with and without log transformation)

-   Liveness (with and without log transformation)

-   Time Signature

-   Pop

-   Country

-   Blues

-   R&B

-   Rap and Hip-Hop

-   Indie

-   Rock

-   OtherGenre

-   Mode

-   Famous Collaborator

-   Number of Collaborators

\newpage

```{r, echo=FALSE, warning=FALSE, message=FALSE, out.height="100%", fig.width=8, fig.height=15, fig.cap="Distribution of songs by duration, danceability, energy, key, loudness, valence, tempo, and time signature."}
library(ggplot2)
library(patchwork)

filtered_data <- data %>% filter(!is.na(Duration) & !is.nan(Duration))
duration_dist <- ggplot(filtered_data, aes(x = Duration)) +
  geom_histogram(binwidth = 30000, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Duration (ms)", x = "Duration (ms)", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Danceability) & !is.nan(Danceability))
danceability_dist <- ggplot(filtered_data, aes(x = Danceability)) +
  geom_histogram(binwidth = 0.05, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Danceability", x = "Danceability Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Energy) & !is.nan(Energy))
energy_dist <- ggplot(filtered_data, aes(x = Energy)) +
  geom_histogram(binwidth = 0.05, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Energy", x = "Energy Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Key) & !is.nan(Key))
filtered_data$Key <- factor(filtered_data$Key)
key_dist <- ggplot(filtered_data, aes(x = Key)) +
  geom_bar(fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Key", x = "Key", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Loudness) & !is.nan(Loudness))
loudness_dist <- ggplot(filtered_data, aes(x = Loudness)) +
  geom_histogram(binwidth = 1, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Loudness (dB)", x = "Loudness (dB)", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Valence) & !is.nan(Valence))
valence_dist <- ggplot(filtered_data, aes(x = Valence)) +
  geom_histogram(binwidth = 0.05, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Valence", x = "Valence Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Tempo) & !is.nan(Tempo))
tempo_dist <- ggplot(filtered_data, aes(x = Tempo)) +
  geom_histogram(binwidth = 10, fill = "orange", color = "black") +#
  labs(title = "Distribution of Songs by Tempo (bpm)", x = "Tempo (beats per min)", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(TimeSignature) & !is.nan(TimeSignature))
ts_dist <- ggplot(filtered_data, aes(x = TimeSignature)) +
  geom_histogram(binwidth = 1, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Time Signature", x = "Time Signature (beats per measure)", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

grid.arrange(duration_dist, danceability_dist, energy_dist, key_dist, loudness_dist, valence_dist, tempo_dist, ts_dist, ncol = 2)
```

\newpage

```{r, echo=FALSE, warning=FALSE, message=FALSE, out.height="100%", fig.width=8, fig.height=15, fig.cap="Distribution of songs by speechiness, acousticness, instrumentalness, and liveness with the associated log transformations."}
library(ggplot2)
library(patchwork)

path <- "/Users/ericlee/Desktop/Spring 2024/DATASCI 120/"
data <- read.csv(paste0(path, "DATASCI 120 - Processed Data.csv"))
names(data)[names(data) == 'Peak.Position'] <- 'PeakPosition'
names(data)[names(data) == 'Weeks.on.Chart'] <- 'WeeksOnChart'
names(data)[names(data) == 'Genre.List'] <- 'GenreList'
names(data)[names(data) == 'Duration..ms.'] <- 'Duration'
names(data)[names(data) == 'Time.Signature'] <- 'TimeSignature'
names(data)[names(data) == 'R.B'] <- 'R&B'
names(data)[names(data) == 'Rap.Hip.Hop'] <- 'RapHipHop'
names(data)[names(data) == 'Other.Genre'] <- 'OtherGenre'
names(data)[names(data) == 'No.Genre'] <- 'NoGenre'
names(data)[names(data) == 'Collaborators.List'] <- 'CollaboratorsList'
names(data)[names(data) == 'Famous.Collaborator'] <- 'FamousCollaborator'
data <- data %>%
  filter(Year != 2020)

# This code chunk includes the four audio features that were log transformed 
filtered_data <- data %>% filter(!is.na(Speechiness) & !is.nan(Speechiness))
speechiness_dist <- ggplot(filtered_data, aes(x = Speechiness)) +
  geom_histogram(binwidth = 0.05, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Speechiness", x = "Speechiness Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

log_speechiness_dist <- ggplot(filtered_data, aes(x = log(Speechiness))) +
  geom_histogram(binwidth = 0.25, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by log(Speechiness)", x = "log(Speechiness) Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Acousticness) & !is.nan(Acousticness))
acousticness_dist <- ggplot(filtered_data, aes(x = Acousticness)) +
  geom_histogram(binwidth = 0.05, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Acousticness", x = "Acousticness Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Acousticness) & !is.nan(Acousticness))
log_acousticness_dist <- ggplot(filtered_data, aes(x = log(Acousticness))) +
  geom_histogram(binwidth = 0.5, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by log(Acousticness)", x = "log(Acousticness) Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Instrumentalness) & !is.nan(Instrumentalness))
instrumentalness_dist <- ggplot(filtered_data, aes(x = Instrumentalness)) +
  geom_histogram(binwidth = 0.05, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Instrumentalness", x = "Instrumentalness Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Instrumentalness) & !is.nan(Instrumentalness))
log_instrumentalness_dist <- ggplot(filtered_data, aes(x = log(Instrumentalness))) +
  geom_histogram(binwidth = 0.5, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by log(Instrumentalness)", x = "log(Instrumentalness) Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Liveness) & !is.nan(Liveness))
liveness_dist <- ggplot(filtered_data, aes(x = Liveness)) +
  geom_histogram(binwidth = 0.05, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by Liveness", x = "Liveness Score", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Liveness) & !is.nan(Liveness))
log_liveness_dist <- ggplot(filtered_data, aes(x = log(Liveness))) +
  geom_histogram(binwidth = 0.3, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by log(Liveness)", x = "log(Liveness Score)", y = "Frequency") + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

grid.arrange(speechiness_dist, log_speechiness_dist, acousticness_dist, log_acousticness_dist, instrumentalness_dist, log_instrumentalness_dist, liveness_dist, log_liveness_dist, ncol = 2)
```

\newpage

```{r, echo=FALSE, warning=FALSE, message=FALSE, out.height="100%", fig.width=8, fig.height=15, fig.cap="Boxplots of peak position and chart longevity against time signature, pop genre, country genre, and blues genre."}
ts_filtered_data <- data %>% filter(!is.na(TimeSignature) & !is.nan(TimeSignature) & TimeSignature != 0)
ts_labels <- c("1", "3", "4", "5")

# Check unique levels of TimeSignature
unique_ts <- unique(ts_filtered_data$TimeSignature)

ts_plot1 <- ggplot(ts_filtered_data, aes(x = factor(TimeSignature, labels = ts_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs Time Signature", x = 'Time Signature (beats per measure)', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() + 
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))
  
ts_plot2 <- ggplot(ts_filtered_data, aes(x = factor(TimeSignature, labels = ts_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Time Signature", x = 'Time Signature (beats per measure)', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))


filtered_data <- data %>% filter(!is.na(Pop) & !is.nan(Pop))
pop_labels <- c('Not Pop', 'Pop')
pop_plot1 <- ggplot(filtered_data, aes(x = factor(Pop, labels = pop_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs Pop Genre", x = 'Pop Genre', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

pop_plot2 <- ggplot(filtered_data, aes(x = factor(Pop, labels = pop_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Pop Genre", x = 'Pop Genre', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Country) & !is.nan(Country))
country_labels <- c('Not Country', 'Country')
country_plot1 <- ggplot(filtered_data, aes(x = factor(Country, labels = country_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs Country Genre", x = 'Country Genre', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

country_plot2 <- ggplot(filtered_data, aes(x = factor(Country, labels = country_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Country Genre", x = 'Country Genre', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Blues) & !is.nan(Blues))
blues_labels <- c('Not Blues', 'Blues')
blues_plot1 <- ggplot(filtered_data, aes(x = factor(Blues, labels = blues_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs Blues Genre", x = 'Blues Genre', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

blues_plot2 <- ggplot(filtered_data, aes(x = factor(Blues, labels = blues_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Blues Genre", x = 'Blues Genre', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

grid.arrange(ts_plot1, ts_plot2, pop_plot1, pop_plot2, country_plot1, country_plot2, blues_plot1, blues_plot2, ncol = 2)
```

\newpage

```{r, echo=FALSE, warning=FALSE, message=FALSE, out.height="100%", fig.width=8, fig.height=15, fig.cap="Boxplots of peak position and chart longevity against R&B genre, rap/hip-hop genre, indie genre, and rock genre."}
filtered_data <- data %>% filter(!is.na(`R&B`) & !is.nan(`R&B`))
rnb_labels <- c('Not R&B', 'R&B')
rnb_plot1 <- ggplot(filtered_data, aes(x = factor(`R&B`, labels = rnb_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs R&B Genre", x = 'R&B Genre', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))


rnb_plot2 <- ggplot(filtered_data, aes(x = factor(`R&B`, labels = rnb_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs R&B Genre", x = 'R&B Genre', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(RapHipHop) & !is.nan(RapHipHop))
rhh_labels <- c('Not Rap/HipHop', 'Rap/HipHop')
rhh_plot1 <- ggplot(filtered_data, aes(x = factor(RapHipHop, labels = rhh_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs Rap/Hip-Hop Genre", x = 'Rap and Hip-Hop Genre', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

rhh_plot2 <- ggplot(filtered_data, aes(x = factor(RapHipHop, labels = rhh_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Rap/Hip-Hop Genre", x = 'Rap and Hip-Hop Genre', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Indie) & !is.nan(Indie))
indie_labels <- c('Not Indie', 'Indie')
ind_plot1 <- ggplot(filtered_data, aes(x = factor(Indie, labels = indie_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs Indie Genre", x = 'Indie Genre', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

ind_plot2 <- ggplot(filtered_data, aes(x = factor(Indie, labels = indie_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Indie Genre", x = 'Indie Genre', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Rock) & !is.nan(Rock))
rock_labels <- c('Not Rock', 'Rock')

rock_plot1 <- ggplot(filtered_data, aes(x = factor(Rock, labels = rock_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots in Peak Position vs Rock Genre", x = 'Rock Genre', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

rock_plot2 <- ggplot(filtered_data, aes(x = factor(Rock, labels = rock_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Rock Genre", x = 'Rock Genre', y = 'Longevity (Weeks on Chart)') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

grid.arrange(rnb_plot1, rnb_plot2, rhh_plot1, rhh_plot2, ind_plot1, ind_plot2, rock_plot1, rock_plot2, ncol = 2)
```

\newpage

```{r, echo=FALSE, warning=FALSE, message=FALSE, out.height="100%", fig.width=8, fig.height=15, fig.cap="Boxplots of peak position and chart longevity against less popular genres (OtherGenre), modality, and presence of famous collaborators."}

filtered_data <- data %>% filter(!is.na(OtherGenre) & !is.nan(OtherGenre))
othergenre_labels <- c("Known Genre", "Other Genre")
othergenre_plot1 <- ggplot(filtered_data, aes(x = factor(OtherGenre, labels = othergenre_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs OtherGenre", x = 'Genre Label', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

othergenre_plot2 <- ggplot(filtered_data, aes(x = factor(OtherGenre, labels = othergenre_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Longevity vs OtherGenre", x = 'Presence of OtherGenre', y = 'Longevity (Weeks on Chart)') +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))


filtered_data <- data %>% filter(!is.na(FamousCollaborator) & !is.nan(FamousCollaborator))
famous_labels <- c("None", "Famous Collaborators")
faco_plot1 <- ggplot(filtered_data, aes(x = factor(FamousCollaborator, labels = famous_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'skyblue') + 
  labs(title = "Boxplots of Peak Position vs Famous Collaborator", x = 'Presence of Famous Collaborator', y = 'Peak Position') + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  scale_y_reverse() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

faco_plot2 <- ggplot(filtered_data, aes(x = factor(FamousCollaborator, labels = famous_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Longevity vs Famous Collaborator", x = 'Presence of Famous Collaborator', y = 'Longevity (Weeks on Chart)') +
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Collaborators) & !is.nan(Collaborators))
collab_dist <- ggplot(filtered_data, aes(x = Collaborators)) +
  geom_histogram(binwidth = 1, fill = "orange", color = "black") +
  labs(title = "Distribution of Songs by # of Collaborators", x = "# of Collaborators", y = "Frequency") + 
  theme_bw() + 
  scale_x_continuous(breaks = seq(min(filtered_data$Collaborators), max(filtered_data$Collaborators), by = 1)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

mode_colors <- c("0" = "lightblue", "1" = "gold")
mode_labels <- c("0" = "Minor", "1" = "Major")
mode_counts <- table(data$Mode)
mode_data <- data.frame(Mode = factor(names(mode_counts)), Frequency = mode_counts)
mode_dist <- ggplot(mode_data, aes(x = "", y = Frequency.Freq, fill = Mode)) +
  scale_fill_manual(values = mode_colors, labels = mode_labels) +
  geom_bar(stat = "identity") +
  coord_polar("y") +
  labs(title = "Modality of Songs on Hot 100 Charts in 2000-2020", fill = "Mode") +
  theme_void() +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

filtered_data <- data %>% filter(!is.na(Mode) & !is.nan(Mode))

mode_labels <- c("Minor", "Major")

mode_plot1 <- ggplot(filtered_data, aes(x = factor(Mode, labels = mode_labels), y = PeakPosition)) +
  geom_boxplot(fill = 'lightblue') +
  labs(title = "Boxplots of Peak Position by Song Modality", x = "Song Modality", y = "Peak Position") +
  scale_y_reverse() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

mode_plot2 <- ggplot(filtered_data, aes(x = factor(Mode, labels = mode_labels), y = WeeksOnChart)) +
  geom_boxplot(fill = 'tan') + 
  labs(title = "Boxplots of Chart Longevity vs Song Modality", x = 'Song Modality', y = 'Longevity (Weeks on Chart)') + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(size = 11), axis.title.x = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))

grid.arrange(othergenre_plot1, othergenre_plot2, mode_plot1, mode_plot2, faco_plot1, faco_plot2, collab_dist, ncol = 2)
```
