Introduction

As social media consultants, our objective was to identify what drives YouTube engagement. For our midterm project, we utilized the datasnaek/youtube-new Kaggle dataset. For our final project, we scraped all the data directly from YouTube every single day using the tuber package, which allowed us to include additional variables like video length, description length, and title punctuation.

Begin Midterm Project Summary

Data Preparation Summary (Midterm Phases 1-3) by Keith

library(ggplot2)
# this was kaggle data we downloaded and would load from local
youtube_model_df <- read.csv("data_set/youtube_model.csv")

cols_to_convert <- c("high_engagement", "high_views", "high_likes", "high_dislikes", "high_comments")

youtube_model_df[cols_to_convert] <- lapply(youtube_model_df[cols_to_convert], function(x) {
  as.numeric(x == "Yes")
})

Model 1: Logistic Regression by Nico

Logistic regression was used to estimate the odds that key video metrics (views, likes, dislikes, comments) drive high engagement, using a binary-encoded version of the cleaned data (youtube_model.csv) that retains duplicate video snapshots to capture engagement decay over time.

engagement_model <- glm(
  high_engagement ~ high_views + high_likes + high_dislikes + high_comments,
  data = youtube_model_df,
  family = "binomial"
)

summary(engagement_model)

## 
## Call:
## glm(formula = high_engagement ~ high_views + high_likes + high_dislikes + 
##     high_comments, family = "binomial", data = youtube_model_df)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.38322    0.01518  -91.11   <2e-16 ***
## high_views    -2.52515    0.06208  -40.67   <2e-16 ***
## high_likes     2.36164    0.05558   42.49   <2e-16 ***
## high_dislikes -0.91792    0.04582  -20.03   <2e-16 ***
## high_comments  1.59291    0.04770   33.39   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 46029  on 40925  degrees of freedom
## Residual deviance: 40006  on 40921  degrees of freedom
## AIC: 40016
## 
## Number of Fisher Scoring iterations: 5

model_summary <- summary(engagement_model)$coefficients[-1, ]


plot_data <- data.frame(
  Metric = rownames(model_summary),
  Impact = model_summary[, "Estimate"]
)

plot_data$Metric <- gsub("high_", "High ", plot_data$Metric)
plot_data$Metric <- gsub("yes", "", plot_data$Metric) 

ggplot(plot_data, aes(x = reorder(Metric, Impact), y = Impact, fill = Impact > 0)) +
  geom_bar(stat = "identity", width = 0.7, color = "black") +
  coord_flip() +  # This flips the chart sideways!

 
  scale_fill_manual(values = c("TRUE" = "#2ecc71", "FALSE" = "#e74c3c"), guide = "none") +
  
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray40", linewidth = 0.8) +
  
  labs(
    title = "What Drives YouTube Engagement?",
    subtitle = "Model Coefficients (Log-Odds Impact on High Engagement)",
    x = "Video Metric",
    y = "Direction & Strength of Impact"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = "bold", size = 16),
    axis.text = element_text(color = "black", face = "bold")
  )

Bonus Model: Engagement Decay Over Time by Nico

Adding days_since_published tests whether engagement odds decay after a video is posted.

engagement_time_model <- glm(
  high_engagement ~ high_views + high_likes + high_dislikes + high_comments + days_since_published,
  data = youtube_model_df,
  family = "binomial"
)

summary(engagement_time_model)

## 
## Call:
## glm(formula = high_engagement ~ high_views + high_likes + high_dislikes + 
##     high_comments + days_since_published, family = "binomial", 
##     data = youtube_model_df)
## 
## Coefficients:
##                       Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.935814   0.021697  -43.13   <2e-16 ***
## high_views           -2.420017   0.063418  -38.16   <2e-16 ***
## high_likes            2.445976   0.056581   43.23   <2e-16 ***
## high_dislikes        -0.913524   0.046647  -19.58   <2e-16 ***
## high_comments         1.641240   0.048677   33.72   <2e-16 ***
## days_since_published -0.081131   0.003093  -26.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 46029  on 40925  degrees of freedom
## Residual deviance: 39084  on 40920  degrees of freedom
## AIC: 39096
## 
## Number of Fisher Scoring iterations: 8

time_range_60 <- seq(1, 60, length.out = 150)

curve_data_60 <- data.frame(
  days_since_published = time_range_60,
  high_views = 1,
  high_likes = 1,
  high_dislikes = 0,
  high_comments = 1
)

curve_data_60$Probability <- predict(engagement_time_model, newdata = curve_data_60, type = "response")

ggplot(curve_data_60, aes(x = days_since_published, y = Probability)) +
  geom_line(color = "#e74c3c", linewidth = 1.5) +  # Decay line
  scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
  
  scale_x_continuous(breaks = c(1, 7, 14, 30, 45, 60), 
                     labels = c("Day 1", "Week 1", "Week 2", "Month 1", "Day 45", "Month 2")) +
  coord_cartesian(xlim = c(1, 60)) + 
  
  labs(
    title = "The 60-Day Engagement Decay Effect",
    subtitle = "Predicted Probability of High Engagement Over a Video's First Two Months",
    x = "Time Since Published",
    y = "Probability of High Engagement Rate"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(face = "bold", size = 16),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(face = "bold")
  )

Model 2: Random Forest By Lennin

A random forest was used to test which video attributes (if any) drive engagement rate. This is panel data (the same video appears across multiple trending days), so videos were deduplicated to one snapshot each, and view counts were log-transformed to manage splits.

library(tidyverse)
library(knitr)

