How I'm Watching the MCU, With Stats

Guides

A quick Startpage search turned up eleven guides to the MCU up through Endgame.

Sources ranged from geek blogs to media sites to general news outlets to some Redditor:

• 13th Dimension
• A Place To Hang Your Cape
• BGR
• Collider
• Insider
• Polygon
• Screenrant
• Slashfilm
• Vox
• Vulture (contains link to a later revision; I averaged the two versions)
• zomghax92 on Reddit

Each guide recommends some movies and not others, and they don’t agree on which.

• That is, each one’s picks partly capture some true signal of how essential each movie is.
• But they also contain some random noise due to its own author’s tendencies.

So how do we pull out that common signal to rate each movie statistically?

A first idea would be to average the guides, scoring each movie by how many lists picked it.

• This would treat each guide as containing the same balance of signal and noise.
• But some guides might just be better at choosing essential MCU movies than others!
• It’d be better to weight the average of the votes by the quality of the guides.
• However, that quality isn’t a manifest feature – not a number that can just be read off it.

Instead, a more involved model is needed to jointly rate the movies and the guides.

Item-response theory

Load the data

I copied each guide’s picks into a spreadsheet and loaded them.

• I also captured each movie’s Rotten Tomatoes score, but I’m not using it here.

vengr <- suppressMessages(readr::read_csv(here::here("static", "data", "vengr.csv")))
vengr <- dplyr::select(vengr, -.data$tomatoes)

venga <- tidyr::pivot_longer(vengr, -.data$title, names_to = "item")

# also calculate the naive total votes and set aside
venge <- dplyr::group_by(venga, .data$title)
venge <- suppressMessages(dplyr::summarise(venge, votes =  sum(.data$value)))

Impute the data

Where guides didn’t cast hard up/down votes on movies, I assigned partial votes.

• This mostly affects Vulture, who published two revisions of a tier ranking.
• Values for partial categories were evenly spaced within each guide (e.g. in steps of 1/6).

But the item-response theory model above doesn’t allow for partial votes.

So I used multiple imputation to randomly create several plausible, usable data sets:

1. For each guide, roll a d6 to get a random value \(p\).
2. Treat each guide as voting to watch any movie with \(p/6\) partial votes or above.

Results can then be pooled across imputations by averaging the estimates from each one.

An alternative would be to roll a separate die for each vote, not a single die per guide.

• But it’s weird to think guides suppose movies are required independently of each other.

foo <- dplyr::tibble(i = 1:128)
foo <- dplyr::group_by(foo, .data$i)
foo <- dplyr::mutate(foo, venga = list(venga))
foo <- tidyr::unnest(foo, .data$venga)
foo <- dplyr::group_by(foo, .data$i, .data$item)
foo <- dplyr::mutate(foo, draw = runif(1))
foo <- dplyr::mutate(foo, value = .data$value >= .data$draw)
foo <- dplyr::mutate(foo, value  = as.numeric(.data$value))
foo <- dplyr::select(foo, -.data$draw)

Filter the data

In each imputed data set, I filter out movies on which voting was unanimous.

• Parameters for consensus entries (either way) can’t be estimated.
• They can be simply set to +/- infinity as appropriate, without fitting.

This affected a few movies across all imputations, and two in a fraction of data sets.

• Every guide said to watch The Avengers (along with Infinity War and Endgame).
• All eleven said to skip The Incredible Hulk and Ant-Man.
• Captain America was a consensus watch choice except for one partial vote.
• Iron Man 3 received only partial votes to watch.

foo <- dplyr::group_by(foo, .data$i, .data$title)
foo <- dplyr::mutate(foo,
                     OK    = any(.data$value > 0) & any(.data$value < 1),
                     count = sum(.data$value))

quux <- dplyr::filter(foo, !.data$OK)
quux <- dplyr::group_by(quux, .data$i, .data$title, .data$count, .data$OK)
quux <- dplyr::summarise(quux)
quux <- dplyr::mutate(quux, delta = ifelse(.data$count > 0, Inf, -Inf))

foo <- dplyr::filter(foo, .data$OK)
foo <- dplyr::select(foo, -.data$count, -.data$OK)
foo <- dplyr::group_by(foo, .data$i)
foo <- tidyr::nest(foo)

Fit the model

There are R packages to fit IRT models, but the interfaces confuse me, so I’m rolling my own.

• I wouldn’t use the below for peer-reviewed work without more care, but here it’ll be fine.

The \(\beta\) and \(\delta\) parameters cannot both be fitted in one logistic regression.

