largescaleobserved.qmd

---
fig-width: 4
fig-height: 3
fig-dpi: 150
fig-format: png
engine: julia
julia:
  exeflags: 
    - --project
---

\newcommand{\bbtheta}{{\boldsymbol\theta}}

# A large-scale observational study {#sec-largescaleobserved}

Load the packages to be used,

```{julia}
#| code-fold: true
#| output: false
using AlgebraOfGraphics
using CairoMakie
using CategoricalArrays
using DataFrames
using EmbraceUncertainty: dataset
using GLM                  # for the lm function
using MixedModels
using MixedModelsMakie
using SparseArrays         # for the nnz function
using Statistics           # for the mean function
using TidierPlots
using TypedTables
```

and define some constants and a utility function

```{julia}
#| code-fold: true
#| output: false
#| label: constants05
optsumdir(paths::AbstractString...) =
  joinpath(@__DIR__, "optsums", paths...)
@isdefined(contrasts) || const contrasts = Dict{Symbol,Any}()
@isdefined(progress) || const progress = false
```

In the previous chapter we explored and fit models to data from a large-scale designed experiment.
Here we consider an observational study - ratings of movies by users of [movielens.org](https://movielens.org) made available at the [grouplens.org download site](https://grouplens.org/datasets/movielens) [@Harper2016].

We analyze the *MovieLens 25M Dataset* of roughly 25 million ratings of movies by over 162,000 users of about 59,000 movies.

Because of constraints on the redistribution of these data, the first time that `dataset(:ratings)` or `dataset(:movies)` is executed, the 250 Mb zip file is downloaded from the grouplens data site, expanded, and the tables created as Arrow files.

This operation can take a couple of minutes.

## Structure of the data

One of the purposes of this chapter is to study which dimensions of the data have the greatest effect on the amount of memory used to represent the model and the time required to fit a model.

The two datasets examined in @sec-largescaledesigned, from the English Lexicon Project [@Balota_2007], consist of `trials` or `instances` associated with `subject` and `item` factors.
The `subject` and `item` factors are "incompletely crossed" in that each item occurred in trials with several subjects, but not all subjects, and each subject responded to several different items, but not to all items.

Similarly in the movie ratings, each instance of a rating is associated with a user and with a movie, and these factors are incompletely crossed.

```{julia}
ratings = dataset(:ratings)
```

Convert this Arrow table to a `Table` and drop the `timestamp` column, which we won't be using.
(For small data sets dropping such columns is not important but with over 25 million ratings it does help to drop unnecessary columns to save some memory space.)

```{julia}
ratings =
  Table(getproperties(ratings, (:userId, :movieId, :rating)))
```

Information on the movies is available in the `movies` dataset.

```{julia}
movies = Table(dataset(:movies))
```

In contrast to data from a designed experiment, like the English Lexicon Project, the data from this observational study are extremely unbalanced with respect to the observational grouping factors, `userId` and `movieId`.
The `movies` table includes an `nrtngs` column that gives the number of ratings for each movie, which varies from 1 to over 80,000.

```{julia}
extrema(movies.nrtngs)
```

The number of ratings per user is also highly skewed

```{julia}
users = Table(
  combine(
    groupby(DataFrame(ratings), :userId),
    nrow => :urtngs,
    :rating => mean => :umnrtng,
  ),
)
```

```{julia}
extrema(users.urtngs)
```

This selection of ratings was limited to users who had provided at least 20 ratings.

One way of visualizing the imbalance in the number of ratings per movie or per user is as an [empirical cumulative distribution function](https://en.wikipedia.org/wiki/Empirical_distribution_function) (ecdf) plot, which is a "stair-step" plot where the vertical axis is the proportion of observations less than or equal to the corresponding value on the horizontal axis.
Because the distribution of the number of ratings per movie or per user is so highly skewed in the low range we use a logarithmic horizontal axis in @fig-nrtngsecdf.

```{julia}
#| code-fold: true
#| fig-cap: "Empirical distribution plots of the number of ratings per movie and per user.  The horizontal axes are on a logarithmic scale."
#| label: fig-nrtngsecdf
#| warning: false
let
  f = Figure(; size=(700, 350))
  xscale = log10
  xminorticksvisible = true
  xminorgridvisible = true
  yminorticksvisible = true
  xminorticks = IntervalsBetween(10)
  ylabel = "Relative cumulative frequency"
  nrtngs = sort(movies.nrtngs)
  ecdfplot(
    f[1, 1],
    nrtngs;
    npoints=last(nrtngs),
    axis=(
      xlabel="Number of ratings per movie (logarithmic scale)",
      xminorgridvisible,
      xminorticks,
      xminorticksvisible,
      xscale,
      ylabel,
      yminorticksvisible,
    ),
  )
  urtngs = sort(users.urtngs)
  ecdfplot(
    f[1, 2],
    urtngs;
    npoints=last(urtngs),
    axis=(
      xlabel="Number of ratings per user (logarithmic scale)",
      xminorgridvisible,
      xminorticks,
      xminorticksvisible,
      xscale,
      yminorticksvisible,
    ),
  )
  f
end
```

In this collection of about 25 million ratings, nearly 20% of the movies are rated only once and over half of the movies were rated 6 or fewer times.

```{julia}
count(≤(6), movies.nrtngs) / length(movies.nrtngs)
```

The ecdf plot of the number of ratings per user shows a similar pattern to that of the movies --- a few users with a very large number of ratings and many users with just a few ratings.

For example, the minimum number or movies rated is 20 (due to the inclusion constraint); the median number of movies rated is around 70; but the maximum is over 32,000 (which is a lot of movies - over 30 years of 3 movies per day every day - if this user actually watched all of them).

Movies with very few ratings provide little information about overall trends or even about the movie being rated.
We can imagine that the "shrinkage" of random effects for movies with just a few ratings pulls their adjusted rating strongly towards the overall average.

Furthermore, the distribution of ratings for movies with only one rating is systematically lower than the distribution of ratings for movies rated at least five times.

First, add `nrtngs` and `urtngs` as columns of the `ratings` table.

```{julia}
ratings = Table(
  disallowmissing!(
    leftjoin!(
      leftjoin!(
        DataFrame(ratings),
        select!(DataFrame(movies), :movieId, :nrtngs);
        on=:movieId,
      ),
      select!(DataFrame(users), :userId, :urtngs);
      on=:userId,
    ),
  ),
)
```

then create a bar chart of the two distributions

```{julia}
#| code-fold: true
#| fig-cap: "Distribution of ratings for movies with only one rating compared to movies with at least 5 ratings"
#| label: fig-ratingsbarcharts
#| warning: false
let
  fiveplus = zeros(Int, 10)
  onlyone = zeros(Int, 10)
  for (i, n) in
      zip(round.(Int, Float32(2) .* ratings.rating), ratings.nrtngs)
    if n > 4
      fiveplus[i] += 1
    elseif isone(n)
      onlyone[i] += 1
    end
  end
  fiveprop = fiveplus ./ sum(fiveplus)
  oneprop = onlyone ./ sum(onlyone)
  draw(
    data((;
      props=vcat(fiveprop, oneprop),
      rating=repeat(0.5:0.5:5.0; outer=2),
      nratings=repeat(["≥5", "only 1"]; inner=10),
    )) *
    mapping(
      :rating => nonnumeric,
      :props => "Proportion of ratings";
      color=:nratings => "Ratings/movie",
      dodge=:nratings,
    ) *
    visual(BarPlot);
    figure=(; size=(600, 400)),
  )
end
```

Similarly, users who rate very few movies add little information, even about the movies that they rated, because there isn't sufficient information to contrast a specific rating with a typical rating for the user.

One way of dealing with the extreme imbalance in the number of observations per user or per movie is to set a threshold on the number of observations for a user or a movie to be included in the data used to fit the model.

To be able to select ratings according to the number of ratings per user and the number of ratings per movie, we left-joined the `movies.nrtngs` and `users.urtngs` columns into the `ratings` data frame.

```{julia}
describe(DataFrame(ratings), :mean, :min, :median, :max, :nunique, :nmissing, :eltype)
```

::: {.callout-note collapse="true"}

### Seemingly inconsistent medians of "nrtngs" and "urtngs"

The medians in this table of `nrtngs` and `urtngs` are much higher than the values from the `movies` and `users` tables because a movie with 98,000 ratings occurs 98,000 times in this table whereas it occurs only once in the `movies` table.
:::

## Models fit with lower bounds on ratings per user and per movie {#sec-lrgobsmods}

We fit a simple model to this dataset using different thresholds on the number of ratings per movie and the number of ratings per user.
These fits were performed on compute servers with generous amounts of memory (128 GiB/node) and numbers of compute cores (48/node).
A sample fit is shown in @sec-lrgobsmemprint.

The results are summarized in the following table

```{julia}
#| code-fold: true
sizespeed = DataFrame(dataset(:sizespeed))
sizespeed.ucutoff = collect(sizespeed.ucutoff)  # temporary fix to get TidierPlots to work
sizespeed
```

In this table, `mc` is the "movie cutoff" (i.e. the threshold on the number of ratings per movie); `uc` is the user cutoff (threshold on the number of ratings per user); `nratings`, `nusers` and `nmvie` are the number of ratings, users and movies in the resulting trimmed data set; `modelsz` is the size (in GiB) of the model fit; `L22sz` is the size of the [2,2] block of the `L` matrix in that model; `fittime` is the time (in seconds) required to fit the model; `nev` is the number of function evaluations until convergence; and `evtime` is the time (s) per function evaluation.

The "[2,2] block of the `L` matrix" is described in @sec-lrgobsmemprint.

### Dimensions of the model versus cut-off values

As shown if @fig-nratingsbycutoff, the number of ratings varies from a little over 21 million to 25 million, and is mostly associated with the user cutoff.

```{julia}
#| code-fold: true
#| label: fig-nratingsbycutoff
#| fig-cap: "Number of ratings in reduced table by movie cutoff value"
ggplot(sizespeed, aes(; x=:mc, y=:nratings, color=:ucutoff)) +
    geom_point() +
    geom_line() +
    labs(x="Minimum number of ratings per movie", y="Number of ratings")
```

For this range of choices of cutoffs, the user cutoff has more impact on the number of ratings in the reduced dataset than does the movie cutoff.

A glance at the table shows that the number of users, `nusers`, is essentially a function of only the user cutoff, `uc` (the one exception being at `mc = 50` and `uc=20`).

@fig-nusersbycutoff shows the similarly unsurprising result that the number of movies in the reduced table is essentially determined by the movie cutoff.

```{julia}
#| code-fold: true
#| label: fig-nusersbycutoff
#| fig-cap: "Number of users in reduced table by movie cutoff value"
ggplot(sizespeed, aes(; x=:mc, y=:nmvie, color=:ucutoff)) +
    geom_point() +
    geom_line() +
    labs(x="Minimum number of ratings per movie", y="Number of movies in table")
```

```{julia}
#| code-fold: true
describe(DataFrame(sizespeed), :mean, :min, :median, :max, :nunique, :eltype)
```

### Memory footprint of the model representation {#sec-lrgobsmemprint}

To explain what "the [2,2] block of the `L` matrix" is and why its size is important, we provide a brief overview of the evaluation of the "profiled" log-likelihood for a `LinearMixedModel` representation.

To make the discussion concrete we consider one of the models represented in this table, with cut-offs of 80 ratings per user and 50 ratings per movie.
This, and any of the models shown in the table, can be restored in a few minutes from the saved `optsum` values, as opposed to taking up to two hours to perform the fit.

```{julia}
#| code-fold: true
function ratingsoptsum(
  mcutoff::Integer,
  ucutoff::Integer;
  data=Table(ratings),
  form=@formula(rating ~ 1 + (1 | userId) + (1 | movieId)),
  contrasts=contrasts,
)
  optsumfnm = optsumdir(
    "mvm$(lpad(mcutoff, 2, '0'))u$(lpad(ucutoff, 2, '0')).json",
  )
  model = LinearMixedModel(
      form,
      filter(
        r -> (r.nrtngs ≥ mcutoff) & (r.urtngs ≥ ucutoff),
        data,
      );
      contrasts,
    )
  isfile(optsumfnm) && return restoreoptsum!(model, optsumfnm)

  @warn "File $optsumfnm is not available, fitting model."
  model.optsum.initial .= 0.5
  fit!(model; thin=1)
  saveoptsum(optsumfnm, model)
  return model
end
mvm50u80 = ratingsoptsum(50, 80)
print(mvm50u80)
```

Creating the model representation and restoring the optimal parameter values can take a couple of minutes because the objective is evaluated twice --- at the initial parameter values and at the final parameter values --- during the call to `restoreoptsum!`.

Each evaluation of the objective, which requires setting the value of the parameter $\bbtheta$ in the numerical representation of the model, updating the blocked Cholesky factor, $\mathbf{L}$, and evaluating the scalar objective value from this factor, takes a little over a minute (71 seconds) on a server node and probably longer on a laptop.

The lower triangular `L` factor is large but sparse.
It is stored in six blocks of dimensions and types as shown in

```{julia}
BlockDescription(mvm50u80)
```

This display gives the types of two blocked matrices: `A` which is derived from the data and does not depend on the parameters, and `L`, which is derived from `A` and the $\bbtheta$ parameter.
The only difference in their structures is in the [2,2] block, which is diagonal in `A` and a dense, lower triangular matrix in `L`.

The memory footprint (bytes) of each of the blocks is

```{julia}
#| code-fold: true
let
  block = String[]
  for i in 1:3
    for j in 1:i
      push!(block, "[$i,$j]")
    end
  end
  Table((;
    block,
    Abytes=Base.summarysize.(mvm50u80.A),
    Lbytes=Base.summarysize.(mvm50u80.L),
  ))
end
```

resulting in total memory footprints (GiB) of

```{julia}
#| code-fold: true
NamedTuple{(:A, :L)}(
  Base.summarysize.(getproperty.(Ref(mvm50u80), (:A, :L))) ./ 2^30,
)
```

That is, `L` requires roughly 10 times the amount of storage as does `A`, and that difference is entirely due to the different structure of the [2,2] block.

This phenomenon of the Cholesky factor requiring more storage than the sparse matrix being factored is described as [fill-in](https://en.wikipedia.org/wiki/Sparse_matrix).

Note that although the dimensions of the [2,1] block are larger than those of the [2,2] block its memory footprint is smaller because it is a sparse matrix.
The matrix is about 98% zeros or, equivalently, a little over 2% nonzeros,

```{julia}
let
  L21 = mvm50u80.L[2]  # blocks are stored in a one-dimensional array
  nnz(L21) / length(L21)
end
```

which makes the sparse representation much smaller than the dense representation.

This fill-in of the [2,2] block leads to a somewhat unintuitive conclusion.
The memory footprint of the model representation depends strongly on the number of movies, less strongly on the number of users and almost not at all on the number of ratings [@fig-memoryfootprint].

```{julia}
#| code-fold: true
#| fig-cap: Memory footprint of the model representation by minimum number of ratings per user and per movie."
#| label: fig-memoryfootprint
ggplot(sizespeed, aes(x=:mc, y=:modelsz, color=:ucutoff)) +
    geom_point() +
    geom_line() +
    labs(x="Minimum number of ratings per movie", y="Size of model (GiB)")
```

@fig-memoryvsl22 shows the dominance of the `[2, 2]` block of `L` in the overall memory footprint of the model

```{julia}
#| code-fold: true
#| fig-cap: Memory footprint of the model representation (GiB) versus the size of the [2, 2] block of L (GiB)
#| label: fig-memoryvsl22
ggplot(sizespeed, aes(; x=:L22sz, y=:modelsz, color=:ucutoff)) +
    geom_point() +
    geom_line() +
    labs(y="Size of model representation (GiB)", x="Size of [2,2] block of L (GiB)")
```

@fig-memoryfootprint shows that when all the movies are included in the data to which the model is fit (i.e. `mc == 1`) the total memory footprint is over 20 GiB, and nearly 90% of that memory is that required for the `[2,2]` block of `L`.
Even when requiring a minimum of 50 ratings per movie, the `[2,2]` block of `L` is over 30% of the memory footprint.

In a sense this is good news because the amount of storage required for the `[2,2]` block can be nearly cut in half by taking advantage of the fact that it is a triangular matrix.
The [rectangular full packed format](https://netlib.org/lapack/lawnspdf/lawn199.pdf) looks especially promising for this purpose.

In general, for models with scalar random effects for two incompletely crossed grouping factors, the memory footprint depends strongly on the smaller of the number of levels of the grouping factors, less strongly on the larger number, and almost not at all on the number of observations.

### Speed of log-likelihood evaluation

The time required to fit a model to large data sets is dominated by the time required to evaluate the log-likelihood during the optimization of the parameter estimates.
The time for one evaluation is given in the `evtime` column of `sizespeed`.
Also given is the number of evaluations to convergence, `nv`, and the time to fit the model, `fittime`
The reason for considering `evtime` in addition to `fittime` and `nev` is because the `evtime` for one model, relative to other models, is reasonably stable across computers whereas `nev`, and hence, `fittime`, can be affected by seemingly trivial variations in function values resulting from different implementations of low-level calculations, such as the BLAS (Basic Linear Algebra Subroutines).

That is, we can't expect to reproduce `nv` exactly when fitting the same model on different computers or with slightly different versions of software but the pattern in `evtime` with respect to `uc` and `mc` can be expected to reproducible.

As shown in @fig-evtimevsl22 the evaluation time for the objective is predominantly a function of the size of the `[2, 2]` block of `L`.

```{julia}
#| code-fold: true
#| fig-cap: "Evaluation time for the objective (s) versus size of the [2, 2] block of L (GiB)"
#| label: fig-evtimevsl22
ggplot(sizespeed, aes(x=:L22sz, y=:evtime, color=:ucutoff)) +
    geom_point() +
    geom_line() +
    labs(x="Size of [2,2] block of L (GiB)", y="Time for one evaluation of objective (s)")
```

However the middle panel shows that the number of iterations to convergence is highly variable.
Most of these models required between 20 and 25 evaluations but some required almost 50 evaluations.

The derivation of the log-likelihood for linear mixed-effects models is given in @sec-lmmtheory, which provides a rather remarkable result:
the profiled log-likelihood for a linear mixed-effects model can be evaluated from Cholesky factor of a blocked, positive-definite symmetric matrix.

*This page was rendered from git revision {{< git-rev short=true >}}.*