# modelling (random forest sections)
library(lubridate)
library(ranger)
library(caret)
library(vip)
library(pdp)
library(patchwork)

# text mining / word cloud
library(tidytext)
library(widyr)
library(igraph)
library(ggraph)
library(ggwordcloud)

category_map <- c(
  "1"  = "Film & Animation", "2"  = "Autos & Vehicles", "10" = "Music",
  "15" = "Pets & Animals",   "17" = "Sports",            "19" = "Travel & Events",
  "20" = "Gaming",           "22" = "People & Blogs",    "23" = "Comedy",
  "24" = "Entertainment",    "25" = "News & Politics",   "26" = "Howto & Style",
  "27" = "Education",        "28" = "Science & Technology",
  "29" = "Nonprofits & Activism", "43" = "Shows"
)

raw <- read_csv("data_set/USvideos.csv", show_col_types = FALSE)
cat("Raw rows:", nrow(raw), "| Unique videos:", n_distinct(raw$video_id), "\n")

## Raw rows: 40949 | Unique videos: 6351

df <- raw %>%
  group_by(video_id) %>%
  slice_max(views, n = 1, with_ties = FALSE) %>%
  ungroup()

cat("Unique videos after dedup:", nrow(df), "\n")

## Unique videos after dedup: 6351

# feature engineering
df <- df %>%
  mutate(
    engagement_rate    = (likes + dislikes + comment_count) / views,
    title_length       = nchar(title),
    title_word_count   = str_count(title, "\\S+"),
    has_caps           = str_detect(title, "[A-Z]{2,}"),
    has_question       = str_detect(title, "\\?"),
    has_exclaim        = str_detect(title, "!"),
    tag_count          = if_else(tags == "[none]", 0L, str_count(tags, "\\|") + 1L),
    publish_dt         = ymd_hms(publish_time),
    publish_hour       = hour(publish_dt),
    publish_dow        = wday(publish_dt, label = FALSE),
    trend_dt           = as.Date(trending_date, format = "%y.%d.%m"),
    days_to_trend      = as.numeric(trend_dt - as.Date(publish_dt)),
    comments_disabled  = (comments_disabled == "True"),
    ratings_disabled   = (ratings_disabled  == "True"),
    description_length = nchar(coalesce(description, "")),
    log_views          = log1p(views),
    category           = factor(category_map[as.character(category_id)])
  )

df_model <- df %>%
  select(
    engagement_rate,
    title_length, title_word_count, has_caps, has_question, has_exclaim,
    tag_count, publish_hour, publish_dow, days_to_trend,
    comments_disabled, ratings_disabled, description_length,
    log_views, category
  ) %>%
  mutate(
    has_caps          = as.factor(has_caps),
    has_question      = as.factor(has_question),
    has_exclaim       = as.factor(has_exclaim),
    comments_disabled = as.factor(comments_disabled),
    ratings_disabled  = as.factor(ratings_disabled),
    publish_dow       = as.factor(publish_dow)
  ) %>%
  filter(
    engagement_rate >= 0,
    engagement_rate <= quantile(engagement_rate, 0.99, na.rm = TRUE),
    days_to_trend >= 0,
    !is.na(engagement_rate)
  ) %>%
  drop_na()

cat("Final model rows:", nrow(df_model), "\n")

## Final model rows: 6280

summary(df_model$engagement_rate)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.01566 0.02882 0.03404 0.04646 0.13691

# Typical 80/20 split, stratified on the target variable via `caret::createDataPartition` to ensure the engagement rate distribution is balanced across both sets.
set.seed(42)
train_idx <- createDataPartition(df_model$engagement_rate, p = 0.8, list = FALSE)
train <- df_model[ train_idx, ]
test  <- df_model[-train_idx, ]

cat("Train rows:", nrow(train), "| Test rows:", nrow(test), "\n")

## Train rows: 5024 | Test rows: 1256

# Train Random Forest

rf_model <- ranger(
  engagement_rate ~ .,
  data       = train,
  num.trees  = 500,
  importance = "permutation",
  seed       = 42
)

cat("OOB R²:", round(rf_model$r.squared, 4), "\n")

## OOB R²: 0.2665

# Model Evaluation on Test Set

preds  <- predict(rf_model, data = test)$predictions
rmse   <- sqrt(mean((test$engagement_rate - preds)^2))
mae    <- mean(abs(test$engagement_rate - preds))
ss_res <- sum((test$engagement_rate - preds)^2)
ss_tot <- sum((test$engagement_rate - mean(test$engagement_rate))^2)
r2     <- 1 - ss_res / ss_tot

cat("Test R²:", round(r2, 4), "| RMSE:", round(rmse, 5), "| MAE:", round(mae, 5), "\n")

## Test R²: 0.2624 | RMSE: 0.02088 | MAE: 0.01532

Random Forest Data Viz

tibble(actual = test$engagement_rate, predicted = preds) %>%
  ggplot(aes(x = actual, y = predicted)) +
  geom_point(alpha = 0.3, color = "#2196F3", size = 1.2) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +
  labs(
    title    = "Actual vs. Predicted Engagement Rate",
    subtitle = paste0("Test R² = ", round(r2, 4)),
    x = "Actual", y = "Predicted"
  ) +
  theme_minimal(base_size = 13)

vip(rf_model, num_features = 15, aesthetics = list(fill = "#2196F3")) +
  labs(
    title    = "What Drives YouTube Engagement Rate?",
    subtitle = "Random Forest — permutation importance (US trending videos)",
    x = "Feature", y = "Importance"
  ) +
  theme_minimal(base_size = 13)

