Update manuscript for RFP storage of large diagonal blocks beyond [1,1]

[dmbates]

This version uses the more compact storage available in the db/FasterrankUpdate branch of MixedModels.jl

[dmbates]

@palday and @ajinkya-k Could you, at your convenience, take a look at this re-arrangment of the JDS version of the paper. I moved the discussion of the RFP format into section 3.2 because it is now the default for dense triangular matrices of dimension > 1000. To reproduce the timings you will need to use the db/FasterrankUpdate branch of the MixedModels.jl package.

I can't guarantee that this version is error-free, just that it is the version I had when I got tired of re-reading the text over and over again.

I guess we should comment on the interpretation of results for some of the movielens fits. One possible approach is to look at the "most popular" 20 movies as determined by raw mean scores and BLUPS and how that changes according the movie-cutoff value in the model fit. (I.e. are there movies that have very few but very high ratings?) Any other suggestions?

[dmbates]

Regarding how the mixed model helps with analysis of the ml32m data set, I did find that movies with very few ratings had both the highest and the lowest mean ratings

julia> sort!(combine(groupby(DataFrame(tbl), [:movieId, :mnrtngs]), :rating => mean), :rating_mean; rev=true)
84432×3 DataFrame
   Row │ movieId  mnrtngs  rating_mean 
       │ Int32    Int32    Float32     
───────┼───────────────────────────────
     1 │  257777        1          5.0
     2 │  257983        1          5.0
     3 │  240154        1          5.0
     4 │  200016        3          5.0
     5 │  196787        1          5.0
     6 │  223908        1          5.0
     7 │  221334        1          5.0
     8 │  279542        2          5.0
     9 │  205593        1          5.0
    10 │  292047        1          5.0
    11 │  226640        2          5.0
    12 │  270540        1          5.0
    13 │  117656        1          5.0
    14 │  189587        1          5.0
    15 │  206168        1          5.0
    16 │  141064        1          5.0
    17 │  172071        1          5.0
    18 │  254544        1          5.0
    19 │  276277        1          5.0
   ⋮   │    ⋮        ⋮          ⋮
 84414 │  160108        1          0.5
 84415 │  160498        1          0.5
 84416 │  160523        1          0.5
 84417 │  160742        1          0.5
 84418 │  160784        1          0.5
 84419 │  160994        1          0.5
 84420 │  161000        1          0.5
 84421 │  161004        1          0.5
 84422 │  203529        1          0.5
 84423 │  203589        1          0.5
 84424 │  203707        1          0.5
 84425 │  155577        1          0.5
 84426 │  234239        1          0.5
 84427 │  244112        1          0.5
 84428 │  131288        1          0.5
 84429 │  180659        1          0.5
 84430 │  244546        1          0.5
 84431 │  254997        1          0.5
 84432 │  150562        1          0.5
                     84394 rows omitted

[dmbates]

I have added the formulation of the gradient suggested by Claude fable and cited the Murray (2016) arxiv.org paper on using symbolic differentiation for the Cholesky factor. Unfortunately the citation style used by JDS reduces this to just the author and title without mention of the arxiv identifier. @palday or @ajinkya-k if you can work out how to fix this please do. I feel this is ready for resubmission after addition of some brief discussion of the results of the ml-32m fits. My thought was to concentrate on the most popular movies and to show that those at the extremes for simple averages have very few ratings whereas the extremes of conditional modes for the movies have many ratings. Also to show that the conditional modes for movies with many ratings are not affected strongly by changing the movie cutoff or the user cutoff.

@palday and @ajinkya-k If possible I would appreciate your comments over the weekend.

[dmbates]

I did check on the number of ratings for movies with extreme average ratings (i.e. average rating of 0.5 or 5.0) and the vast majority have only a single rating

julia> freqtable(r5, :mnrtngs)
4-element Named Vector{Int64}
mnrtngs  │ 
─────────┼─────
1        │ 1359
2        │   80
3        │    5
4        │    1

julia> size(r5)
(1445, 3)

julia> freqtable(filter(x -> x.rating_mean == 0.5, mnrtng), :mnrtngs)
5-element Named Vector{Int64}
mnrtngs  │ 
─────────┼─────
1        │ 1167
2        │   98
3        │   15
4        │    5
9        │    1

Now I need to go back to a fit of the whole data set to get the conditional modes of the movie random effects.

Does this seem like a suitable direction to respond to the question of what we learn about the movie ratings by fitting these models?