10M MovieLens dataset

We worked until now with a small sample of the total GroupLens dataset of movie ratings. As presented in the beginning, we worked with 100,004 entries, with rating from 671 users to 9,066 movies.

We will work now with a bigger set of the data, the 10M version, which approximates more to the full dataset. The 10M version contains 10 million entries, with rating from 72,000 users to 10,000 movies. The size of the dataset poses a new challenge, since it requires much more processing power from the computer.

Our goal will be applying an algorithm to study the internal structure of the data and predict ratings for movies, based on available data of other users and movies ratings and minimize the residual mean squared error of the predictions.
We start by downloading this set and doing some preprocessing.

dl <- tempfile()
download.file("http://files.grouplens.org/datasets/movielens/ml-10m.zip", dl)

ratings <- fread(text = gsub("::", "\t", readLines(unzip(dl, "ml-10M100K/ratings.dat"))),
                 col.names = c("userId", "movieId", "rating", "timestamp"))

movies <- str_split_fixed(readLines(unzip(dl, "ml-10M100K/movies.dat")), "\\::", 3)
colnames(movies) <- c("movieId", "title", "genres")

movies <- as.data.frame(movies) %>% mutate(movieId = as.numeric(movieId),
                                           title = as.character(title),
                                           genres = as.character(genres))

movielens <- left_join(ratings, movies, by = "movieId")

Creation of the training and testing set. The test set will be 10% of the MovieLens data.

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

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

# Add rows removed from validation set back into training set
removed <- anti_join(temp, test_set)
train_set <- rbind(train_set, removed)

rm(dl, ratings, movies, test_index, temp, removed)

Joining, by = c("userId", "movieId", "rating", "timestamp", "title", "genres")

We can look the size of the training dataset.

dim(train_set)

1. 9000060
2. 6

From previous analysis, the Movie and User Regularized model achieved the highest accuracy and will be applied to the 10M dataset.
Our first step will be applying 10-fold cross validation to tune the best $\lambda$ parameter for our regularization model.

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)]

4.75

The optimal tunned parameter is $\lambda = 4.75$. We apply this $\lambda$ to the entire training set to fit the model and predict it on the test set.

lambda <- lambdas[which.min(rmses)]
mu_hat <- mean(train_set$rating)
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 - mu_hat - b_i)/(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_rmse <- RMSE(test_set$rating, predicted_ratings)
rmse_results <- tibble(method = "Regularized Movie + User Effect Model",
                       RMSE = model_reg_rmse)
rmse_results %>% knitr::kable()

|method                                |      RMSE|
|:-------------------------------------|---------:|
|Regularized Movie + User Effect Model | 0.8640345|

We achieved a very good RMSE. In comparison to the small dataset, the 10M set presented a better result.

Movie, User and Genre effect

In previous analysis we verified that there is correlation between the rating and movie genres. This was reinforced by looking at the matrix factorization.
We now try to improve our prediction model by adding the $b_{g}$ genre effect.
The size of the dataset prevent us from using time and memory consuming actions, so we will proceed with the genres as their original info. We group together each different genre or combination of genres and look for their average values, applying to the model.
At the end we also apply regularization to the data.
We start by tunning the $\lambda$ parameter with 10-fold cross validation.

rm(rmses, lambda, b_i, b_u, predicted_ratings)

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))
    b_g <- train_set_star %>%
      left_join(b_i, by = "movieId") %>%
      left_join(b_u, by = "userId") %>%
      group_by(genres) %>%
      summarize(b_g = sum(rating - mu - b_i - b_u)/(n()+l))

    predicted_ratings <- test_set_star %>%
      left_join(b_i, by = "movieId") %>%
      left_join(b_u, by = "userId") %>%
      left_join(b_g, by = "genres") %>%
      mutate(pred = mu + b_i + b_u + b_g) %>%
      .$pred
    return(RMSE(test_set_star$rating, predicted_ratings))
  })
  return(mean(rmse))
})

qplot(lambdas, rmses)
lambdas[which.min(rmses)]

4.75

The best $\lambda$ tunned is 4.75. We now test our prediction on the test set.

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 - mu_hat - b_i)/(n()+lambda))
b_g <- train_set %>%
  left_join(b_i, by = "movieId") %>%
  left_join(b_u, by = "userId") %>%
  group_by(genres) %>%
  summarize(b_g = sum(rating - mu_hat - b_i - b_u)/(n()+lambda))

predicted_ratings <- test_set %>% 
  left_join(b_i, by = "movieId") %>%
  left_join(b_u, by = "userId") %>%
  left_join(b_g, by = "genres") %>%
  mutate(pred = mu_hat + b_i + b_u + b_g) %>%
  .$pred

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

|method                                        |      RMSE|
|:---------------------------------------------|---------:|
|Regularized Movie + User Effect Model         | 0.8640345|
|Regularized Movie + User + Genre Effect Model | 0.8637239|

Movies recommendation application

Including the genres as predictors in the model, we can try some direct approach.
One application is analyzing a user profile, by finding which genres a specific user watched more and applying the model only to those movies that include the most watched genre(s).
Other application is to apply the model in the search tool, for each individual genre. Here we will present an example, with a user searching for romance movies to watch. We look at the dataset and show those with the highest prediction.

movielens %>%
  filter(str_detect(genres, "Romance")) %>%
  left_join(b_i, by = "movieId") %>%
  left_join(b_u, by = "userId") %>%
  left_join(b_g, by = "genres") %>%
  mutate(pred = mu_hat + b_i + b_u + b_g) %>%
  arrange(desc(pred)) %>%
  distinct(title) %>%
  head(10)

A data.table: 10 × 1

title
<chr>
Eternal Sunshine of the Spotless Mind (2004)
Annie Hall (1977)
Beautiful Thing (1996)
Postman, The (Postino, Il) (1994)
Léon: The Professional (Léon) (Professional, The) (1994)
Good Will Hunting (1997)
Sense and Sensibility (1995)
Forrest Gump (1994)
Shrek (2001)
Casablanca (1942)