importance_df <- tibble(
  feature    = names(rf_model$variable.importance),
  importance = rf_model$variable.importance
) %>% arrange(desc(importance))

knitr::kable(importance_df, digits = 6, caption = "Feature importance (permutation)")

Feature importance (permutation)
feature	importance
title_length	0.000135
category	0.000113
title_word_count	0.000109
description_length	0.000062
publish_hour	0.000059
log_views	0.000051
tag_count	0.000041
days_to_trend	0.000031
has_caps	0.000019
publish_dow	0.000008
has_exclaim	0.000008
has_question	0.000003
comments_disabled	0.000000
ratings_disabled	0.000000

pdp_plot <- function(feature, label) {
  partial(rf_model, pred.var = feature, train = train, plot = FALSE) %>%
    as_tibble() %>%
    rename(x = 1, yhat = 2) %>%
    ggplot(aes(x = x, y = yhat)) +
    geom_line(color = "#2196F3", linewidth = 1.1) +
    geom_rug(data = train, aes_string(x = feature, y = NULL),
             alpha = 0.15, sides = "b") +
    labs(title = label, x = feature, y = "Predicted engagement rate") +
    theme_minimal(base_size = 12)
}

p_title_len  <- pdp_plot("title_length",       "Title Length")
p_word_count <- pdp_plot("title_word_count",    "Title Word Count")
p_desc       <- pdp_plot("description_length",  "Description Length")
p_hour       <- pdp_plot("publish_hour",        "Publish Hour")
p_tags       <- pdp_plot("tag_count",           "Tag Count")
p_views      <- pdp_plot("log_views",           "Log Views (channel size proxy)")

(p_title_len | p_word_count | p_desc) /
(p_hour      | p_tags       | p_views) +
  plot_annotation(
    title    = "How Each Feature Shapes Engagement Rate",
    subtitle = "Partial dependence plots. all other features held at their mean",
    theme    = theme_minimal(base_size = 13)
  )

df_model %>%
  mutate(category = fct_reorder(category, engagement_rate, .fun = median)) %>%
  ggplot(aes(x = category, y = engagement_rate)) +
  geom_boxplot(fill = "#2196F3", alpha = 0.6, outlier.alpha = 0.2) +
  coord_flip() +
  labs(
    title    = "Engagement Rate by Category",
    subtitle = "Median-sorted",
    x = NULL, y = "Engagement Rate"
  ) +
  theme_minimal(base_size = 12)

df_model %>%
  mutate(publish_dow = as.integer(as.character(publish_dow))) %>%
  group_by(publish_hour, publish_dow) %>%
  summarise(avg_engagement = mean(engagement_rate), .groups = "drop") %>%
  mutate(
    day_label = factor(publish_dow, levels = 1:7,
                       labels = c("Sun","Mon","Tue","Wed","Thu","Fri","Sat"))
  ) %>%
  ggplot(aes(x = publish_hour, y = day_label, fill = avg_engagement)) +
  geom_tile(color = "white") +
  scale_fill_gradient(low = "#e3f2fd", high = "#1565C0", name = "Avg engagement") +
  labs(
    title    = "Best Time to Post",
    subtitle = "Average engagement rate by publish hour × day of week",
    x = "Hour of day (UTC)", y = NULL
  ) +
  theme_minimal(base_size = 12)

Begin Final Project

Data Scraping by Keith

We setup a daily scraper to pull the top 200 trending videos in the US, along with their comments. The data is stored in CSV files and at the end, we combine all the daily snapshots into a single dataset for analysis. The scraper uses the YouTube API and is designed to run automatically every day.

Code Shown for reference only. This chunk is not executed. It ran on a schedule from 2026-07-06 to 2026-08-01 to collect the data.

## Daily scraper: top 200 "most popular" YouTube videos in the US + their comments

library(httr)
library(jsonlite)
library(dplyr)
library(readr)
library(purrr)
library(tibble)

## ---- Config -----------------------------------------------------------

readRenviron(".env")
api_key <- Sys.getenv("youtube_api_key")
if (nchar(api_key) == 0) stop("API key not found. Check your .env file (youtube_api_key=...).")

REGION_CODE          <- "US"
TARGET_VIDEO_COUNT   <- 200
COMMENTS_PER_VIDEO   <- 20     # keep modest to conserve daily quota
REQUEST_DELAY_SEC    <- 0.5    # politeness delay between any two API calls
MAX_RETRIES          <- 5

today_str  <- format(Sys.Date(), "%Y-%m-%d")
out_dir    <- file.path("final/scrapers/data", today_str)
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)

videos_url   <- "https://www.googleapis.com/youtube/v3/videos"
comments_url <- "https://www.googleapis.com/youtube/v3/commentThreads"

## ---- Helper: polite GET with exponential backoff ----------------------

