Rating children's books with empirical Bayes estimation (1 of 2)

Distribution of star ratings

We already looked at the distribution of the overall rating, but to determine our prior, we’ll need to understand the distribution of ratings for each star level. Specifically, we’re going to look at the percent of a book’s ratings that are at each level.

Why is that important? It turns out that the distribution of ratings at each star level fits well with the Dirichlet distribution, which describes the range of possible probabilities for multiple outcomes. In our case, these outcomes are 1-star rating, 2-star rating, etc.

If you’re familiar with the beta distribution, the Dirichlet is conceptually similar, except that the Dirichlet can apply to more than two outcomes, whereas the beta distribution applies only to situations with two outcomes. (In fact, the Dirichlet is just the more-than-two-outcomes generalization of the beta distribution. A two-outcome Dirichlet is a beta distribution.)

Again, why is this important? If we look at books with enough ratings that their overall rating has “settled”, fitting a Dirichlet distribution will give us a sense of what a good estimate of the percent of a book’s ratings at each star level will be. Visualizing it will help. These are the distributions for books with at least 500 ratings:

rating_pct <- books %>%
  filter(rating_count > 0) %>%
  mutate_at(vars(rating_5:rating_1), ~ . / rating_count) %>%
  # Filter for books with at least 500 reviews
  filter(rating_count > 500) %>%
  pivot_longer(cols = rating_5:rating_1,
               names_to = "rating_level",
               values_to = "pct_of_ratings") %>%
  # Use Unicode symbol for full star (U+2605) and empty star (U+2606) as facet titles
  mutate(rating_num = parse_number(rating_level),
         stars_label = paste0(strrep("\U2605", rating_num),
                              strrep("\u2606", 5 - rating_num)),
         stars_label = fct_reorder(stars_label, rating_num)) %>%
  select(-rating_num)

rating_pct %>%
  ggplot(aes(pct_of_ratings)) +
  geom_histogram(binwidth = 0.025, fill = "#461220", alpha = 0.8) +
  facet_wrap(~ stars_label, nrow = 1) +
  scale_x_continuous(labels = label_percent()) +
  labs(title = "Proportions of ratings at each star level have their own distributions",
       x = "% of ratings",
       y = "Number of books") +
  theme(text = element_text(family = "Bahnschrift"),
        axis.text = element_text(size = 9),
        strip.text = element_text(size = 20, colour = "black"),
        strip.background = element_blank(),
        panel.grid.minor = element_blank())

The percent of a book’s ratings at each star level tends toward certain values. For example, 2-star ratings are likely to constitute 0-15% of ratings, while 4-star ratings are likely to constitute 20-50% of ratings. So, if a book got a new rating, we could look at these distributions to get an idea of how likely it is that this new rating will be for 1, 2, 3, 4, or 5 stars (1%, 4%, 20%, 30%, and 45%, respectively).

We can now use this data to fit a Dirichlet-multinomial distribution that will give us an actual number of ratings at each star level instead of a probability for one rating. This, in effect, gives us “starter ratings” to use for each book. So we will treat a book with zero real ratings as having, say, 40 or 50 ratings that start the book off with an average overall rating. This is how our shrinkage works: a book’s first few ratings will be tempered by our “starter ratings”, but as a book gets more real ratings, those real ratings will have more of an effect on the overall rating.
 

Maximum likelihood estimation of the Dirichlet-Multinomial distribution

We use maximum likelihood estimation (MLE) to estimate the parameters of the Dirichlet-Multinomial distribution using the DirichletMultinomial package. (For an introduction to maximum likelihood estimation, I recommend this StatQuest video by Josh Starmer. I’ll also note that I’m indebted to David Robinson’s already-mentioned book Introduction to Empirical Bayes, which I relied on heavily for the following chunk of code.)

# Create a matrix to feed our MLE of Dirichlet-Multinomial
rating_matrix <- books %>%
  filter(rating_count > 500) %>%
  select(rating_5:rating_1) %>%
  as.matrix()

# Fit a Dirichlet Multinomial distribution
# Input is a matrix of rating counts at each star level
dm_fit <- DirichletMultinomial::dmn(rating_matrix, 1)

# Write a function to tidy DMN object (which dm_fit is)
# (broom's tidy() doesn't work on DMN objects)
# Function from http://varianceexplained.org/r/empirical-bayes-book/
tidy.DMN <- function(x, ...) {
  ret <- as.data.frame(x@fit)
  as_tibble(fix_data_frame(ret, c("conf.low", "estimate", "conf.high")))
}

# Tidy the DMN fit
tidy_dm <- tidy(dm_fit)

# Get parameters into a useful format
par <- tidy_dm %>%
  separate(term, into = c("constant", "rating_stars"), sep = "_", convert = TRUE) %>%
  select(rating_stars,
         prior_ratings = estimate)

# Calculate total prior ratings and mean (used for graphing)
par_total <- sum(par$prior_ratings)
par_mean <- par %>%
  summarise(prior_rating = sum(rating_stars * prior_ratings) / sum(prior_ratings)) %>%
  pull(prior_rating)

The result is a fairly simple dataframe with the number of starter ratings each star level gets:

par %>%
  kable(format = "html") %>%
  kable_styling()

So a book with zero ratings, using our empirical Bayes estimate, will start with 21.29 5-star ratings, 16.50 4-star ratings, etc. (I know it doesn’t make sense to have fractions of ratings, but we’re not talking about real ratings – we can still calculate an overall rating with fractions.)

