STAT0030 Assessment 3 — Instructions

For this assessment you should submit online – on the course Moodle page using the link “ICA3: Click here to submit your assignment”.  Make sure none of the files contains your surname, as the marking must be anonymous. You must submit two files:

● An electronic copy of your StudentNumber.rmd file, containing your R markdown code. For example, if your student number is 18239004, your R markdown script should be saved in the file 18239004.rmd.

● A single  PDF file  named  StudentNumber.pdf  containing the  knitted output of the Rmarkdown file.   This should correspond exactly to what is produced when knitting the submitted .rmd file.

Any output within your pdf should be clearly presented and structured according to the question parts.

You are allowed to use R libraries for this assessment. However, you need to make sure that you fully understand what each function is doing and briefly explain it either in your writeup or in your code. Use of any libraries beyond what we have used in STAT0030 (or ones mentioned in these instructions) will result in a lower mark.

Answers to each question must run independently.  This implies repeating code for library and data loading in different questions, where necessary.

STAT0030 Assessment 3 – Marking guidelines

The assessment is marked out of 100, with Questions 1, 2 and 3 having a total of 40, 20 and 40 marks, respectively. The marks for each question are subdivided into the following components.

1.  For Question 1, items (a), (b), (c) and (d) have a breakdown of 5, 5, 15 and 15 marks, respectively.

2.  For Question 2, items (a) and (b) have a breakdown of 10 and 10 marks, respectively.

3.  For Question 3, items (a), (b) and (c) have a breakdown of 25, 10 and 5 marks, respec- tively.

For each coding part, graphical presentation (appropriate choice of graphs and formatting), quality of printed output (appropriate messages printed) and quality of the code (your code should be clean, readable – with sufficient commenting for the user – and efficient) are factors to be considered.

STAT0030 Assessment 3 — Questions

Background. Recommendation systems are statistical models that rank items of interest to the user of a particular service. For instance, Netflix, Spotify and Amazon use those to recommend new movies, music and products in general.  One key idea in recommendation systems is to generate rating predictions, so that high-rating items are given priority (to be combined with other factors, such as diversity of items).

For this assignment, we will consider one MovieLens database, taken out of a public reposi- tory of movie ratings across a sample of pseudo-anonymised users of https://movielens.org . The particular dataset we will use is the MovieLens 100K Dataset, which should be downloaded from https://grouplens.org/datasets/movielens/ (it’s under the “older datasets” sec- tion).   It contains 100,000 ratings from about 1000 users and 1700 movies.   Read the ap- propriately named README file for an on overview of its contents.  The main file is u.data, containing four columns, of which the relevant ones are the first three:  a  “user id”,  “movie id” and “rating”. The first two columns point to information to be found in files u.user and u.item, respectively. Ratings are given in a 1 to 5 “stars” scale.

Questions.   You will  perform several  unsupervised  learning tasks with this data,  including clustering and dimensionality reduction.  You will see how the classical principal component analysis (PCA) idea interacts with clustering and prediction, and how to build extensions of PCA.

1.  In this section, we will perform cluster analysis of users and movies in the MovieLens 100K dataset using a variety of methods.

(a)  Let’s start by loading the data.  In this question, write code to create three data frames, to be called ratings, movies and users, by loading files u.data, u.item and u.user, respectively.  Choose appropriate column names to be given to those data frames, based on the README file.  Notice that u.item will pose problems to read.table. You will need to figure out why.  Consult the documentation of that function (i.e. type ?read.table in the console and read it) in order to use a variant or particular choice of arguments that will allow you to read that file, documenting what you did.

(b)  Create a function to be called triplets to matrix that does the following. While data frame ratings describes the ratings by triplets (“user id”, “movie id”, “number of stars”), as explained in the README file, it is sometimes useful to transform this information into a matrix of all possible potential pairs of users and movies. Function triplets to matrix must take two arguments, default value and triplets, and return a matrix m as follows: m[i,  j] contains the respective rating (number of starts) given by user with id i to the movie with id j if this rating is contained in data frame triplets.  Data frame triplets is assumed to follow the same format as the first three columns of data frame ratings, constructed in the previous question. If user i exists and movie j exists, but movie pair (i, j) is not in triplets, then m[i,  j] should be assigned the value default value, which can be arbitrary.

(c)  Construct a matrix m zero which is the result of calling triplets to matrix with a default value of 0, and argument triplets set to be a data frame formed from the first three columns of ratings . Run kmeans on m zero to get a clustering of users. The choice of number of clusters is yours. Provide your own interpretation of the results, using any information of your choice out of data frames movies and/or users.  For instance, you may want to investigate whether there are meaningful patterns of gender or occupation (for users) or genre (for movies) that vary across clusters.  Your code should provide all steps of your analysis, which should be ac- companied by a short report (recommended length no longer than 300 words, plus tables and figures you may find useful).   In your report, you should discuss how meaningful you think the distance function used in kmeans is for this application, which shortcomings it may have and one suggestion of how to improve it.