polite_get <- function(url, query) {
  for (attempt in seq_len(MAX_RETRIES)) {
    resp <- GET(url = url, query = query)
    status <- status_code(resp)

    if (status == 200) {
      Sys.sleep(REQUEST_DELAY_SEC)
      return(content(resp, as = "text", encoding = "UTF-8") |> fromJSON(flatten = TRUE))
    }

    parsed <- tryCatch(
      content(resp, as = "text", encoding = "UTF-8") |> fromJSON(flatten = TRUE),
      error = function(e) NULL
    )
    reason <- parsed$error$errors$reason[1]

    # Comments disabled / not found — don't retry, just signal caller to skip
    if (status == 403 && !is.na(reason) && reason %in% c("commentsDisabled", "forbidden")) {
      return(list(error = list(message = reason, skip = TRUE)))
    }

    # Quota exceeded — stop the whole run, retrying won't help today
    if (status == 403 && !is.na(reason) && reason == "quotaExceeded") {
      stop("Daily YouTube API quota exceeded. Stopping run.")
    }

    # Rate limited or transient server error — back off and retry
    if (status %in% c(429, 500, 502, 503)) {
      wait <- REQUEST_DELAY_SEC * (2 ^ attempt)
      message("Got status ", status, " — backing off for ", round(wait, 1), "s (attempt ", attempt, "/", MAX_RETRIES, ")")
      Sys.sleep(wait)
      next
    }

    # Any other error — surface it
    warning("Unexpected status ", status, ": ", parsed$error$message %||% "unknown error")
    Sys.sleep(REQUEST_DELAY_SEC)
    return(list(error = list(message = parsed$error$message %||% "unknown error")))
  }
  list(error = list(message = paste("Max retries exceeded for", url)))
}

`%||%` <- function(a, b) if (is.null(a) || is.na(a)) b else a

## ---- Step 1: fetch top 200 most-popular US videos ---------------------

fetch_top_videos <- function() {
  all_videos <- list()
  page_token <- NULL

  repeat {
    query <- list(
      part        = "snippet,statistics,contentDetails",
      chart       = "mostPopular",
      regionCode  = REGION_CODE,
      maxResults  = 50,
      key         = api_key
    )
    if (!is.null(page_token)) query$pageToken <- page_token

    page <- polite_get(videos_url, query)

    if (!is.null(page$error)) {
      warning("Error fetching video list: ", page$error$message)
      break
    }

    items <- page$items
    if (is.null(items) || NROW(items) == 0) break

    batch <- tibble(
      video_id      = items$id,
      title         = items$snippet.title,
      channel_title = items$snippet.channelTitle,
      published_at  = items$snippet.publishedAt,
      category_id   = items$snippet.categoryId,
      view_count    = as.numeric(items$statistics.viewCount),
      like_count    = as.numeric(items$statistics.likeCount),
      comment_count = as.numeric(items$statistics.commentCount),
      duration      = items$contentDetails.duration,
      scrape_date   = today_str
    )

    all_videos <- append(all_videos, list(batch))
    message("Fetched ", sum(map_int(all_videos, nrow)), " videos so far...")

    page_token <- page$nextPageToken
    total_so_far <- sum(map_int(all_videos, nrow))
    if (is.null(page_token) || total_so_far >= TARGET_VIDEO_COUNT) break
  }

  bind_rows(all_videos) |> slice_head(n = TARGET_VIDEO_COUNT)
}

top_videos <- fetch_top_videos()
message("\nTotal videos collected: ", nrow(top_videos))

videos_out_path <- file.path(out_dir, "top200_videos.csv")
write_csv(top_videos, videos_out_path)
message("Saved video list to ", videos_out_path)

## ---- Step 2: fetch comments for each video -----------------------------

fetch_comments_for_video <- function(video_id) {
  query <- list(
    part       = "snippet",
    videoId    = video_id,
    maxResults = COMMENTS_PER_VIDEO,
    textFormat = "plainText",
    order      = "relevance",
    key        = api_key
  )

  page <- polite_get(comments_url, query)

  if (!is.null(page$error)) {
    if (isTRUE(page$error$skip)) {
      message("  Comments disabled/unavailable for ", video_id, " — skipping.")
    } else {
      warning("  Error fetching comments for ", video_id, ": ", page$error$message)
    }
    return(NULL)
  }

  items <- page$items
  if (is.null(items) || NROW(items) == 0) return(NULL)

  tibble(
    video_id     = video_id,
    comment_id   = items$id,
    author       = items$snippet.topLevelComment.snippet.authorDisplayName,
    text         = items$snippet.topLevelComment.snippet.textOriginal,
    like_count   = as.integer(items$snippet.topLevelComment.snippet.likeCount),
    published_at = items$snippet.topLevelComment.snippet.publishedAt,
    scrape_date  = today_str
  )
}

message("\nFetching up to ", COMMENTS_PER_VIDEO, " comments for each of ", nrow(top_videos), " videos...")

all_comments <- map(seq_len(nrow(top_videos)), function(i) {
  vid <- top_videos$video_id[i]
  message("[", i, "/", nrow(top_videos), "] ", vid)
  fetch_comments_for_video(vid)
})

comments_df <- bind_rows(all_comments)
message("\nTotal comments collected: ", nrow(comments_df))

comments_out_path <- file.path(out_dir, "top200_comments.csv")
write_csv(comments_df, comments_out_path)
message("Saved comments to ", comments_out_path)

message("\nDone. Run for ", today_str, " complete.")

Word Cloud by Nico

# Libraries already loaded above

# 1. Define GitHub URLs
comments_url <- "https://raw.githubusercontent.com/ok-jaime/msba_580_final_project/54d4d92874b958593b09c1aab3403d1f588d7a75/all_comments.csv"
videos_url <- "https://raw.githubusercontent.com/ok-jaime/msba_580_final_project/54d4d92874b958593b09c1aab3403d1f588d7a75/all_videos.csv"

# 2. Read Data
video_data <- read_csv(videos_url, show_col_types = FALSE)
comment_data <- read_csv(comments_url, show_col_types = FALSE)

# 3. Preprocess and Clean Text
data(stop_words)

