analysis.Rmd

---
title: "Rate This Interruption"
output:
  html_document:
    df_print: paged
---


```{r, echo = FALSE, message = FALSE}
#Load Materials
library(googlesheets4)
library(tidyverse)
library(ThemePark)
library(scales)
library(janitor)
library(broom)

couples <- read_sheet(ss = "https://docs.google.com/spreadsheets/d/1R8WJhW0E-5_SnKRsT66lw4UuY45dezcGFVhQgJG2LIA/edit?usp=sharing", 
                      sheet = "Unique_Couples")
couples <- clean_names(couples)
couples <- couples |> 
  filter(!exclude) |> 
  mutate(show = unlist(show)) |> 
  dplyr::select(-starts_with(c("make_it_official", "exclude")), -already_together, -imdb_id)
```

```{r}
couples |> pull(show) |> unique() |> length() #108 shows
couples |> pull(couple) |> unique() |> length() #137 couples
```

```{r}
couples |>
  ggplot(aes(x = mentions)) +
  geom_bar(fill = "purple") +
  theme_barbie() +
  labs(title = "Who are we talking about?",
       x = "Number of Mentions",
       y = "Number of Couples")
```

Who are our power couples? These guys are in the tail of the distribution.

```{r}
couples |> 
  filter(mentions > 10) |> 
  dplyr::select(show, couple, mentions)
```

```{r}
couples |>
  group_by(show) |>
  summarize(year = mean(year_of_start)) |>
  ggplot(aes(x = year)) +
  geom_bar(fill = "#573164") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(color = "#573164"),
        axis.text = element_text(color = "#573164",
                                 size = 24),
        axis.title = element_text(color = "#573164",
                                  size = 26,
                                  face = "bold"),
        plot.title = element_text(color = "#573164", 
                                  size = 32,
                                  face = "bold",
                                  hjust = 0.5)) +
  labs(title = "Distribution of Show Premieres",
       x = "Year of Start",
       y = "Number of Shows") 

ggsave("show_premiere_year.png", 
       path = "/Users/ashleymullan/Documents/Grad School/Wake Forest/M.S. Coursework/Research/TV-Modeling/tv_analysis_its/",
       width = 11.5,
       height = 6,
       unit = "cm")
```

Who are our top 3 oldies and top three new shows?

```{r}
rbind(couples |> 
        group_by(show) |>
        summarize(year = first(year_of_start)) |>
        arrange(year) |> 
        head(3),
      couples |> 
        group_by(show) |>
        summarize(year = first(year_of_start)) |>
        arrange(year) |>
        tail(3))

```


```{r}
couples |>
  group_by(show) |>
  summarize(year = as.character(first(year_of_start))) |>
  mutate(decade = case_when(
    grepl("197", year) ~ "1970s",
    grepl("198", year) ~ "1980s",
    grepl("199", year) ~ "1990s",
    grepl("200", year) ~ "2000s",
    grepl("201", year) ~ "2010s",
    grepl("202", year) ~ "2020s"
  )) |>
  group_by(decade) |>
  summarize(count = n())
```

```{r}
couples[1:20,] |>
  mutate(decade = case_when(
    grepl("197", year_of_start) ~ "1970s",
    grepl("198", year_of_start) ~ "1980s",
    grepl("199", year_of_start) ~ "1990s",
    grepl("200", year_of_start) ~ "2000s",
    grepl("201", year_of_start) ~ "2010s",
    grepl("202", year_of_start) ~ "2020s"
  )) |> 
  mutate(perc_seasons_to_go = (num_seasons - first_kiss_season)/num_seasons) |>
  ggplot(aes(x = decade, y = perc_seasons_to_go)) +
  geom_boxplot(fill = "#573164") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(color = "#573164"),
        axis.text = element_text(color = "#573164",
                                 size = 24),
        axis.title = element_text(color = "#573164",
                                  size = 26,
                                  face = "bold"),
        plot.title = element_text(color = "#573164", 
                                  size = 32,
                                  face = "bold",
                                  hjust = 0.5),
        plot.subtitle = element_text(color = "#573164",
                                     size = 24,
                                     face = "bold",
                                     hjust = 0.5)) +
  scale_y_continuous(labels = scales::percent) +
  labs(title = "Distribution of Kiss Times",
       x = "",
       subtitle = "(Top 20 Couples)",
       y = "Percentage of Show to Go") 

ggsave("show_premiere_decade.png", 
       path = "/Users/ashleymullan/Documents/Grad School/Wake Forest/M.S. Coursework/Research/TV-Modeling/tv_analysis_its/",
       width = 11.5,
       height = 6,
       units = "cm")
```

```{r}
couples |>
  group_by(show) |>
  summarize(pairs = n()) |>
  ggplot(aes(x = pairs)) +
  geom_bar(fill = "purple") +
  theme_barbie() +
  labs(title = "How many couples do we have per show?",
       x = "Number of Couples",
       y = "Number of Shows")
```

Who are our super romantic shows?
```{r}
couples |>
  group_by(show) |>
  summarize(pairs = n()) |>
  filter(pairs > 2)
```