(d)  Run prcomp on m zero and choose a number k of principal components, using any criteria of your choice, to form a matrix where the number of rows is the number of movies and the columns are the top k eigenvectors found by PCA. Run kmeans on this matrix to get a clustering of movies. Provide your own interpretation of the results following the template of the previous item, a report with a recommended length no longer than 300 words.

2.  In the previous question, we pretended we could put a default value into the missing entries of our user-movie matrix.  As you may imagine, there are better ways of doing that.  In fact, matrix completion is a task upon which part of the business models of Netflix and others are built: by estimating which ratings an user would give to a movie that they have not seen yet, recommendations can be made. Often simplifying assumptions are adopted, such as the missingness pattern being completely at random (that is, there are no factors that would correlate having a missing entry with the hypothetical rating it would get).

Given its widespread applications, many methods exist for matrix completion.  Here, we will use as an example the package softImpute, which you can install in the usual way. For the sake of this question, we won’t ask you to understand how the methods in this package work. However, if you are curious, the main reference is

Rahul Mazumder, Trevor Hastie and Rob Tibshirani (2010).  “Spectral Regu-     larization Algorithms for Learning Large Incomplete Matrices”, https://web.     stanford.edu/~hastie/Papers/mazumder10a.pdf , Journal of Machine Learn- ing Research 11 (2010) 2287-2322.

The R vignette for this package1  provides further background and examples.

(a) Write code for the following pipeline.   Run softImpute in your matrix m zero, now constructed using argument default value set to NA (so that softImpute knows which values to complete and which values are from the data).  Please use

softImpute with the argument type="als", which is a computationally efficient

algorithm that will save you a lot of time.   Function  complete, from the same package, will generate a completed (estimated) user-movie matrix.   Re-run your kmeans analysis on this completed matrix, and contrast results against what you obtained in the first question, writing some comments on the difference.  Which results do you think are more sensible, if any, and why?

(b) Another way of understanding matrix completion is as a prediction task where the existing labels ratings are training labels, and the missing ones are test labels. We can assess how well it does using data-splitting techniques like cross-validation: we get the full matrix completion using a method like softImpute2 , but only compute the errors for those user-movie pairs which are in the test set (so that we know which values they actually had).  Write code for assessing prediction errors using the 5 splits provided by the MovieLens 100K Dataset (they split u.item into the 5 pairs u1.train, u1.test,  . . . , u5.train, u5.test).   Perform an analysis of your choice by partitioning the test sets according to categories of interest (by user attributes like gender and occupation, or by movie genre and popularity - which you can estimate by the number of ratings they get), commenting on how errors vary among subcategories.  Take into consideration that some strata may be too small, and perhaps not suitable for this analysis (uncertainty measures would help).

3.  Now, we are going to create from scratch our own matrix completion method.   It is a highly simplified one,  but it shares some of the core ideas of many other methods. In particular, it is defines a probabilistic model, so that we can do maximum likelihood estimation of it.   For simplicity, we will consider binary  ratings only (think of it as a “thumbs up” vs.  “thumbs down” classification), but this model could easily be extended to a more complex rating system such as the one in MovieLens.

Our model consists of having parameter ui  for each user i, and parameter vj  for each movie j. Data is given as a set of binary rating variables Rij  such that

puivj   三 P (Rij  = 1; ui , vj ) =            1           

Let e be the set of all  recorded  ratings rij    ∈  R0, 1},  with n  being the  number of users and m the number of movies.  The likelihood function for the set of parameters θ 三 Ru1 , . . . , un , v1 , . . . , vm } is

L(θ; e) =

p (1 - puivj )1_rij .

(i,j)-e

Is it crazy to have the number of parameters growing with n and m? It depends how we regularise them, and we in practice expect that the size of e will be considerably higher than n + m anyway (in our MovieLens 100K example, it’s 100,000 data points against

< 3000 parameters). We will ignore regularisation for this question though.

You can see why matrix completion is better interpreted as an unsupervised learning task: there is no input-output pairing, just a set of measurements rij  for which we can model the joint distribution. We can still use unsupervised learning for prediction: for any pair (i, j), we use the estimates i  and vˆj  so that we can make a prediction for Rij  based on these estimated parameters. As a matter of fact, ui  and vj  play a role analogous to the principal component eigenvectors. There is a whole literature on the connections between PCA variations and matrix completion, and you can see Mazumder et al. for an example.

Maximum likelihood can be done by a variety of methods.  Perhaps the simplest one is gradient ascent. Verify that the partial derivative of the log-likelihood function l(θ; e) 三 log L(θ; e), with respect to some ui , is given by

∂l(θ; e)

∂ui

