options(warn=-1)
library(dslabs)
library(tidyverse)
library(caret)
library(ggrepel)
library(data.table)
data("movielens")

Data Analysis

We will investigate some of machine learning recommendation systems approaches to create specific recommendations for users, based on previous rating analysis. These are based on the Netflix challenge, which occurred in 2006.
The GroupLens research lab has a database of over 27 millions ratings for over 58,000 movies by more than 280,000 users.
We use here the movielens dataset, which is a small sample of the GroupLens dataset. The movielens has 100,004 entries, with rating from 671 users in a total of 9,066 movies.
We will proceed with some data analysis first.

dim(movielens)
movielens %>%
  summarize(n_users = n_distinct(userId),
            n_movies = n_distinct(movieId))
head(movielens)

1. 100004
2. 7

A data.frame: 1 × 2

n_users	n_movies
<int>	<int>
671	9066

A data.frame: 6 × 7

	movieId	title	year	genres	userId	rating	timestamp
	<int>	<chr>	<int>	<fct>	<int>	<dbl>	<int>
1	31	Dangerous Minds	1995	Drama	1	2.5	1260759144
2	1029	Dumbo	1941	Animation|Children|Drama|Musical	1	3.0	1260759179
3	1061	Sleepers	1996	Thriller	1	3.0	1260759182
4	1129	Escape from New York	1981	Action|Adventure|Sci-Fi|Thriller	1	2.0	1260759185
5	1172	Cinema Paradiso (Nuovo cinema Paradiso)	1989	Drama	1	4.0	1260759205
6	1263	Deer Hunter, The	1978	Drama|War	1	2.0	1260759151

We look here a small subset of ratings from 6 users over 5 movies. There are many movies without ratings from users.

movieIds <- movielens %>%
  dplyr::count(movieId) %>%
  top_n(5) %>%
  pull(movieId)
tab <- movielens %>%
  filter(userId %in% c(13:20)) %>% 
  filter(movieId %in% movieIds) %>% 
  select(userId, title, rating) %>% 
  spread(title, rating)
tab

Selecting by n

A data.frame: 7 × 6

userId	Forrest Gump	Pulp Fiction	Shawshank Redemption, The	Silence of the Lambs, The	Star Wars: Episode IV - A New Hope
<int>	<dbl>	<dbl>	<dbl>	<dbl>	<dbl>
13	5.0	3.5	4.5	NA	NA
15	1.0	5.0	2.0	5.0	5.0
16	NA	NA	4.0	NA	NA
17	2.5	5.0	5.0	4.5	3.5
18	NA	NA	NA	NA	3.0
19	5.0	5.0	4.0	3.0	4.0
20	2.0	0.5	4.5	0.5	1.5

This code shows how the dataset matrix is sparse, with a random sample of 100 movies and 100 users, by showing with colored dots where there are ratings.

users <- sample(unique(movielens$userId), 100)
rafalib::mypar()
movielens %>% filter(userId %in% users) %>% 
  select(userId, movieId, rating) %>%
  mutate(rating = 1) %>%
  spread(movieId, rating) %>% select(sample(ncol(.), 100)) %>% 
  as.matrix() %>% t(.) %>%
  image(1:100, 1:100,. , xlab="Movies", ylab="Users")
abline(h=0:100+0.5, v=0:100+0.5, col = "grey")

The following distribution shows that some movies are more rated than others, since there are blockbusters watched by millions and independent movies watched by just a few.

movielens %>% 
  dplyr::count(movieId) %>% 
  ggplot(aes(n)) + 
  geom_histogram(bins = 30, color = "black") + 
  scale_x_log10() + 
  ggtitle("Movies")

In similar way, there are users that are more active than others at rating movies.

movielens %>%
  dplyr::count(userId) %>% 
  ggplot(aes(n)) + 
  geom_histogram(bins = 30, color = "black") + 
  scale_x_log10() +
  ggtitle("Users")

By plotting the number of ratings by year the movie came out, we notice there is an increase on the average around thes 90s.

options(repr.plot.width=12, repr.plot.height=7)
movielens %>%
  group_by(movieId) %>%
  summarize(n_rating = n(), year = first(year)) %>%
  ggplot(aes(year, n_rating, group=year)) + 
  geom_boxplot() +
  coord_trans(y = "sqrt") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Among the movies that came out since 1990, we show the top 10 with highest average number of ratings per year.

movielens %>% filter(year >= 1990) %>%
  group_by(movieId) %>%
  summarize(n = n(),
            years = 2018 - first(year),
            title = title[1],
            rating = mean(rating),
            nrates = n/years) %>%
  arrange(desc(nrates)) %>%
  head(10)

A tibble: 10 × 6

