What Drives YouTube Engagement?
Lennin, Keith, Nico
2026-07-06
Introduction
As social media consultants, our team’s objective was to see what type of YouTube content drives the highest engagement. This R markdown will show and covers Phase 1 & 2 (Exploratory Data Analysis) and Phase 3 (Data Validation) using the datasnaek YouTube trending dataset sourced from Kaggle, which contains real trending video statistics from the United States.
Data Source
The dataset (USvideos.csv) is sourced from Mitchell JY’s
Trending YouTube Scraper via Kaggle
(datasnaek/youtube-new). It contains 40,949 observations
across 16 variables including views, likes, dislikes, comment count,
publish time, and trending date.
Phase 1 & 2: Loading Data and Understanding the Shape
library(tidyverse)
youtube <- read_csv("data_set/USvideos.csv")
glimpse(youtube)## Rows: 40,949
## Columns: 16
## $ video_id <chr> "2kyS6SvSYSE", "1ZAPwfrtAFY", "5qpjK5DgCt4", "p…
## $ trending_date <chr> "17.14.11", "17.14.11", "17.14.11", "17.14.11",…
## $ title <chr> "WE WANT TO TALK ABOUT OUR MARRIAGE", "The Trum…
## $ channel_title <chr> "CaseyNeistat", "LastWeekTonight", "Rudy Mancus…
## $ category_id <dbl> 22, 24, 23, 24, 24, 28, 24, 28, 1, 25, 17, 24, …
## $ publish_time <dttm> 2017-11-13 17:13:01, 2017-11-13 07:30:00, 2017…
## $ tags <chr> "SHANtell martin", "last week tonight trump pre…
## $ views <dbl> 748374, 2418783, 3191434, 343168, 2095731, 1191…
## $ likes <dbl> 57527, 97185, 146033, 10172, 132235, 9763, 1599…
## $ dislikes <dbl> 2966, 6146, 5339, 666, 1989, 511, 2445, 778, 11…
## $ comment_count <dbl> 15954, 12703, 8181, 2146, 17518, 1434, 1970, 34…
## $ thumbnail_link <chr> "https://i.ytimg.com/vi/2kyS6SvSYSE/default.jpg…
## $ comments_disabled <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ ratings_disabled <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ video_error_or_removed <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ description <chr> "SHANTELL'S CHANNEL - https://www.youtube.com/s…dim(youtube)## [1] 40949 16The dataset contains 40,949 rows and 16 columns. Key variables
include views, likes, dislikes,
comment_count, publish_time, and
trending_date. Notably, the derived variables required for
modeling (engagement_rate, like_per_view,
video_length_category, days_since_published)
are not present in the raw data and does have to be engineered.
Phase 3: Data Validation
Missing Value Check
colSums(is.na(youtube))## video_id trending_date title
## 0 0 0
## channel_title category_id publish_time
## 0 0 0
## tags views likes
## 0 0 0
## dislikes comment_count thumbnail_link
## 0 0 0
## comments_disabled ratings_disabled video_error_or_removed
## 0 0 0
## description
## 578Only the description column contains missing values (578
NAs). All key numeric columns like views, likes, dislikes, and
comment_count are complete with zero missing values.
Variance Check
youtube %>%
select(views, likes, dislikes, comment_count) %>%
summarise(across(everything(), list(mean = mean, sd = sd, var = var)))## # A tibble: 1 × 12
## views_mean views_sd views_var likes_mean likes_sd likes_var dislikes_mean
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2360785. 7394114. 5.47e13 74267. 228885. 52388498047. 3711.
## # ℹ 5 more variables: dislikes_sd <dbl>, dislikes_var <dbl>,
## # comment_count_mean <dbl>, comment_count_sd <dbl>, comment_count_var <dbl>All four numeric columns show substantial rea world variance. Views range from 549 to over 225 million with a standard deviation of 7.39 million, confirming this isa real data which is a stark contrast to the near zero variance found in the original synthetic dataset we were using.
Derived Variable Engineering
youtube_clean <- youtube %>%
filter(video_error_or_removed == FALSE) %>%
mutate(
engagement_rate = (likes + dislikes + comment_count) / views,
like_per_view = likes / views,
days_since_published = as.numeric(
as.Date(trending_date, format = "%y.%d.%m") - as.Date(publish_time)
),
video_length_category = NA # not available in raw data
)
glimpse(youtube_clean)## Rows: 40,926
## Columns: 20
## $ video_id <chr> "2kyS6SvSYSE", "1ZAPwfrtAFY", "5qpjK5DgCt4", "p…
## $ trending_date <chr> "17.14.11", "17.14.11", "17.14.11", "17.14.11",…
## $ title <chr> "WE WANT TO TALK ABOUT OUR MARRIAGE", "The Trum…
## $ channel_title <chr> "CaseyNeistat", "LastWeekTonight", "Rudy Mancus…
## $ category_id <dbl> 22, 24, 23, 24, 24, 28, 24, 28, 1, 25, 17, 24, …
## $ publish_time <dttm> 2017-11-13 17:13:01, 2017-11-13 07:30:00, 2017…
## $ tags <chr> "SHANtell martin", "last week tonight trump pre…
## $ views <dbl> 748374, 2418783, 3191434, 343168, 2095731, 1191…
## $ likes <dbl> 57527, 97185, 146033, 10172, 132235, 9763, 1599…
## $ dislikes <dbl> 2966, 6146, 5339, 666, 1989, 511, 2445, 778, 11…
## $ comment_count <dbl> 15954, 12703, 8181, 2146, 17518, 1434, 1970, 34…
## $ thumbnail_link <chr> "https://i.ytimg.com/vi/2kyS6SvSYSE/default.jpg…
## $ comments_disabled <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ ratings_disabled <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ video_error_or_removed <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE…
## $ description <chr> "SHANTELL'S CHANNEL - https://www.youtube.com/s…
## $ engagement_rate <dbl> 0.102150796, 0.047972059, 0.049994141, 0.037835…
## $ like_per_view <dbl> 0.076869319, 0.040179297, 0.045757800, 0.029641…
## $ days_since_published <dbl> 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 1, 1, 1, 2, 2, 1,…
## $ video_length_category <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…After filtering out 23 error rows, the cleaned dataset contains
40,926 observations and 20 variables. Three of the four required derived
variables were successfully engineered.
video_length_category does remains NA as video duration is
not included in this dataset and would need to have separate YouTube
Data API lookup.
Distribution Visualizations
ggplot(youtube_clean, aes(x = views)) +
geom_histogram(bins = 50, fill = "steelblue") +
scale_x_log10(labels = scales::comma) +
labs(title = "Distribution of Views (log scale)",
x = "Views", y = "Count")
ggplot(youtube_clean, aes(x = engagement_rate)) +
geom_histogram(bins = 50, fill = "darkgreen") +
labs(title = "Distribution of Engagement Rate",
x = "Engagement Rate", y = "Count")
ggplot(youtube_clean, aes(x = days_since_published)) +
geom_histogram(bins = 50, fill = "coral") +
labs(title = "Days Since Published When Trending",
x = "Days", y = "Count")
### Key Validation Findings
The views distributiondoes follow a near normal curve on a log scale, being around around 1 million views, ranging from ~10,000 to over 100 million. This is very consistent with real YouTube trending behavior.
The engagement rate distribution is righ -skewed, with most videos falling between 0.01 and 0.05, and a long tail of high engagement outliers. This is showing real genuine user behavior.
The days since published distribution shows the vast majority of videos trend within days of publishing, with a small number of outliers trending years later. This is consistent with the way viral patterns are
Phase 1-3 Conclusion by Keith
The datasnaek YouTube trending dataset passes all data validation checks. Variance is real, missing values are minimal and derived variables have been successfully engineered for three of the four required modeling inputs. The dataset is ready to support Phase 4 model development (Random Forest and Logistic Regression) by Lennin and Nico respectively.
video_length_category remains a limitation — the team
should discuss whether to drop this variable or supplement with a
YouTube Data API lookup for video duration.
Model 1: Logit Regression - by Nico
Our consultation roll play needed probabilities to provided to those seeking to increase their engagement rate on YouTube. A Logistic Regression allows us to receive the odds that each variable will increase the engagement of a YouTube video. I kept the duplicates of video within my data set to account for the time frame of repeated interactions with videos. This provided useful insite into the “Engagement Decay” over time.
log_regr set up
library(ggplot2)1. Load data: I used Keith’s cleaned data and converted the model into one called youtube_model.csv for my specific purposes.
youtube_model_df <- read.csv("data_set/youtube_model.csv")2. Variable selection: The logistic regression required Binary variables in order for the results to be easily inheritable. The variables “high_engagement”, “high_views”, “high_likes”, “high_dislikes”, “high_comments” had a string of “Yes” and “No” that needed to be converted.
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")
})3. Logistic Regression: After the conversion the regression was run with high_engagement as the dependent variable and the rest were considered independent variables.
engagement_model <- glm(
high_engagement ~ high_views + high_likes + high_dislikes + high_comments,
data = youtube_model_df,
family = "binomial"
)4. Summary Analysis: All of our selected variables were deemed highly significant from this model which was nice to see. Our estimate or log odds gave us some insite into the probability of increasing engagement.
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: 55. YouTube Engagement Drivers Using Log Odds: This graph provides a intuitive guide as to which variables directly effect the engagement odds of a video. High likes and comments increase the odds, while High dislikes and view decrease the odds. Interesting views decrease the odds of engagements which can be explained by our original formula of likes + dislike + comments / views. Views increase exponentially compared to the the other variables which causes a large shift in the odds. So yes high view do increase your odds of disproportional engagement.
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")
)
6. Second Regression Model For Analyzing Engagement Decay: I wanted to run a second logistic regression but include the continuous variable “days_since_published”. My hopes were this would provide me with a curve that would display the odds of engagement decay overtime for a YouTube video after it has been published.
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: 87. Engagement Decay Over Time : This graph depicts the 60-day projected Engagement decline. A typical youtube video with begin having around 60% engament after it has first been published. There we will see the engaments drop consistantly as the week go by until we platue at the 2 month mark.
engagement_time_model <- glm(
high_engagement ~ high_views + high_likes + high_dislikes + high_comments + days_since_published,
data = youtube_model_df,
family = "binomial"
)
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 1 Conclusion by Nico: Our logistic regression on 40,000 US trending videos successfully mapped out the interaction dynamics that dictate video performance. Provides evidence that high likes and comments directly drive engagement. In terms of impact strength, hitting high likes threshold increases a videos engagement odds by 10x, while high comments increases those odds by 5x. Additionally, high views reduces the odds of high engagement by ~ 90% because the virality of a video introduces more passive viewers. Negative viewer sentiment plays a much smaller role in this performance, as hitting high dislikes only decreases our odds of engagement by a minor 0.4x.
Model 2: Random Forest - By Lennin
One of the questions key questions in our “social media consultant” role play we wanted to answer via modeling was “do video attributes drive engagement? if so, which?” we decided a random forest was a good approach to find out. honestly, because I have a hard time feature engineering for other algos, but for random forest/decision trees its very easy since it isn’t scale dependent, it’s just about managing splits and purity. I deduped (this is panel data and the same video can appear over many days) and just did log-transform for the view counts to manage splitting.
Setup
library(tidyverse)
library(lubridate)
library(ranger)
library(caret)
library(vip)
library(pdp)
library(patchwork)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"
)1. Load Data
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: 63512. Deduplicate
The raw dataset has ~40,000 rows but only ~6,300 unique videos. I kept the most mature snapshot per video (highest view count) so engagement rates reflect an accumulated audience rather than a single early-trending moment.
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: 63513. Feature Engineering
I took my best guess to derive the target variable and all predictors from scratch. Key decisions:
- engagement_rate = (likes + dislikes + comments) / views — drop all three component columns to prevent data leakage
- log_views the raw range for views (549 to 225M) is huge; log-transforming helps the model find clean split points
- days_to_trend — how quickly the video hit trending,
parsed from the raw
YY.DD.MMdate format - Title signals (length, word count, punctuation) extracted from the raw title string
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.136914. Train / Test Split
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: 12565. Train Random Forest
I use claude code for most all this code and it advised the
following: “use the ranger package — a fast implementation
of random forest with identical results to the base
randomForest package. For feature importance go with
permutation importance, which is more honest than
impurity-based importance: it measures how much model accuracy drops
when a feature is randomly shuffled, rather than how often a feature
appears in splits.”
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.26656. Model Evaluation
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.01532The model achieved a test R² of 0.2624, meaning it explains roughly 26% of the variance in engagement rate!! So awesome since this is only using video metadata. While not a perfect predictor, this is a meaningful signal because we have no information about content quality, thumbnail design, or creator-audience relationships, or other factors that go into the UGC ecosystem and certainly drive a ton of variance. RMSE of 0.0209 is reasonable relative to the typical engagement rate range in the dataset.
Actual vs. Predicted
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)
7. Variable Importance
This chart tells us which features the model relies on most. It does not tell us direction, so I asked claude to help me visualize that with pdp below.
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 |
8. Partial Dependence Plots
Variable importance tells us that a feature matters. Partial dependence plots (PDPs) tell us how by holding all other features constant and showing how the model’s prediction changes as one feature varies across its range.
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)
)
9. Engagement Rate by Category
PDPs for categorical features get cluttered, so we use a box plot of the actual data instead. Categories are sorted by median engagement rate.
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 — actual values from dataset",
x = NULL, y = "Engagement Rate"
) +
theme_minimal(base_size = 12)
10. Best Time to Post
Combining publish hour and day of week into a heatmap reveals optimal posting windows.
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)
Model 2 Conclusion by Lennin
Using a random forest model on ~6,300 unique US YouTube trending videos, we achieved a test R² of 0.2624, meaning our model explains about 26% of the variance in engagement rate using only video metadata (no content quality, thumbnail appeal, or audience loyalty data.) The most important predictors were title-related features (title length and word count), followed by content category, description length, publish hour, and tag count. Partial dependence plots let us go beyond just knowing that these features matter and actually see how they affect engagement. The category breakdown confirmed that some niches consistently outperform others, which is a strong signal for any creator deciding what kind of content to make.
Script Summary
The script uses tidyverse for data loading and
manipulation, lubridate for parsing messy date formats,
ranger for training the random forest, caret
for the stratified train/test split, vip for variable
importance plots, pdp for partial dependence plots, and
patchwork for combining multiple plots. The raw CSV had
~40,000 rows but only ~6,300 unique videos, so we deduplicated by
keeping each video’s highest-view snapshot. All features were engineered
from scratch. title length, word count, punctuation flags, tag count,
publish hour and day of week, days from publish to trending, description
length, and a log-transformed view count as a channel-size proxy. Likes,
dislikes, and comment count were excluded entirely since they are direct
inputs to the target variable and including them would constitute data
leakage.