• Instead, I alternate between fitting models for the item and respondent parameters.
• Each is fitted conditional on the last estimate of the other (i.e. as if the other were data).
• The alternating steps are repeated until all the parameters jointly converge.

Technical details that I could say more about, but I won’t here:

• For convenience, I fit mixed models that shrink estimates toward zero.
• The parameters are rescaled at each step for purposes of identification.

# Helper functions
huh <- function(x) {dplyr::select(x, tidyselect:::where(is.numeric))}
hah <- function(new, old) {
  y <- huh(new)
  abs((y - huh(old)) / y)
}

#' one-dimensional item-response theory model via lme4::glmer and an EM algorithm
#' @param u `data.frame` with columns `value` (the response), `item` (the item), and `title` (the respondent)
#' @param iter `numeric`. Maximum iterations. Default `128`.
#' @param tol  `numeric`. Relative tolerance for convergence. Default `sqrt(sqrt(.Machine$double.eps))`.
#' @return a list with elements:
#'  * `B` (the final item model) and 
#'  * `D` (the final respondent model) 
#' and attributes:
#'  * `z` (convergence) and 
#'  * `i` (number of iterations).
#' @export
irt <- function(u, iter = 128, tol = sqrt(sqrt(.Machine$double.eps))) {
  z <- 1 # convergence flag
  v <- u # working copy of data (not really needed)
  
  ## initialize alpha and *also* initialize delta *assuming beta = 1 for all items*
  A <- lme4::glmer(value ~ (1 | item) + # alpha
                     (1 | title) + # delta
                     0, # no intercept -- think this is why many fits are singular
                   family = binomial, data = v)
  
  ## extract delta estimates and join to v
  d <- lme4::ranef(A)$title
  d <- dplyr::as_tibble(d, rownames = "title")
  d <- dplyr::rename(d, beta = .data$`(Intercept)`) # yes, this name is as desired
  v <- suppressMessages(dplyr::left_join(v, d))
  
  ## initialize alpha and beta assuming delta equals initial values
  B <- lme4::glmer(value ~ (beta | item) + 0, family = binomial, data = v)
  
  ## extract alpha and beta estimates and join to v
  b <- lme4::ranef(B)$item
  b <- dplyr::as_tibble(b, rownames = "item")
  b <- dplyr::rename(b,
                     alpha = .data$`(Intercept)`,
                     delta = .data$beta)
  v <- suppressMessages(dplyr::left_join(v, b))
  
  for(i in 1:iter) {
    d_1 <- d
    b_1 <- b
  
    # estimate delta given alpha and beta
    D <- lme4::glmer(value ~ (delta + 0 | title) + offset(alpha) + 0, family = binomial, data = v)
  
    ## extract delta estimates and join to v
    d <- lme4::ranef(D)$title
    d <- dplyr::as_tibble(d, rownames = "title")
    d <- dplyr::rename(d, beta = .data$delta)
    d <- dplyr::mutate(d, beta = .data$beta / mean(abs(.data$beta))) # normalize
    v <- dplyr::select(v, -.data$beta)
    v <- suppressMessages(dplyr::left_join(v, d))
    
    ## estimate alpha and beta given delta
    B <- lme4::glmer(value ~ (beta | item) + 0, family = binomial, data = v)
  
    ## extract alpha and beta estimates and join to v
    b <- lme4::ranef(B)$item
    b <- dplyr::as_tibble(b, rownames = "item")
    b <- dplyr::rename(b,
                       alpha = .data$`(Intercept)`,
                       delta = .data$beta)
    v <- dplyr::select(v, -.data$alpha, -.data$delta)
    v <- suppressMessages(dplyr::left_join(v, b))
    
    if(max(hah(d, d_1)) < tol & 
       max(hah(b, b_1)) < tol &
       !lme4::isSingular(D)   &
       !lme4::isSingular(B)   &
       TRUE) {
      z <- 0
      break
    }
  }
  
  out <- list(B = B, D = D)
  attr(out, "z") <- z
  attr(out, "i") <- i
  
  out 
}

system.time(foo <- dplyr::mutate(foo, bar = list(irt(dplyr::first(.data$data)))))
saveRDS(foo, here::here("static", "data", "mcu_foo.rds"))

Pool the results

Finally, I average over the imputations to get final ratings.

• There are standard error estimates, too, and the rules for pooling these are known.
• However, it’s not clear those standard error estimates mean much in the mixed-model setting.

# here's one I made earlier
foo  <- readRDS(here::here("static", "data", "mcu_foo.rds"))
quux <- readRDS(here::here("static", "data", "mcu_quux.rds"))