custom_stop_words <- tibble(
  word = c(
    # Bare Minimum Web & YouTube Boilerplate
    "https", "http", "www", "com"
  ), 
  lexicon = "custom"
)

all_stop_words <- bind_rows(stop_words, custom_stop_words)

youtube_words <- comment_data %>%
  mutate(id = row_number()) %>%
  select(id, text) %>%                
  filter(!is.na(text)) %>%             
  unnest_tokens(word, text) %>%
  anti_join(all_stop_words, by = "word") %>%
  filter(str_detect(word, "^[a-z]+$")) %>% 
  filter(str_length(word) > 2)

# Filter out low-frequency words (crucial to prevent pairwise_cor from crashing)
frequent_words <- youtube_words %>%
  count(word) %>%
  filter(n >= 70) %>%
  pull(word)

youtube_words_filtered <- youtube_words %>%
  filter(word %in% frequent_words)

# --- WORD CLOUD VISUALIZATION ---

# 1. Calculate word frequencies from cleaned tokens
youtube_words_top100 <- youtube_words_filtered %>%
  count(word, sort = TRUE) %>%
  slice_max(n, n = 100) # Take top 100 most frequent words

# 2. Calculate word frequencies from cleaned tokens
youtube_words_raw_top50 <- youtube_words %>%
  count(word, sort = TRUE) %>%
  slice_max(n, n = 50) # Take top 100 most frequent words

# 2. Render Word Cloud
ggplot(youtube_words_top100, aes(label = word, size = n, color = n)) +
  geom_text_wordcloud_area(shape = "circle") +
  scale_size_area(max_size = 15) + # Adjust size scaling
  scale_color_gradient(low = "darkgray", high = "firebrick") +
  theme_minimal() +
  labs(
    title = "Top 100 Most Frequent Words in YouTube Comments",
  )

Word Co-occurrence

# --- STEP 4: RAW PAIRWISE COUNT ---
youtube_pairs <- youtube_words_filtered %>%
  pairwise_count(word, id, sort = TRUE, upper = FALSE)

knitr::kable(head(youtube_pairs, 20),
             caption = "Top 20 word pairs appearing in the same comment")

Top 20 word pairs appearing in the same comment
item1	item2	n
play	game	268
love	videos	266
playing	game	177
love	video	168
love	game	167
time	game	143
people	game	133
video	game	128
played	game	119
happy	birthday	119
hour	gang	116
fun	game	108
video	watch	107
love	time	105
love	song	103
game	games	100
video	watching	99
love	watching	96
video	time	94
love	play	91

# --- STEP 5: CO-OCCURRENCE SPARSE MATRIX ---
# Mirror the pairs so the matrix is symmetric, then cast to sparse form.
youtube_sparse <- youtube_pairs %>%
  bind_rows(youtube_pairs %>% rename(item1 = item2, item2 = item1)) %>%
  cast_sparse(item1, item2, n)

dim(youtube_sparse)

## [1] 608 608

# Inspect matrix subset
youtube_sparse[1:10, 1:10]

## 10 x 10 sparse Matrix of class "dgCMatrix"

##                                              
## play    268  29  51   2   5  56  84   8 80 50
## love    167 266 168  13   4  59 105 103 60 96
## playing 177  18  58   .   1  17  47   9 66 27
## time    143  59  94   5   9  38   .  28 38 87
## people  133  24  39   2   2  32  65  10 55 28
## video   128  58   .   9   6 107  94  32 39 99
## played  119   7  11   .   1   7  38  12 23  9
## happy    43  23  29 119   1  16  29   6 10 20
## hour     15   9  40   1 116   9   9   4  1  7
## fun     108  32  36   1   1  68  39   9 35 23

# --- STEP 6: NETWORK GRAPH (using pairwise counts) ---
# The "fr" layout is stochastic, so seed it to keep the figure stable per knit.
set.seed(42)

youtube_pairs %>%
  filter(n > 5) %>%                            # drop pairs seen 5 or fewer times
  slice_max(n, n = 50, with_ties = FALSE) %>%   # keep the 50 strongest pairs
  graph_from_data_frame() %>%
  ggraph(layout = "fr") +
  geom_edge_link(aes(edge_alpha = n), show.legend = FALSE) +
  geom_node_point(color = "firebrick", size = 4) +
  geom_node_text(aes(label = name), repel = TRUE) +
  theme_void() +
  labs(title = "Word Co-occurrence Network — 50 Strongest Pairs")

Two structures show up. A dense core of gaming and viewing vocabulary: game, play, video, watch, love, time, reflecting that Gaming is roughly 60% of this scrape.

Random Forest on New 2026 Scraped Data

Same methodology as the midterm. dedupe to one snapshot per video, log-transformed views, 80/20 split, 500 trees, permutation importance, seed 42. The feature set had to change because the YouTube API no longer returns what Kaggle had in 2018.

Lost: dislikes (YouTube removed them in 2021, so the target is now (likes + comments) / views), tag_count and description_length (our scraper didn’t capture them).
Gained: duration_seconds — video length (we didn’t have this in midterm)
Dropped on purpose: days_to_trend. The midterm tracked videos over months; this scrape ran 19 days, so this would be heavily skewed.

videos_url <- "https://raw.githubusercontent.com/ok-jaime/msba_580_final_project/54d4d92874b958593b09c1aab3403d1f588d7a75/all_videos.csv"