movieId	n	years	title	rating	nrates
<int>	<int>	<dbl>	<chr>	<dbl>	<dbl>
356	341	24	Forrest Gump	4.054252	14.20833
79132	111	8	Inception	4.049550	13.87500
2571	259	19	Matrix, The	4.183398	13.63158
296	324	24	Pulp Fiction	4.256173	13.50000
318	311	24	Shawshank Redemption, The	4.487138	12.95833
58559	121	10	Dark Knight, The	4.235537	12.10000
4993	200	17	Lord of the Rings: The Fellowship of the Ring, The	4.182500	11.76471
5952	188	16	Lord of the Rings: The Two Towers, The	4.061170	11.75000
7153	176	15	Lord of the Rings: The Return of the King, The	4.127841	11.73333
2858	220	19	American Beauty	4.236364	11.57895

By analyzing the data it is possible to notice that the most frequently rated movies tend to have higher average ratings, since more people watch popular movies. This trend can be verified in the following graphic.

movielens %>% filter(year >= 1990) %>%
  group_by(movieId) %>%
  summarize(n = n(),
            years = 2018 - first(year),
            title = title[1],
            rating = mean(rating),
            nrates = n/years) %>%
  mutate(nrates = round(nrates)) %>%
  group_by(nrates) %>%
  summarize(avg_rating = mean(rating))%>%
  qplot(nrates, avg_rating, data = .) +
  geom_smooth(formula = 'y ~ x', method = 'loess')

The movies dataset has a column with information of the movie genres. However this column includes every genre that applies to the movie and some may contain more than one genre.
Since a movie with only one genre may differ from other that incorporates several genres, we consider every combination of genres as unique and compute the average and standard error for each one with more than 1,000 ratings.
It is possible to see that there is a strong evidence of genre effect over the ratings.

options(repr.plot.width=12, repr.plot.height=10)
movielens %>%
  group_by(genres) %>%
  filter(n() > 1000) %>%
  summarize(avg = mean(rating),
            se = sd(rating)/sqrt(n())) %>%
  mutate(genres = reorder(genres, avg)) %>%
  ggplot(aes(x = genres, y = avg, ymin = avg - 2*se, ymax = avg + 2*se)) +
  geom_point() +
  geom_errorbar() +
  theme(axis.text.x = element_text(size=14, angle = 90, hjust = 1))

For the prediction algorithms we start by creating the training and test sets. We include in the test set only movies and users that also appear in the training set.

test_index <- createDataPartition(y = movielens$rating, times = 1,
                                  p = 0.2, list = FALSE)
train_set <- movielens[-test_index,]
test_set <- movielens[test_index,]

test_set <- test_set %>% 
  semi_join(train_set, by = "movieId") %>%
  semi_join(train_set, by = "userId")

We need to quantify a loss function, to compare the different machine learning algorithms. The Netflix challenge used the typical error, based on the residual mean squared error: (user $u$ and movie $i$)
$RMSE = \sqrt{\frac{1}{N}\sum_{u,i} (\hat{y}_{u,i} - y_{u,i})^{2}}$

RMSE <- function(true_ratings, predicted_ratings){
  sqrt(mean((true_ratings - predicted_ratings)^2))
}

Average classifier

Same average

We start with the simplest model. We predict the same rating for all movies, regardless of user or movie, with all differences explained by random variation:
$Y_{u,i} = \mu + \varepsilon_{u,i}$
In this model, the estimate that minimizes the residual mean squared error is the least squares estimate of $\mu$, which in this case, is the average of all ratings over all users.

mu_hat <- mean(train_set$rating)
mu_hat

3.54285517862053

We can compute the residual mean squared error for this simple model:

naive_rmse <- RMSE(test_set$rating, mu_hat)
rmse_results <- tibble(method = "Average", RMSE = naive_rmse)
rmse_results %>% knitr::kable()

|method  |     RMSE|
|:-------|--------:|
|Average | 1.047911|

Movie effects

We can improve our model by assigning different average for each movie $i$ by adding the movie effect (or bias) $b_{i}$:
$Y_{u,i} = \mu + b_{i} + \varepsilon_{u,i}$
The movie bias can be calculated by subtracting the overall average from each movie average rating.

options(repr.plot.width=7, repr.plot.height=7)
movie_avgs <- train_set %>% 
  group_by(movieId) %>% 
  summarize(b_i = mean(rating - mu_hat))
movie_avgs %>% qplot(b_i, geom ="histogram", bins = 10, data = ., color = I("black"))

Many movies have substantial $b_{i}$ estimates. Movies with $b_{i} = 1.5$ implies a perfect five-star rating, since the average is $\hat{\mu} = 3.5$
By applying the movie average to each movie, we can compute the new residual mean squared error of this model on the test set.

predicted_ratings <- mu_hat + test_set %>% 
  left_join(movie_avgs, by = "movieId") %>%
  .$b_i

model_mov_rmse <- RMSE(test_set$rating, predicted_ratings)
rmse_results <- bind_rows(rmse_results,
                          tibble(method="Movie Effect Model",
                                 RMSE = model_mov_rmse))
rmse_results %>% knitr::kable()

|method             |      RMSE|
|:------------------|---------:|
|Average            | 1.0479106|
|Movie Effect Model | 0.9900482|