```{r}
couples |>
  #compute percentage of show to go as difference in seasons over total
  mutate(perc_seasons_to_go = (num_seasons - first_kiss_season)/num_seasons) |>
  ggplot(aes(x = perc_seasons_to_go)) +
  #bin at 10%
  geom_histogram(binwidth = 0.10, fill = "purple") +
  scale_x_continuous(labels = scales::percent) +
  coord_flip() +
  labs(x = "Percentage of Show to Go", y = "Number of Couples",
       title = "How early are these kisses happening?",
       subtitle = "Early. Much earlier than I thought.") +
  theme_barbie()
```

Who's a last season kisser?
```{r}
couples |>
  filter(num_seasons == first_kiss_season)
```

```{r}
couples |>
  #compute percentage of show to go as difference in seasons over total
  mutate(perc_seasons_to_go = (num_seasons - first_kiss_season)/num_seasons) |>
  mutate(decade = case_when(
    grepl("197", year_of_start) ~ "1970s",
    grepl("198", year_of_start) ~ "1980s",
    grepl("199", year_of_start) ~ "1990s",
    grepl("200", year_of_start) ~ "2000s",
    grepl("201", year_of_start) ~ "2010s",
    grepl("202", year_of_start) ~ "2020s"
  )) |>
  ggplot(aes(x = perc_seasons_to_go, group = decade)) +
  #bin at 10%
  geom_histogram(binwidth = 0.10, fill = "purple") +
  scale_x_continuous(labels = scales::percent) +
  coord_flip() +
  facet_wrap(vars(decade)) +
  labs(x = "Percentage of Show to Go", y = "Number of Couples",
       title = "How early are these kisses happening?",
       subtitle = "(and does it vary by decade)",
       caption = "(still much earlier than I thought)") +
  theme_barbie()
```

Now, to read in the ratings data.
```{r}
ratings <- read_sheet(ss = "https://docs.google.com/spreadsheets/d/1R8WJhW0E-5_SnKRsT66lw4UuY45dezcGFVhQgJG2LIA/edit?usp=sharing", 
                      sheet = "Episode_Ratings")
ratings <- clean_names(ratings)

# Exclude ratings rows with missing data
ratings <- ratings[complete.cases(ratings), ]
```

Now, to join the two.

```{r, warning = FALSE}
# Left-join episodes and ratings into the couples
ratings_couples <- couples[1:20, ] |> 
  left_join(ratings, 
            by = join_by(show == show))
            

# Create an indicator variable for episodes that aired after first kiss
ratings_couples <- ratings_couples |> 
  group_by(show, couple) |> 
  mutate(show_couple_id = 1:n(), 
         is_first_kiss_ep = season == first_kiss_season & episode == first_kiss_episode,
         first_kiss_id = min(ifelse(test = is_first_kiss_ep, yes = show_couple_id, no = 9999)),
         after_first_kiss = show_couple_id > first_kiss_id) |> 
  dplyr::select(-is_first_kiss_ep) |>
  mutate(eps_since_kiss = show_couple_id - first_kiss_id) |>
  ungroup()
```

```{r}
ratings_couples |> 
  group_by(couple) |>
  summarize(start = first(year_of_start)) |> 
  ggplot(aes(x = start)) +
  geom_bar(fill = "#573164") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(color = "#573164"),
        axis.text = element_text(color = "#573164",
                                 size = 24),
        axis.title = element_text(color = "#573164",
                                  size = 26,
                                  face = "bold"),
        plot.title = element_text(color = "#573164", 
                                  size = 32,
                                  face = "bold",
                                  hjust = 0.5),
        plot.subtitle = element_text(color = "#573164",
                                     size = 28,
                                     face = "bold",
                                     hjust = 0.5)) +
  labs(title = "Distribution of Show Premieres",
       subtitle = "(Top 20 Couples)",
       x = "Year of Start",
       y = "Number of Shows") 

ggsave("show_premiere_year_subset.png", 
       path = "/Users/ashleymullan/Documents/Grad School/Wake Forest/M.S. Coursework/Research/TV-Modeling/tv_analysis_its/",
       width = 11.5,
       height = 6,
       units = "cm")
```