v26 <- read_csv(videos_url, col_types = cols(
  video_id = col_character(), title = col_character(),
  channel_title = col_character(), published_at = col_character(),
  category_id = col_character(), view_count = col_double(),
  like_count = col_double(), comment_count = col_double(),
  duration = col_character(), scrape_date = col_character(),
  composite_key = col_character()
))

df26 <- v26 %>%
  filter(!is.na(view_count), view_count > 0) %>%
  mutate(
    like_count    = coalesce(like_count, 0),
    comment_count = coalesce(comment_count, 0),
    # 53 video_ids starting with "-" were corrupted to "#NAME?" by a spreadsheet
    # round-trip; fall back to a title+channel key so those rows still dedup.
    vid_key = if_else(video_id == "#NAME?",
                      paste(title, channel_title, published_at), video_id)
  ) %>%
  group_by(vid_key) %>%
  slice_max(view_count, n = 1, with_ties = FALSE) %>%
  ungroup() %>%
  mutate(
    engagement_rate  = (like_count + comment_count) / view_count,
    title_length     = nchar(title),
    title_word_count = str_count(title, "\\S+"),
    has_caps         = str_detect(title, "[A-Z]{2,}"),
    has_question     = str_detect(title, "\\?"),
    has_exclaim      = str_detect(title, "!"),
    publish_dt       = ymd_hms(published_at, tz = "UTC"),
    publish_hour     = hour(publish_dt),
    publish_dow      = wday(publish_dt, label = FALSE),
    log_views        = log1p(view_count),
    duration_seconds = suppressWarnings(as.numeric(lubridate::duration(duration))),
    days_to_trend    = as.numeric(mdy(scrape_date) - as.Date(publish_dt)),
    category = fct_lump_min(factor(category_map[category_id]), 20,
                            other_level = "Other")
  )

SHARED26   <- c("title_length", "title_word_count", "has_caps", "has_question",
                "has_exclaim", "publish_hour", "publish_dow", "log_views",
                "category")
FEATURES26 <- c(SHARED26, "duration_seconds")

dm26 <- df26 %>%
  select(engagement_rate, days_to_trend, all_of(FEATURES26)) %>%
  mutate(across(c(has_caps, has_question, has_exclaim, publish_dow), as.factor)) %>%
  filter(
    engagement_rate >= 0,
    engagement_rate <= quantile(engagement_rate, 0.99, na.rm = TRUE),
    days_to_trend >= 0,
    !is.na(engagement_rate)
  ) %>%
  select(-days_to_trend) %>%
  drop_na()

cat("Videos:", nrow(dm26), "| Scrape days:", n_distinct(v26$scrape_date), "\n")

## Videos: 2818 | Scrape days: 19

set.seed(42)
idx26    <- createDataPartition(dm26$engagement_rate, p = 0.8, list = FALSE)
train26  <- dm26[idx26, ]
test26   <- dm26[-idx26, ]

rf26 <- ranger(engagement_rate ~ ., data = train26, num.trees = 500,
               importance = "permutation", seed = 42)

pred26 <- predict(rf26, data = test26)$predictions
r2_26  <- 1 - sum((test26$engagement_rate - pred26)^2) /
              sum((test26$engagement_rate - mean(test26$engagement_rate))^2)

cat("OOB R²:", round(rf26$r.squared, 4),
    "| Test R²:", round(r2_26, 4),
    "| RMSE:", round(sqrt(mean((test26$engagement_rate - pred26)^2)), 5), "\n")

## OOB R²: 0.2098 | Test R²: 0.2527 | RMSE: 0.02632

Results

vip(rf26, num_features = 10, aesthetics = list(fill = "#2196F3")) +
  labs(title = "What Drives Engagement Rate? (2026 scraped data)",
       subtitle = "Random Forest — permutation importance",
       x = "Feature", y = "Importance") +
  theme_minimal(base_size = 13)

imp26 <- tibble(feature = names(rf26$variable.importance),
                importance = rf26$variable.importance) %>%
  arrange(desc(importance)) %>%
  mutate(share = importance / sum(importance))

knitr::kable(imp26, digits = 5, caption = "Feature importance (2026 data)")

Feature importance (2026 data)
feature	importance	share
duration_seconds	0.00026	0.34126
log_views	0.00013	0.16851
category	0.00010	0.12548
title_length	0.00010	0.12538
title_word_count	0.00006	0.07608
publish_hour	0.00006	0.07491
has_caps	0.00005	0.06710
has_exclaim	0.00001	0.01446
publish_dow	0.00000	0.00555
has_question	0.00000	0.00126

Video length is the top driver: a variable the midterm couldn’t test at all. Channel size (log_views) and category follow. Title features still matter but rank lower than they did in the midterm.

Are these rankings real? (noise-floor check)

Every feature gets a shuffled copy added to the model. A shuffled copy keeps its exact distribution but has its link to engagement destroyed, so it’s meaningless by construction — whatever score the best shadow reaches is what a feature can earn by pure chance. Repeated over 10 train/test splits.

stab26 <- map_dfr(1:10, function(s) {
  set.seed(s)
  shadows <- dm26 %>% select(all_of(FEATURES26)) %>%
    mutate(across(everything(), ~ sample(.x)))
  names(shadows) <- paste0("shadow_", names(shadows))
  aug <- bind_cols(dm26, shadows)

  i  <- createDataPartition(aug$engagement_rate, p = 0.8, list = FALSE)
  rf <- ranger(engagement_rate ~ ., data = aug[i, ], num.trees = 500,
               importance = "permutation", seed = s)
  im <- rf$variable.importance
  real <- !str_starts(names(im), "shadow_")
  tibble(feature = names(im), imp = im, is_real = real,
         share = 100 * im / sum(im[real]))
})

