What Drives YouTube Engagement?
Lennin, Keith, Nico
2026-08-04
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: 5model_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: 8time_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: 6351df <- 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: 6280summary(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.01532Random 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 |
|---|---|
| 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
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",
)
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_countanddescription_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: 19set.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.02632Results
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 | 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"))| 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")| 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)")| 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.