with an analogous formula, by symmetry, for the partial derivative ∂l(θ; e)/∂vj . We can then maximize the log-likelihood by starting from a random choice θ (0)  of parameters Ru, . . . , u, v , . . . , v} and repeatedly applying the rule

θ (k+1)  - θ (k) + η (k)  × Rθ =θ7k(l(θ; e),

for k = 0, 1, 2, 3, . . . until a maximum number of iterations is reached or the log-likelihood doesn’t change  “much” (i.e., as k increases, the value of the log-likelihood fluctuates within  pre-specified,  problem dependent,  small  number).   Here,  Rθ =θ7k(l(θ; e)  is the gradient of the  log-likelihood evaluated at  θ (k),  which we get  by stacking  its  partial derivatives with respect to u1 , . . . , un , v1 , . . . , vm . The so-called learning rate η (k) follows a schedule which is pre-defined. In what follows, we will ask of you the simplest possible implementation. The final solution will be our estimate θˆ.

(a)  Create a function, to be called learn uv, to perform maximum likelihood by gradient ascent in the model given above.  This function takes as input a data frame dat where dat[k, 1], dat[k, 2] and dat[k, 3] are the user id, movie id and (binary) rating, respectively, in the k-th triplet in a collection of triplets. Three other inputs should be passed:  eta, the learning rate that will be kept constant through all iterations; tol, a tolerance level that will stop the gradient optimisation method at iteration k if the absolute value of the difference between the log-likelihoods at iteration k and iteration k - 1 is less that tol; and max iter, the maximum number of gradient iterations that we can perform.   This function must return a list of two entries, u and v, corresponding to the respective parameters (you should also name the entries in the list as u and v).  Notice that, when testing your function, the choice of eta can either make your optimisation very slow, or eventually crash if the log- likelihood spirals to infinite. It’s not your problem to select eta in this question but. when testing it, you will probably want to make sure it works by a trial-and-error fiddling of its value.  Feel free to implement (and document!)  heuristics, such as not updating the parameters if the log-likelihood decreases at any particular update, instead halving the value of eta and trying again.

(b) Write a function predict uv that takes (estimated) parameters u, v and a set of test pairs pairs, and returns a binary vector r hat of predictions.  Entries pairs[k, 1], pairs[k,  2], correspond to a  user id and  movie id,  respectively.   We can assume that the user ids and movie ids in u, v are the indices in those vectors (e.g, u[10] is the u-feature of the user with id 10). If either pairs[k,  1] or pairs[k, 2] is larger than the length of either u or v, respectively, just set r hat[k] to NA. To binarise r hat, return 1 if only if the corresponding pˆuivj   is greater than 0.5. Apply this to the train/test split defined by u1.base, u1.test and report the results. To binarise the ratings, convert 1-3 starts to zero and 4-5 stars to 1.

(c) The induced vectors  and vˆ can actually be used to cluster users and movies. As a matter of fact, softImpute returns a list of three items, one of which is called u and another called v (let’s ignore the third one).  If you inspect what you got from Question 2, you will see that u and v are  matrices of two columns each: that’s because in general we can define an inner product ui(T)vi  that generalises the simple uivj  and allows for richer models (the dimensionality to be chosen is itself a hyperparameter, and the default in softImpute is 2).  Re-run softImpute on the complete 100K data with argument rank.max  =  1 and compare it against what you get by running your learn uv output, with comparisons of your choice.  They could include the correlation coefficient between the corresponding u (or v) vectors of each method, or to which extent kmeans applied to the respective u (or v) have similarities and differences.

HOW TO  PREPARE YOUR SUBMISSION. You must provide a single  .rmd file including all of your code and answers.  The file MUST contain exactly one chunk of code for each of the three questions,  which  implements and  runs all the required steps in that question in a self-contained way.  This means that code for one question cannot depend on code for a different question.  The clarity of the printed output will be taken into consideration.   It is up to your judgement to decide how the printed output should be organised. Answers to explanatory points within each question should appear  in a single  block of text,  after the  respective code and output.  Explanatory points should be presented at a level that can be understood easily by somebody with a MSc in Statistics. Code comments should clarify which parts of each question are being solved.  When submitting, knit the corresponding  .rmd into a  .pdf, including all code and output. Both the  .rmd and the .pdf must be submitted.

For instance, for Question 1, prepare a single chunk of code which a) implements the data loading pipeline;  b) implements function triplets to matrix;  c) performs the data processing and necessary steps for generating the outputs of the cluster analysis; iv) generates the output of the cluster analysis combined with principal component analysis. Use your own judgement of what should be plotted and printed. The code will be followed by a block of text with the reports required for Questions 1(c) and 1(d), referring to the print-outs and plots generated by the code.  For Question 1, the  .pdf file will show all of that question’s code, followed by its output, followed by its textual explanations. This will then be followed by the respective steps for Questions 2 and 3.