Overall, our prior gives us about 51 “starter ratings” and an overall rating of 4.07 stars. Again, this will change as a book gets real ratings, but we’ll be less likely to find ourselves in a Pickles and Cake situation where a few good ratings skyrockets a book up to a perfect 5-star rating.
 

Here’s what the Dirichlet fits from MLE look like (light blue curves) compared to the actual distribution of ratings (red histograms):

dirichlet_density <- tidy_dm %>%
  select(rating_level = term, estimate) %>%
  crossing(pct_of_ratings = seq(0, 0.8, by = 0.0125)) %>%
  mutate(density = dbeta(pct_of_ratings, estimate, par_total - estimate)) %>%
  # Use Unicode symbol for full star (U+2605) and empty star (U+2606) as facet titles
  mutate(rating_num = parse_number(rating_level),
         stars_label = paste0(strrep("\U2605", rating_num),
                              strrep("\U2606", 5 - rating_num)),
         stars_label = fct_reorder(stars_label, rating_num)) %>%
  select(-rating_num)

rating_pct %>%
  ggplot(aes(pct_of_ratings)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.025, fill = "#461220", alpha = 0.8) +
  geom_area(data = dirichlet_density, aes(y = density), fill = "#778DA9", alpha = 0.7) +
  facet_wrap(~ stars_label, nrow = 1) +
  scale_x_continuous(labels = label_percent()) +
  labs(title = "Maximum-Likelihood Estimates fit the data fairly well",
       subtitle = "Red histogram = Actual distribution | Blue curve = MLE fit",
       x = "% of ratings",
       y = "Density") +
  theme(text = element_text(family = "Bahnschrift"),
        axis.text = element_text(size = 9),
        strip.text = element_text(size = 20, colour = "black"),
        strip.background = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_blank())

The MLE fit tracks closely to the empirical distribution, so the overall fit is pretty good. One area where the fit isn’t great is the distribution of 5-star ratings: our MLE fit has under-estimated the variance. You can tell this by the MLE blue curve being narrower and taller than the empirical distribution. Even though it’s not perfect, I’m making the call that it’s good enough to use for our shrinkage purposes.
 

Incorporating the prior into ratings

We can finally calculate the shrunk ratings. We just need need to add the prior “starter ratings” to each book’s existing real ratings.

# Calculate empirical Bayes rating using our prior
books_eb <- books %>%
  pivot_longer(rating_5:rating_1,
               names_to = c("rating_stars"),
               names_pattern = "rating_(.*)",
               names_transform = list(rating_stars = as.integer),
               values_to = "ratings") %>%
  left_join(par, by = "rating_stars") %>%
  group_by(isbn, title, author, rating_count) %>%
  # Calculate original rating from scratch because rating field truncates after two decimal places
  summarise(rating_calc = sum(rating_stars * ratings) / sum(ratings),
            rating_eb = sum(rating_stars * (ratings + prior_ratings)) / sum(ratings + prior_ratings)) %>%
  ungroup()

And we end up with an adjusted rating for each book (rating_eb). Here are the top 10 books, according to our empirical Bayes estimate:

books_eb %>%
  select(-isbn, -rating_calc) %>%
  top_n(10, wt = rating_eb) %>%
  arrange(desc(rating_eb)) %>%
  kable(format = "html") %>%
  kable_styling()

Immediately, I’m more confident in these ratings than the originals, which had a bunch of books with 2 or 3 ratings that just happened to be 5 stars. Visually, this is what shrinkage did:

# Format data properly for animation
books_shrunk <- books_eb %>%
  filter(rating_count > 0) %>%
  select(isbn:author, rating_count, rating_calc, rating_eb) %>%
  pivot_longer(rating_calc:rating_eb, names_to = "rating_type", values_to = "rating") %>%
  mutate(rating_type = ifelse(rating_type == "rating_calc", "Original Rating", "Empirical Bayes Rating"))

# Create base plot
p <- books_shrunk %>%
  filter(rating_count <= 200) %>%
  ggplot(aes(rating_count, rating, col = rating_count)) +
  geom_point(alpha = 0.4, size = 0.5) +
  geom_hline(yintercept = par_mean,
             lty = 2,
             size = 0.8,
             col = "red") +
  scale_colour_fish(option = "Ostracion_whitleyi") +
  scale_x_continuous(breaks = seq(0, 200, by = 25)) +
  labs(subtitle = paste0("Dashed line shows Bayesian prior for a book with 0 ratings",
                        "\n",
                        "(Only books with 200 ratings or less shown)"),
       x = "Number of ratings",
       y = "Rating") +
  theme(legend.position = "none",
        text = element_text(family = "Bahnschrift", size = 6),
        plot.title = element_text(size = 8),
        panel.grid.minor.x = element_blank())

# Set animation parameters
anim <- p +
  transition_states(rating_type,
                    transition_length = 1.5,
                    state_length = 1.5) +
  ease_aes("cubic-in-out") +
  ggtitle("{closest_state}")

# Create animation
anim

The books with few ratings – less than 25 – were impact a lot by incorporating a Bayesian prior. But as we get more rating, the prior has less of an effect. You can imagine a book that has 1,000 or 10,000 ratings: our prior will barely move its score.

Before wrapping up, we should note that still need to be careful when interpreting adjusted ratings. They are still susceptible to other forms of statistical bias. For example, the top-rated book is 100 Fun Bible Facts by Ginger Baum. I’m sure it’s an excellent book, but not everyone might consider it excellent. People who took the time to read and rate this book are likely pre-disposed to enjoy books with biblical themes, making this a likely illustration of self-selection bias.