foo <- dplyr::mutate(foo, bar = list(lapply(dplyr::first(.data$bar), list)))
foo <- dplyr::mutate(foo, bar = list(dplyr::as_tibble(dplyr::first(.data$bar))))
foo <- suppressMessages(dplyr::select(foo, .data$bar))
foo <- tidyr::unnest(foo, .data$bar)

# pool the movie ratings

baz <- suppressMessages(dplyr::select(foo, mod = .data$D))
baz <- dplyr::mutate(baz, mod = list(lme4::ranef(dplyr::first(.data$mod))))
baz <- dplyr::mutate(baz, mod = list(dplyr::as_tibble(dplyr::first(.data$mod))))
baz <- tidyr::unnest(baz, .data$mod)

# prepare to join in the ratings that were infinite
baz <- tidyr::pivot_longer(baz, .data$condval:.data$condsd, "quantity", values_to = "value")
baz <- tidyr::pivot_wider(baz, names_from = .data$grp, values_from = .data$value)
baz <- tidyr::pivot_longer(baz, .data$`Ant-Man and the Wasp`:.data$`Captain America: Civil War`, names_to = "title")

quux <- dplyr::group_by(quux, .data$i, .data$title)
quux <- suppressMessages(dplyr::select(quux, .data$delta))
quux <- tidyr::pivot_longer(quux, .data$delta, names_to = "term", values_to = "degenval")

baz <- suppressMessages(dplyr::left_join(baz, quux))
baz <- dplyr::mutate(baz, value = ifelse(is.na(.data$value), .data$degenval, .data$value))
baz <- suppressMessages(dplyr::select(baz, -.data$degenval))
baz <- tidyr::pivot_wider(baz, names_from = .data$quantity, values_from = .data$value)

baz <- dplyr::group_by(baz, .data$term, .data$title)
baz <- suppressMessages(dplyr::summarise(baz,
                                         n      = sum(is.finite(.data$condval)),
                                         mean   = mean(.data$condval[is.finite(.data$condval)]),
                                         median = median(.data$condval)))
baz <- dplyr::group_by(baz)

# tack on vote totals for comparison
baz <- suppressMessages(dplyr::left_join(baz, venge))

# pool the guide ratings

bar <- suppressMessages(dplyr::select(foo, mod = .data$B))
bar <- dplyr::mutate(bar, mod = list(lme4::ranef(dplyr::first(.data$mod))))
bar <- dplyr::mutate(bar, mod = list(dplyr::as_tibble(dplyr::first(.data$mod))))
bar <- tidyr::unnest(bar, .data$mod)
bar <- dplyr::rename(bar, guide = .data$grp)
bar <- dplyr::mutate(bar, term = ifelse(.data$term == "(Intercept)", "alpha", as.character(.data$term)))
bar <- suppressMessages(dplyr::select(bar, .data$term:.data$condval))
bar <- tidyr::pivot_wider(bar, names_from = .data$term, values_from = .data$condval)
bar <- dplyr::group_by(bar, .data$guide)
bar <- suppressMessages(dplyr::summarise(bar, dplyr::across(tidyselect:::where(is.double), mean)))

Movies

knitr::kable(dplyr::arrange(dplyr::select(baz, -.data$term), dplyr::desc(.data$mean)), digits = 3)

Compared to the simple vote total, there turn out not to be a lot of surprises in the ratings.

• The vote totals almost perfectly line up in the same order as the average ratings.

The ratings do provide clearer gaps between groups of pictures than the vote totals do, though.

• The group with positive ratings seem to be one group. These are the clear essential films.
• Then there seems to be another group of five with near-equal ratings. These are optional.
• Everything after the first five of the negative ratings should be clearly safe to skip.

Median values are included because the reported mean is biased for Civil War and Iron Man 3.

• The mean can only be taken over the data sets where the rating wasn’t infinite.
• But the median can be taken over all 128 imputations, no matter what.

Sources

Sources ranged from geek sites to media blogs to news outlets.

The best at picking out essential movies to watch was… some Redditor.

knitr::kable(dplyr::arrange(bar, dplyr::desc(.data$beta)), digits = 3)

I’m not informed enough to speak of “surprises” here, but a couple facts stood out to me:

• Vulture’s picks surprised me the most compared to other lists (Doctor Strange?).
• Vox dot com managed to hang with the big kids when it came to geek media criticism.
• A Place To Hang Your Cape’s list seemed well thought out, but the model wasn’t impressed.
• At least no one had a negative beta (positively preferring non-essential movies).