In [21]:
qplot(lambdas, rmses)
lambdas[which.min(rmses)]
4.5
Regularization is a technique that permits us to penalize large estimates that come from small sample sizes, thus shrinking the estimates torwards zero and minimizing the residual mean squared errors.
When applying the movie effect, many obscure movies get high rating prediction, but in fact, they have low ratings, thus showing large mistake residuals in prediction.
test_set %>%
left_join(movie_avgs, by = "movieId") %>%
mutate(residual = rating - (mu_hat + b_i)) %>%
arrange(desc(abs(residual))) %>%
select(title, residual) %>% slice(1:10) %>% knitr::kable()
|title | residual| |:--------------------------------------------|--------:| |In July (Im Juli) | -4.5| |Psycho II | -4.5| |Day of the Beast, The (Día de la Bestia, El) | 4.5| |The Secret Life of Pets | 4.0| |Faust | -4.0| |Walk on the Moon, A | -4.0| |Beginners | 4.0| |High School Musical 3: Senior Year | 4.0| |Extremities | -4.0| |On the Ropes | -4.0|
If we look at the top 10 best and worse predicted movie effect $b_{i}$, we see that both have many obscure movies. Many of them were rated by very few people, in most cases just one. With few users ratings, we have more uncertainty, therefore larger $b_{i}$ are more likely to happen.
movie_titles <- movielens %>%
select(movieId, title) %>%
distinct()
train_set %>% dplyr::count(movieId) %>%
left_join(movie_avgs, by = "movieId") %>%
left_join(movie_titles, by = "movieId") %>%
arrange(desc(b_i)) %>%
select(title, b_i, n) %>%
slice(1:10) %>%
knitr::kable()
train_set %>% dplyr::count(movieId) %>%
left_join(movie_avgs, by = "movieId") %>%
left_join(movie_titles, by = "movieId") %>%
arrange(b_i) %>%
select(title, b_i, n) %>%
slice(1:10) %>%
knitr::kable()
|title | b_i| n| |:-------------------------------------------------------|--------:|--:| |When Night Is Falling | 1.457145| 1| |Lamerica | 1.457145| 1| |Picture Bride (Bijo photo) | 1.457145| 1| |Red Firecracker, Green Firecracker (Pao Da Shuang Deng) | 1.457145| 3| |Faces | 1.457145| 1| |Maya Lin: A Strong Clear Vision | 1.457145| 2| |Heavy | 1.457145| 1| |Gate of Heavenly Peace, The | 1.457145| 1| |Day the Sun Turned Cold, The (Tianguo niezi) | 1.457145| 1| |Vive L'Amour (Ai qing wan sui) | 1.457145| 1|
|title | b_i| n| |:-------------------------------------------------|---------:|--:| |Santa with Muscles | -3.042855| 1| |B*A*P*S | -3.042855| 1| |Indian Summer (a.k.a. Alive & Kicking) | -3.042855| 1| |Merry War, A | -3.042855| 1| |Day of the Beast, The (Día de la Bestia, El) | -3.042855| 1| |Whiteboyz | -3.042855| 1| |Catfish in Black Bean Sauce | -3.042855| 1| |Horrors of Spider Island (Ein Toter Hing im Netz) | -3.042855| 1| |Monkeybone | -3.042855| 1| |Arthur 2: On the Rocks | -3.042855| 1|
These obscure movies are considered noisy estimates (outliers) and should not be considered in our models. We solve this problem by applying regularization, since the data has many outliers.
The new movie effect $b_{i}$ is estimated by minimizing the following equation, which considers the RMSE and a penalty term:
$\frac{1}{N} \sum_{u,i} (y_{u,i} - \mu - b_{i})^{2} + \lambda\sum_{i} b_{i}^{2}$
So with regularization, the estimate $b_{i}$ that minimizes the equation is given by:
$\hat{b}_{i}(\lambda) = \frac{1}{\lambda + n_{i}} \sum_{i=1}^{n_{i}} (Y_{u,i} - \hat{\mu})$
Since $\lambda$ is a tunning parameter of our model, we will apply 10-fold cross validation to find the best value for it.
lambdas <- seq(0, 10, 0.25)
folds <- createFolds(train_set$rating, k = 10, list = TRUE, returnTrain = FALSE)
rmses <- sapply(lambdas, function(l){
rmse <- sapply(folds, function(i){
train_set_star <- train_set[i,]
test_set_star <- train_set[-i,] %>%
semi_join(train_set_star, by = "movieId") %>%
semi_join(train_set_star, by = "userId")
mu <- mean(train_set_star$rating)
sum <- train_set_star %>%
group_by(movieId) %>%
summarize(s = sum(rating - mu), n_i = n())
predicted_ratings <- test_set_star %>%
left_join(sum, by = "movieId") %>%
mutate(b_i = s/(n_i+l)) %>%
mutate(pred = mu + b_i) %>%
.$pred
return(RMSE(test_set_star$rating, predicted_ratings))
})
return(mean(rmse))
})
qplot(lambdas, rmses)
lambdas[which.min(rmses)]
The tunned $\lambda$ with the lowest residual mean squared error is $\lambda = 4.5$. We use this $\lambda$ to fit the regularized model with the entire training set.
lambda <- lambdas[which.min(rmses)]
movie_reg_avgs <- train_set %>%
group_by(movieId) %>%
summarize(b_i = sum(rating - mu_hat)/(n()+lambda), n_i = n())
To see the regularization effect, we plot the regularized estimate versus the least square estimates, with the size of the circle showing the size of n_i.
options(repr.plot.width=12, repr.plot.height=7)
data_frame(original = movie_avgs$b_i,
regularized = movie_reg_avgs$b_i,
n = movie_reg_avgs$n_i) %>%
ggplot(aes(original, regularized, size=sqrt(n))) +
geom_point(shape=1, alpha=0.5)
It is possible to see that for movies that have only 1 rating, the regularized estimate $b_{i}$ is much smaller than the original estimate, shrinking towards zero.
Applying regularization on the data, we look for the new top 10 best movies (highest $b_{i}$) and 10 worst movies (lowest $b_{i}$).
train_set %>%
dplyr::count(movieId) %>%
left_join(movie_reg_avgs, by = "movieId") %>%
left_join(movie_titles, by = "movieId") %>%
arrange(desc(b_i)) %>%
select(title, b_i, n) %>%
slice(1:10) %>%
knitr::kable()
train_set %>%
dplyr::count(movieId) %>%
left_join(movie_reg_avgs, by = "movieId") %>%
left_join(movie_titles, by = "movieId") %>%
arrange(b_i) %>%
select(title, b_i, n) %>%
slice(1:10) %>%
knitr::kable()
|title | b_i| n| |:----------------------------|---------:|---:| |Shawshank Redemption, The | 0.9249504| 238| |Godfather, The | 0.9006195| 169| |All About Eve | 0.8076476| 31| |Usual Suspects, The | 0.8023610| 159| |Maltese Falcon, The | 0.7955457| 50| |On the Waterfront | 0.7602922| 25| |Roger & Me | 0.7543617| 34| |Godfather: Part II, The | 0.7469407| 104| |African Queen, The | 0.7416030| 41| |Best Years of Our Lives, The | 0.7290654| 10|
|title | b_i| n| |:-----------------------|---------:|--:| |Battlefield Earth | -1.789188| 14| |Super Mario Bros. | -1.599999| 14| |Joe's Apartment | -1.547825| 7| |Speed 2: Cruise Control | -1.545188| 20| |Little Nicky | -1.383782| 14| |Wild Wild West | -1.363509| 37| |Mod Squad, The | -1.285713| 5| |10,000 BC | -1.285713| 5| |Anaconda | -1.281703| 24| |Spice World | -1.270935| 10|
By calculating the residual mean squared error over the test set, we can see that we improved our model.
predicted_ratings <- test_set %>%
left_join(movie_reg_avgs, by = "movieId") %>%
mutate(pred = mu_hat + b_i) %>%
.$pred
model_reg_mov_rmse <- RMSE(test_set$rating, predicted_ratings)
rmse_results <- bind_rows(rmse_results,
tibble(method="Regularized Movie Effect Model",
RMSE = model_reg_mov_rmse))
rmse_results %>% knitr::kable()
|method | RMSE| |:------------------------------|---------:| |Average | 1.0479106| |Movie Effect Model | 0.9900482| |Regularized Movie Effect Model | 0.9672495|
The next step is to apply the average for users, with the users effect $b_{u}$, in addition to the movie effect $b_{i}$.
$Y_{u,i} = \mu + b_{i} + b_{u} + \varepsilon_{u,i}$
We select only those that have rated over 100 movies. Below we can see through their histogram that the users also has a significant variability.
options(repr.plot.width=7, repr.plot.height=7)
train_set %>%
group_by(userId) %>%
summarize(user_rating = mean(rating)) %>%
filter(n()>=100) %>%
ggplot(aes(user_rating)) +
geom_histogram(bins = 30, color = "black")
The user bias can be calculated by subtracting the overall average and movie effect from each user average rating.
user_avgs <- train_set %>%
left_join(movie_avgs, by = "movieId") %>%
group_by(userId) %>%
summarize(b_u = mean(rating - mu_hat - b_i))
user_avgs %>% qplot(b_u, geom ="histogram", bins = 10, data = ., color = I("black"))
Below we calculate the residual mean squared error with the two bias.
predicted_ratings <- test_set %>%
left_join(movie_avgs, by = "movieId") %>%
left_join(user_avgs, by = "userId") %>%
mutate(pred = mu_hat + b_i + b_u) %>%
.$pred
model_mov_usr_rmse <- RMSE(test_set$rating, predicted_ratings)
rmse_results <- bind_rows(rmse_results,
tibble(method="Movie + User Effects Model",
RMSE = model_mov_usr_rmse ))
rmse_results %>% knitr::kable()
|method | RMSE| |:------------------------------|---------:| |Average | 1.0479106| |Movie Effect Model | 0.9900482| |Regularized Movie Effect Model | 0.9672495| |Movie + User Effects Model | 0.9102831|
We can apply regularization to the users effect as well. The new equation to be minimized will have both movie and user effect in the penalty term:
$\frac{1}{N} \sum_{u,i} (y_{u,i} - \mu - b_{i} - b_{u})^{2} + \lambda(\sum_{i} b_{i}^{2} + \sum_{u} b_{u}^{2})$
So with regularization, the estimates $b_{u}$ and $b_{i}$ that minimize the equation are given by:
$\hat{b}_{i}(\lambda) = \frac{1}{\lambda + n_{i}} \sum_{i=1}^{n_{i}} (Y_{u,i} -
\hat{\mu})$
$\hat{b}_{u}(\lambda) = \frac{1}{\lambda + n_{u}} \sum_{u=1}^{n_{u}} (Y_{u,i} -
\hat{\mu} - \hat{b}_{i})$
Since $\lambda$ is a tunning parameter of our model, we will apply 10-fold cross validation to find the best value for it.
lambdas <- seq(0, 10, 0.25)
folds <- createFolds(train_set$rating, k = 10, list = TRUE, returnTrain = FALSE)
rmses <- sapply(lambdas, function(l){
rmse <- sapply(folds, function(i){
train_set_star <- train_set[i,]
test_set_star <- train_set[-i,] %>%
semi_join(train_set_star, by = "movieId") %>%
semi_join(train_set_star, by = "userId")
mu <- mean(train_set_star$rating)
b_i <- train_set_star %>%
group_by(movieId) %>%
summarize(b_i = sum(rating - mu)/(n()+l))
b_u <- train_set_star %>%
left_join(b_i, by = "movieId") %>%
group_by(userId) %>%
summarize(b_u = sum(rating - mu - b_i)/(n()+l))
predicted_ratings <- test_set_star %>%
left_join(b_i, by = "movieId") %>%
left_join(b_u, by = "userId") %>%
mutate(pred = mu + b_i + b_u) %>%
.$pred
return(RMSE(test_set_star$rating, predicted_ratings))
})
return(mean(rmse))
})
qplot(lambdas, rmses)
lambdas[which.min(rmses)]
For both movie and user effect, the optimal tunned parameter is $\lambda = 4.75$. Now we use the entire training set to fit the model and predict it on the test set. It can be observed that the residual mean squared error improved a little.
lambda <- lambdas[which.min(rmses)]
b_i <- train_set %>%
group_by(movieId) %>%
summarize(b_i = sum(rating - mu_hat)/(n()+lambda))
b_u <- train_set %>%
left_join(b_i, by = "movieId") %>%
group_by(userId) %>%
summarize(b_u = sum(rating - b_i - mu_hat)/(n()+lambda))
predicted_ratings <- test_set %>%
left_join(b_i, by = "movieId") %>%
left_join(b_u, by = "userId") %>%
mutate(pred = mu_hat + b_i + b_u) %>%
.$pred
model_reg_mov_usr_rmse <- RMSE(test_set$rating, predicted_ratings)
rmse_results <- bind_rows(rmse_results,
tibble(method="Regularized Movie + User Effect Model",
RMSE = model_reg_mov_usr_rmse))
rmse_results %>% knitr::kable()
|method | RMSE| |:-------------------------------------|---------:| |Average | 1.0479106| |Movie Effect Model | 0.9900482| |Regularized Movie Effect Model | 0.9672495| |Movie + User Effects Model | 0.9102831| |Regularized Movie + User Effect Model | 0.8824529|