Statilisticles

In my line of work—and in my social-media universe—I seem to see a lot of ranking of categories: most frequent words in a text corpus, most frequent genres in a bibliographic listing, most frequent nationalities in a group of authors… I’ve often wondered how much uncertainty lies behind rankings but had no idea how to quantify this uncertainty. I caught a link (via @coulmont) to a short essay on the subject: Michael Höhle, “Rank Uncertainty: Why the ‘Most Popular’ Baby Names Might Not Be the Most Popular” (available in OA preprint), the basics of which I was able to follow, thanks to its clear exposition and accompanying R code. Höhle explains that even when one has a full census of the names of all the babies born in a year, there is still uncertainty associated with ranking the names: the true propensity in the population to give a particular name is not fully revealed by the actual total. Höhle suggests we can begin by modeling the uncertainty arising from the fact that (1) even if the birth rate is fixed, there are “chance” variations in the number of births and (2) even if the probabilities of names are fixed, each draw from the bag of names has some randomness (little Jacob might well have been William). The rank uncertainty can then be found by simulating a year’s worth of names many times and finding popularity ranks in each simulation.

This seemed like a useful method to be able to use the next time I look at a “top 10 list” for evidence—a way to replace listicles with statilisticles. And the data Höhle uses, the British Office for National Statistics lists of baby names and sexes for 2015, has recently been made available in a convenient R package by Thomas Leeper. So I decided to try to understand the point by reproducing the uncertainties in Höhle’s essay. Here’s my code:

library(tidyverse)
library(doMC)
library(foreach)
registerDoMC(cores=4)

# data: US and UK
library(babynames)
library(ukbabynames)

library(knitr) # for kable

# leave nothing to pseudo-chance!
set.seed(20170623)

The simulation is a simple hierarchical model: the number of individuals is first drawn from a Poisson distribution with mean N, and then each individual’s category membership is drawn from a multinomial with a vector of category probabilities p.

simulate <- function (N, p) rpois(1, N) %>%
    rmultinom(1, size=., prob=p)

We’ll repeat the simulation many times using doMC::times and then use the percentile method to obtain bootstrapped confidence intervals for the ranks. Here is my function for this, closely modeled on Höhle’s, but with some effort to make this something you could use for a non-baby-name application. One subtlety is that in the case of names, the ranks of interest are actually ranks within sexes, so we have to keep track of this grouping factor. I’ll do that by making the ranking function a parameter r and passing in a function that respects this grouping (the default r parameter does not).¹

bootstrap_ranks <- function (x, # data frame with counts in rows
                             var, # name of column with counts
                             r=function (x) min_rank(desc(x[[var]])),
                             rank_threshold=Inf,
                             level=0.95, # confidence level
                             n_reps=1000 # no. of simulations
                             ) {
    # ignore grouping on x for these parameters
    x <- ungroup(x)
    N <- sum(x[[var]])
    p <- x[[var]] / N

    # data ranks
    ranks <- r(x)
    # rank mask
    keep <- ranks <= rank_threshold

    # parametric bootstrap: sample ranks (then apply mask)
    ranks_boot <- times(n_reps) %dopar% {
        x[[var]] <- simulate(N, p)
        list(r(x)[keep])
    }

    # tack data on as first replicate
    ranks_boot <- c(list(ranks[keep]), ranks_boot)

    # transpose to a list of ranks for each name
    ranks_boot <- ranks_boot %>%
        transpose() %>%
        map(flatten_int)

    # return rows from with simulated ranks and CIs
    x %>% filter(keep) %>%
        mutate(rank_samples=ranks_boot) %>%
        # add percentile CIs; type=3 following Höhle
        mutate(ci=map(rank_samples, quantile,
                      prob=c((1 - level) / 2, 1 - (1 - level) / 2), type=3))
}

Now let’s try to reproduce Höhle’s table with UK names. To set this up we need a gender-specific ranking function and a function for extracting what he calls the “uncertainty-corrected rank,” UCR, which is the lower bound of the bootstrap confidence interval for the rank—that is to say, the upper bound of our estimate for the name’s position in the league table. If two names have the same UCR we can’t reject the hypothesis that their ranks are the same at the given confidence level.²

# function for ranking names within sexes: returns a vector
name_rank <- . %>%
    group_by(sex) %>% 
    mutate(rank=min_rank(desc(n))) %>%
    pull(rank) # dplyr 0.7

name_bootstraps <- . %>%
    filter(year == 2015) %>%
    select(name, sex, n) %>%
    group_by(sex) %>%
    mutate(p=n / sum(n)) %>%
    # simulate
    bootstrap_ranks(var="n", r=name_rank,
                    rank_threshold=16, level=0.95, n_reps=10000)

K_actual <- 697788 # according to Höhle

ucr_table <- . %>%
    # extract CI lower bound
    mutate(ucr=map_int(ci, 1)) %>%
    # format table
    arrange(sex, ucr) %>%
    mutate(p=scales::percent(p), gsp=scales::percent(n / K_actual)) %>%
    filter(ucr <= 14) %>%
    select(sex, name, UCR=ucr, `% of sex`=p, `% of all`=gsp) %>%
    kable()

The simulation with 10000 replicates takes a couple of minutes on my computer.

ukbabynames %>% name_bootstraps() %>% ucr_table()

Now let’s look at analogous results for the US in 2015, making use of the babynames R package:

us_ranks <- babynames %>% name_bootstraps()
us_ranks %>% ucr_table()

rank_plot <- function (x, s) {
    x <- x %>% filter(sex==s) %>% 
        mutate(pop_rank=map_int(rank_samples, 1),
               name=reorder(name, pop_rank))

    csets <- x %>%
        select(name, rank_samples) %>% 
        unnest(rank_samples) %>%
        group_by(name) %>%
        count(rank_samples)

    cis <- x %>%
        select(name, pop_rank, ci) %>%
        # integers, you know?
        mutate(lb=map_int(ci, 1) - 0.5, ub=map_int(ci, 2) + 0.5)

    ggplot(csets) + 
        geom_raster(aes(name, rank_samples, fill=n, alpha=n)) +
        scale_fill_gradient(low="white", high="black", guide="none") +
        # alpha scaling because my website doesn't have a white background
        # and I'm fussy
        scale_alpha_continuous(range=c(0, 1), guide="none") +
        geom_errorbar(data=cis, aes(name, ymin=lb, ymax=ub), color="blue") +
        geom_point(data=cis, aes(name, pop_rank), color="blue") +
        theme(axis.text.x=element_text(angle=45, hjust=1),
              panel.grid=element_blank()) +
        xlab(NULL) + ylab("rank")
}

rank_plot(us_ranks, "F")

rank_plot(us_ranks, "M")

As always when I try to learn something with statistics, what I learn is that I know less than I thought I knew. Use my code, and tell your friends!

Edited, 6/23/17: Set the random seed (reproducibility! the UCR of “Grace” in the UK actually did fluctuate). To get exactly the graphics here you’d also need the custom R markdown format I created to work with Hugo: hugormd. Also replaced an abuse of the term “confidence sets”; I have no confidence when talking about confidence.