NOISE_FLOOR <- max(stab26$share[!stab26$is_real])

stab_sum <- stab26 %>% filter(is_real) %>%
  group_by(feature) %>%
  summarise(mean_share = mean(share), lo = min(share), hi = max(share),
            .groups = "drop") %>%
  mutate(verdict = if_else(mean_share > NOISE_FLOOR,
                           "real", "indistinguishable from noise")) %>%
  arrange(desc(mean_share))

knitr::kable(stab_sum, digits = 3,
             caption = paste0("Noise floor = ", round(NOISE_FLOOR, 2),
                              "% of total importance"))

Noise floor = 1.87% of total importance
feature	mean_share	lo	hi	verdict
duration_seconds	36.504	34.856	39.690	real
log_views	22.763	18.753	27.374	real
title_length	11.659	9.641	13.561	real
category	8.283	7.088	10.304	real
publish_hour	7.645	5.665	9.372	real
title_word_count	6.430	4.776	8.011	real
has_caps	4.732	2.740	6.213	real
publish_dow	1.013	0.416	2.232	indistinguishable from noise
has_exclaim	0.884	0.333	1.634	indistinguishable from noise
has_question	0.087	-0.129	0.234	indistinguishable from noise

stab_sum %>%
  mutate(feature = fct_reorder(feature, mean_share)) %>%
  ggplot(aes(x = feature, y = mean_share)) +
  geom_hline(yintercept = NOISE_FLOOR, linetype = "dashed", color = "#d32f2f") +
  geom_linerange(aes(ymin = lo, ymax = hi), color = "#90A4AE", linewidth = 1) +
  geom_point(size = 2.6, color = "#2196F3") +
  coord_flip() +
  labs(title = "Which drivers survive a noise check?",
       subtitle = "Mean importance over 10 splits (bars = min-max). Below the dashed line = no better than a shuffled feature.",
       x = NULL, y = "Share of total importance (%)") +
  theme_minimal(base_size = 12)

Only the features above the dashed line are worth interpreting. Notably title punctuation (has_question, has_exclaim) falls at or below the noise floor. This means the punctuation angle from our intro doesn’t survive scrutiny on this data.

Testing the Old Model Against the New Data!!

Comparing importance rankings only tells us whether the two eras agree on which features matter. The stronger test: train on the 2017-18 Kaggle data and see if it can actually predict 2026 videos. Both models are scored on the same held-out sets, using the 9 features and the same dislike-free target, so the only thing that changes is which era the model learned from.

if (!exists("raw")) raw <- read_csv("data_set/USvideos.csv", show_col_types = FALSE)

dm_old <- raw %>%
  group_by(video_id) %>% slice_max(views, n = 1, with_ties = FALSE) %>% ungroup() %>%
  mutate(
    engagement_rate  = (likes + comment_count) / views,   # harmonised: no dislikes
    title_length     = nchar(title),
    title_word_count = str_count(title, "\\S+"),
    has_caps         = str_detect(title, "[A-Z]{2,}"),
    has_question     = str_detect(title, "\\?"),
    has_exclaim      = str_detect(title, "!"),
    publish_dt       = ymd_hms(publish_time),
    publish_hour     = hour(publish_dt),
    publish_dow      = wday(publish_dt, label = FALSE),
    days_to_trend    = as.numeric(as.Date(trending_date, format = "%y.%d.%m") -
                                  as.Date(publish_dt)),
    log_views        = log1p(views),
    category = fct_lump_min(factor(category_map[as.character(category_id)]), 20,
                            other_level = "Other")
  ) %>%
  select(engagement_rate, days_to_trend, all_of(SHARED26)) %>%
  mutate(across(c(has_caps, has_question, has_exclaim, publish_dow), as.factor)) %>%
  filter(engagement_rate >= 0,
         engagement_rate <= quantile(engagement_rate, 0.99, na.rm = TRUE),
         days_to_trend >= 0, !is.na(engagement_rate)) %>%
  select(-days_to_trend) %>%
  drop_na()

# One model must score the other era's rows, so factor levels must match exactly.
align_eras <- function(a, b) {
  common <- intersect(names(a), names(b))
  a <- a[, common, drop = FALSE]; b <- b[, common, drop = FALSE]
  for (nm in common) {
    if (is.factor(a[[nm]]) || is.factor(b[[nm]])) {
      la <- unique(as.character(a[[nm]])); lb <- unique(as.character(b[[nm]]))
      shared <- sort(intersect(la, lb))
      lev <- if (length(setdiff(la, shared)) > 0 || length(setdiff(lb, shared)) > 0) {
        union(shared, "Other")   # union, not c() — "Other" may already exist
      } else shared
      to_f <- function(v) factor(ifelse(as.character(v) %in% shared,
                                        as.character(v), "Other"), levels = lev)
      a[[nm]] <- to_f(a[[nm]]); b[[nm]] <- to_f(b[[nm]])
    }
  }
  list(old = a, new = b, common = setdiff(common, "engagement_rate"))
}

al <- align_eras(dm_old, dm26)

split_era <- function(d) {
  set.seed(42)
  i <- createDataPartition(d$engagement_rate, p = 0.8, list = FALSE)
  list(train = d[i, ], test = d[-i, ])
}
s_old <- split_era(al$old)
s_new <- split_era(al$new)