```{r}
#center on first kiss before fitting the model
its_mod <- lm(data = ratings_couples,
              rating ~ eps_since_kiss + after_first_kiss + eps_since_kiss * after_first_kiss)

its_mod |> tidy(conf.int = TRUE) |> 
  dplyr::select(term, estimate, p.value, conf.low, conf.high)

ratings_couples |>
  dplyr::mutate(fit = its_mod$fitted.values,
                cf = its_mod$coefficients[1] + its_mod$coefficients[2] * eps_since_kiss) |>
  ggplot(aes( x = eps_since_kiss, y = rating, group = after_first_kiss)) +
  geom_point(aes(color = after_first_kiss), alpha = 0.5) +
  geom_line(aes(y = fit), color = "#573164", size = 1.2) +
  geom_line(aes(y = cf), linetype = "dotted", 
            color = "#573164", size = 1.2) + ## add dashed line with CF trajectory to plot 
  theme_minimal() +
   theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(color = "#573164"),
        axis.text = element_text(color = "#573164",
                                 size = 24),
        axis.title = element_text(color = "#573164",
                                  size = 26,
                                  face = "bold"),
        plot.title = element_text(color = "#573164", 
                                  size = 32,
                                  face = "bold",
                                  hjust = 0.5),
        plot.subtitle = element_text(color = "#573164",
                                     size = 28,
                                     face = "bold",
                                     hjust = 0.5),
        legend.position = "none") +
  scale_color_manual(values = c("#89c8f5", "#fa9bd1")) +
  labs(x = "Episodes Since the Kiss",
       y = "Rating",
       title = "Rating Trends Before and After the Kiss",
       color = "After The Kiss?")

ggsave("model_all_shows.png", 
       path = "/Users/ashleymullan/Documents/Grad School/Wake Forest/M.S. Coursework/Research/TV-Modeling/tv_analysis_its/",
       width = 11.5,
       height = 6,
       units = "cm")

#delta <- 10
#post <- 0
#ggplot(data |> filter(time > as.Date("2014-06-02") - delta,
                      #time < as.Date("2014-06-02") + delta + post), aes(time, value)) +
  #geom_point(color = "cornflower blue") +
  #geom_ribbon(aes(ymin = m_lb_1, 
                  #ymax = m_ub_1), fill = "red", alpha = 0.25) +
  #geom_line(aes(time, m_pred_1), color = "red", lty = 2) +
  #geom_ribbon(aes(ymin = m_lb, 
        #          ymax = m_ub), fill = "orange", alpha = 0.25) +
  #geom_line(aes(time, m_pred), color = "orange", lty = 1) +
  #geom_vline(xintercept = as.Date("2014-06-02"), lty = 2, color = "grey", lwd = 1) +
  #annotate("label", x = as.Date("2014-06-02"), y = -10, label = "(1)") +
#  scale_x_date(date_labels = "%b %d %Y") +
  #scale_y_continuous(limits = c(min(data$value), max(data$value))) +
#  theme_minimal() + 
#  labs(x = "Date",
 #      y = "Normalized Google Trend Value")
```

**A note about interpreting the model:** 

The ITS model quantifies the altered trajectories of the episode ratings after versus before the will they/won't they couple first kiss. Primary interpretations from the ITS models were two-fold: (1) did the episode ratings change immediately following the first kiss and (2) if so, did they/how quickly did they return to pre-first-kiss levels?

  -   The immediate change in episode ratings after versus before the first kiss, we obtained estimates from the main effect of the post-kiss period.
  
```{r}
its_mod$coefficients["after_first_kissTRUE"] 
```
  
  -   The episode-on-episode change in the episode ratings after versus before the first, we obtained estimates using the interaction effect between the episode-on-episode rate of change and an indicator for the post-kiss period. 
  
```{r}
its_mod$coefficients["eps_since_kiss:after_first_kissTRUE"] 
```
  
Both are presented with their 95% confidence intervals (95% CI). 

```{r}
its_mod_ng <- lm(data = ratings_couples,
              rating ~ eps_since_kiss + after_first_kiss + eps_since_kiss * after_first_kiss, 
              subset = show == "New Girl") ## fit model only to New Girl

its_mod_ng |> tidy(conf.int = TRUE) |>
  dplyr::select(term, estimate, p.value, conf.low, conf.high)

newgirl <- ratings_couples |> 
    filter(show == "New Girl")

newgirl |> 
  dplyr::mutate(fit = its_mod_ng$fitted.values,
                cf = its_mod_ng$coefficients[1] + its_mod_ng$coefficients[2] * eps_since_kiss ## counterfactual, extending pre-kiss model
                ) |>
  ggplot(aes(x = eps_since_kiss, y = rating, group = after_first_kiss)) +
  geom_point(aes(color = after_first_kiss), alpha = 0.5) +
  geom_line(aes(y = fit), size = 1.2, color = "#573164") +
  geom_line(aes(y = cf), linetype = "dotted", 
            size = 1.2, color = "#573164") + ## add dashed line with cf trajectory to plot
  theme_minimal() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(color = "#573164"),
        axis.text = element_text(color = "#573164",
                                 size = 24),
        axis.title = element_text(color = "#573164",
                                  size = 26,
                                  face = "bold"),
        plot.title = element_text(color = "#573164", 
                                  size = 32,
                                  face = "bold",
                                  hjust = 0.5),
        plot.subtitle = element_text(color = "#573164",
                                     size = 28,
                                     face = "bold",
                                     hjust = 0.5),
        legend.position = "none") +
  scale_color_manual(values = c("#89c8f5", "#fa9bd1")) +
  labs(title = "New Girl Rating Trends",
       x = "Episode Since the Kiss",
       y = "Rating",
       color = "After the Kiss?") 

ggsave("model_newgirl.png", 
       path = "/Users/ashleymullan/Documents/Grad School/Wake Forest/M.S. Coursework/Research/TV-Modeling/tv_analysis_its/",
       width = 11.5,
       height = 6,
       units = "cm")
```