m_old <- ranger(engagement_rate ~ ., data = s_old$train, num.trees = 500, seed = 42)
m_new <- ranger(engagement_rate ~ ., data = s_new$train, num.trees = 500, seed = 42)

eval_on <- function(m, nd, trained_on, eval_set) {
  p <- predict(m, data = nd)$predictions; a <- nd$engagement_rate
  tibble(trained_on, eval_set,
         r2_asis = 1 - sum((a - p)^2) / sum((a - mean(a))^2),
         pearson_r = cor(p, a),
         r2_shape = cor(p, a)^2,       # best R² any rescaling could reach
         bias = mean(p) - mean(a))
}

transfer <- bind_rows(
  eval_on(m_old, s_old$test, "Kaggle 2017-18", "Kaggle test"),
  eval_on(m_new, s_old$test, "Scraped 2026",   "Kaggle test"),
  eval_on(m_old, s_new$test, "Kaggle 2017-18", "2026 test"),
  eval_on(m_new, s_new$test, "Scraped 2026",   "2026 test")
)

knitr::kable(transfer, digits = 4, caption = "Cross-era transfer")

Cross-era transfer
trained_on	eval_set	r2_asis	pearson_r	r2_shape	bias
Kaggle 2017-18	Kaggle test	0.1679	0.4104	0.1685	0.0003
Scraped 2026	Kaggle test	-0.6415	-0.0290	0.0008	0.0133
Kaggle 2017-18	2026 test	-0.1138	-0.0260	0.0007	-0.0041
Scraped 2026	2026 test	0.1578	0.3985	0.1588	0.0002

Results

as_num <- function(v) if (is.numeric(v)) as.numeric(v) else
  if (is.factor(v) && nlevels(v) == 2) as.numeric(v) - 1 else NULL

signs <- map_dfr(al$common, function(f) {
  a <- as_num(al$old[[f]]); b <- as_num(al$new[[f]])
  if (is.null(a) || is.null(b)) return(NULL)
  ra <- cor(a, al$old$engagement_rate, method = "spearman")
  rb <- cor(b, al$new$engagement_rate, method = "spearman")
  tibble(feature = f, rho_kaggle = ra, rho_2026 = rb,
         direction = if_else(sign(ra) == sign(rb), "same", "REVERSED"),
         pct_2026_outside_kaggle_range =
           100 * mean(b < quantile(a, .01) | b > quantile(a, .99)))
}) %>% arrange(direction, desc(abs(rho_kaggle)))

knitr::kable(signs, digits = 3,
             caption = "Direction of each relationship, by era (Spearman)")

Direction of each relationship, by era (Spearman)
feature	rho_kaggle	rho_2026	direction	pct_2026_outside_kaggle_range
title_length	-0.233	0.059	REVERSED	3.903
title_word_count	-0.199	0.014	REVERSED	1.845
has_exclaim	0.107	-0.047	REVERSED	0.000
log_views	0.083	-0.254	REVERSED	0.071
publish_hour	0.163	0.108	same	0.000
has_question	0.070	0.053	same	0.000
has_caps	-0.014	-0.081	same	0.000

signs %>%
  pivot_longer(c(rho_kaggle, rho_2026), names_to = "era", values_to = "rho") %>%
  mutate(era = if_else(era == "rho_kaggle", "Kaggle 2017-18", "Scraped 2026"),
         feature = fct_reorder(feature, abs(rho), .fun = max)) %>%
  ggplot(aes(x = rho, y = feature, color = era)) +
  geom_vline(xintercept = 0, color = "grey40") +
  geom_line(aes(group = feature), color = "grey70", linewidth = 0.8) +
  geom_point(size = 3) +
  scale_color_manual(values = c("Kaggle 2017-18" = "#90A4AE",
                                "Scraped 2026" = "#EF6C00"), name = NULL) +
  labs(title = "The relationships reversed, they didn't just weaken",
       subtitle = "Spearman correlation with engagement rate. Points straddling zero flipped direction.",
       x = "Correlation with engagement rate", y = NULL) +
  theme_minimal(base_size = 12)

The old model does not transfer. Trained on 2017-18 and pointed at 2026 videos it scores a negative R² (worse than just guessing the average) and its predictions correlate with reality at roughly zero. The reverse direction fails identically.

Why: the relationships reversed sign. Four of seven testable relationships flipped between eras, most dramatically log_views (bigger videos used to have slightly higher engagement rates; now they have clearly lower ones) and title_length (short titles used to help; now they slightly hurt). At most ~4% of 2026 videos fall outside the Kaggle range, so it saw plenty of comparable videos and simply learned the opposite lesson.

The takeaway, and the reason to run this test. Importance rankings across the two eras correlate at ρ ≈ 0.73, which looks like a successful replication. It isn’t. Permutation importance measures magnitude, not direction, so it cannot see a sign flip: both eras agree title_length is a top-3 driver precisely because it mattered a lot in both, just in opposite directions. Our midterm’s variable selection replicated; its advice inverted.

What Drives YouTube Engagement?

Lennin, Keith, Nico

2026-08-04

Introduction

Begin Midterm Project Summary

Data Preparation Summary (Midterm Phases 1-3) by Keith

Model 1: Logistic Regression by Nico

Bonus Model: Engagement Decay Over Time by Nico

Model 2: Random Forest By Lennin

Random Forest Data Viz

Begin Final Project

Data Scraping by Keith

Word Cloud by Nico

Word Co-occurrence

Random Forest on New 2026 Scraped Data

Results

Are these rankings real? (noise-floor check)

Testing the Old Model Against the New Data!!

Results