Is the EPL getting more unequal?

I recently heard that Manchester City were so far ahead in the English Premier League (EPL) that the race for first was basically over, even though they were still about 6-7 more games to go (out of a total of 38 games). At the other end of the table, I heard that Sheffield United were so far behind that they are all but certain to be relegated. This got me wondering: has the EPL become more unequal over time? I guess there are several ways to interpret what “unequal” means, and this post is my attempt at quantifying it.

The data for this post is available here, while all the code is available in one place here. (I also have a previous blog post from 2018 on the histogram of points scored in the EPL over time.)

Let’s start by looking at a density plot of the points scored by the 20 teams and how that has changed over the last 12 years (note that the latest year is at the top of the chart):

library(tidyverse)
df <- read_csv("../data/epl_standings.csv")

theme_set(theme_bw())

# joy plot
library(ggridges)
ggplot(df, aes(x = Points, y = factor(Season))) +
    geom_density_ridges(scale = 3) +
    labs(y = "Season")

It looks like the spread for the more recent seasons is greater than that in earlier seasons, evidenced by the bulk of the distribution widening as we go from the bottom to the top of the chart.

Another way to view the data is to draw a boxplot for each season. In this view, it’s harder to see the spread that we saw in the joy plot.

# boxplot
ggplot(df, aes(x = Season, y = Points, group = Season)) +
    geom_boxplot() +
    labs(title = "Distribution of points by season")

One possible way to define “unequal” is to look at the difference between the number of points the top team earned vs. that for the bottom team: the larger the difference, the more unequal the EPL is becoming. The plot below shows the maximum and minimum number of points earned by a team across the years, as well as the difference between the two. We also show the linear regression lines for each of these values. With the upward slopes, it looks like the gap between the best and the worst is certainly increasing.

# plot of maximum and minimum
df %>% group_by(Season) %>%
    summarize(max = max(Points), min = min(Points)) %>%
    mutate(range = max - min) %>%
    pivot_longer(max:range) %>%
    ggplot(aes(x = Season, y = value, col = name)) +
    geom_line() +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    scale_x_continuous(breaks = seq(min(df$Season), max(df$Season), by = 2)) +
    labs(title = "Max/min points by season", y = "Points") +
    theme(legend.title = element_blank(), legend.position = "bottom")

The problem with using range is that it only tracks the difference between the best and the worst teams while ignoring all the teams in the middle. A measure that takes the middle into account is variance. Note that the number of points that a team can score lies in the range 0 to . For random variables bounded in the interval , the smallest possible variance is 0 (all teams scoring the same number of points), while the maximum possible variance is (half the teams scoring 114 points, half the teams scoring 0 points). Based on the configurations attaining the extremes, it seems reasonable to interpret the scores having a higher variance as the league being more unequal.

Here is a plot of point variance over time along with the linear regression line:

df %>% group_by(Season) %>%
    summarize(var = var(Points)) %>%
    ggplot(aes(x = Season, y = var)) +
    geom_line() +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    scale_x_continuous(breaks = seq(min(df$Season), max(df$Season), by = 2)) +
    labs(title = "Variance of point distribution by season", y = "Variance")

There’s quite a lot of fluctuation from year to year, but there seems to be an increasing trend. However, notice that the -axis only goes from 250 to 450, while our min and max variances were 0 and 3249. Perhaps it’s better to have the -axis go from 0 to 3249 to have a better sense of scale?

Before doing that, note that it’s actually not possible for half the teams to score maximum points. In the EPL, every team plays every other team exactly twice. Hence, only one team can have maximum points, since it means that everyone else loses to this team. I don’t have a proof for this, but I believe maximum variance happens when team 1 beats everyone, team 2 beats everyone except team 1, team 3 beats everyone except teams 1 and 2, and so on. With this configuration, the variance attained is just 1260.

Here’s the same line plot but with the -axis going from 0 to 1260. With this scale, the change in variance over time looks pretty flat.

max_var_dist <- seq(0, 38 * 3, by = 6)
max_var <- var(max_var_dist)
df %>% group_by(Season) %>%
    summarize(var = var(Points)) %>%
    ggplot(aes(x = Season, y = var)) +
    geom_line() +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    geom_hline(yintercept = c(0, max_var), col = "red", linetype = "dashed") +
    scale_x_continuous(breaks = seq(min(df$Season), max(df$Season), by = 2)) +
    labs(title = "Variance of point distribution by season", y = "Variance")

Which is the correct interpretation? Is the change in variance over time important enough that we can say the EPL has become more unequal, or is it essentially the same over time? This is where I think domain expertise comes into play. 1260 is a theoretical maximum for the variance, but my guess is that the layperson looking at two tables, one with variance 300 and one with variance 900, would be able to tell them apart and say that the latter is more unequal. Can the layperson really tell the difference between variances of 250 and 450? I would generate several tables having these variances and test if people could tell them apart.

Finally, the Gini coefficient is one other measure of inequality, with 0 being the most equal and 1 being the most unequal.

# plot of Gini
library(DescTools)
df %>% group_by(Season) %>%
    summarize(gini = DescTools::Gini(Points)) %>%
    ggplot(aes(x = Season, y = gini)) +
    geom_line() +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    geom_hline(yintercept = c(0, 1), col = "red", linetype = "dashed") +
    scale_x_continuous(breaks = seq(min(df$Season), max(df$Season), by = 2)) +
    labs(title = "Gini coefficient of point dist by season", 
         y = "Gini coefficient")

Here are the plots with different scales for the -axis:

As with variance, the different scales give very different interpretations. It will require some research to figure out if a change of Gini coefficient from 0.17 to 0.22 is perceptible to the viewer.

What is a proper scoring rule?

What is a proper scoring rule?

In the realm of forecasting, a scoring rule is a way to measure how good a probabilistic forecast is.

Mathematically, if is the set of possible outcomes, then a probabilistic forecast is a probability distribution on (i.e. how likely each possible outcome is). A scoring rule (or a score) is a function that takes a probability distribution on (denoted by ) and an outcome as input and returns a real-valued number as output. If is the probabilistic forecast and is the actual outcome, then is interpreted as the reward (loss resp.) to the forecaster if the score is positively-oriented (negatively-oriented resp.).

Assume that the true probability distribution is denoted by . For each probabilistic forecast , we can define the expected score as

A positively oriented scoring rule is said to be proper if for all probability distributions and , we have

In other words, for a proper score, the forecaster maximizes the expected reward if he/she forecasts the true distribution. A strictly proper score is a score such that equality above is achieved uniquely at . To maximize a strictly proper score, a forecaster has every incentive to give an “honest” forecast and has no incentive to “hedge”.

Why would you ever use an improper scoring rule?

Since score propriety (i.e. a score being proper) seems like such a basic requirement, why would anyone use an improper scoring rule? It turns out that there are some scores with nice properties but are not improper.

For example, imagine that we want to compare scores against some baseline forecaster. Given some score , we can convert it to a skill score through normalization:

where is the probabilistic forecast which the baseline makes. If is negatively-oriented (smaller is better) and if the scores it produces are always non-negative, then the associated skill score is positively-oriented and in the range .

Skill scores seem reasonable as they give us a fixed upper bound to aspire to. However, in general skill scores are improper. (See Reference 2 for a proof, and Section 2.3 of Reference 1 for another normalization scheme.)

Two other examples of improper scores that have been used are the naive linear score and mean squared error (MSE) (see Reference 3 for details). Note that mean squared error here is not the usual MSE we use for point forecasts: there is a more general definition for probabilistic forecasts.

Some examples of proper scoring rules

Reference 1 contains a number of examples of proper scoring rules. First, assume that the sample space for the response is categorical (with loss of generality, let it be ), and let the probabilistic forecast be represented by . Here are 4 proper scoring rules for categorical variables (only 1-3 are strictly proper):

1. The Brier score: , where . (I wrote a previous post on the Brier score here.)
2. The spherical score: For some parameter , the (generalized) spherical score is defined as . The traditional spherical score is the special case .
3. The logarithmic score: .
4. The zero-one score: Let be the set of modes of the probabilistic forecast. Then the zero-one score is defined as .

There are similar examples for continuous responses. For simplicity, assume that the probabilistic forecast has a density w.r.t. Lebesgue measure. (See Reference 1 for a more mathematically rigorous description.) Define

1. The quadratic score: .
2. The pseudospherical score: For , .
3. The logarithmic score: . This can be viewed as the limiting case of the pseudospherical score as . It is widely used, and is also known as the predictive deviance or ignorance score.

There are options for proper scoring rules if the probabilistic forecast does not have a density. Every probabilistic forecast will have a cumulative distribution function , and we can construct scores from that. The continuous ranked probability score (CRPS) is one such score, defined as

The CRPS corresponds to the integral of the Brier scores for the associated binary probability forecasts at all real-valued thresholds.

References:

1. Gneiting, T., and Raftery, A. E. (2007). Strictly proper scoring rules, prediction and estimation.
2. Murphy, A. H. (1973). Hedging and skill scores for probability forecasts.

3. Bröcker, J., and Smith, L. A. (2007). Scoring probabilistic forecasts: The importance of being proper.

Estimating pi with Monte Carlo simulation

Happy Pi Day! I don’t encounter very much in my area of statistics, so this post might seem a little forced… In this post, I’m going to show one way to estimate .

The starting point is the integral identity

There are two ways to see why this identity is true. The first is that the integral is simply computing the area of a quarter-circle with unit radius. The second is by explicitly evaluating the integral:

If , then the integral identity means that

Hence, if we take i.i.d. draws from the uniform distribution on , it is reasonable to expect that

The code below shows how well this estimation procedure does for one run as the sample size goes from 1 to 200:

set.seed(1)
N <- 200
x <- runif(N)
samples <- 4 * sqrt(1 - x^2)
estimates <- cumsum(samples) / 1:N
plot(1:N, estimates, type = "l",
     xlab = "Sample size", ylab = "Estimate of pi",
     main = "Estimates of pi vs. sample size")
abline(h = pi, col ="red", lty = 2)

The next plot shows the relative error on the y-axis instead (the red dotted line represents 1% relative error):

rel_error <- abs(estimates - pi) / pi * 100
plot(1:N, rel_error, type = "l",
xlab = "Sample size", ylab = "Relative error (%)",
main = "Relative error vs. sample size")
abline(h = 0)
abline(h = 1, col = "red", lty = 2)

Log-likelihood of probit regression is globally concave

Assume we are in the supervised learning setting with observations, where observation consists of the response and features . A generalized linear model (GLM) consists of 3 components:

1. A random component, or a family of distributions indexed by (usually an exponential family), such that ,
2. A systematic component , and
3. A link function such that .

(See this previous post for more details of the components of a GLM.) The user gets to define the family of distributions and the link function , and is the parameter to be determined by maximum likelihood estimation.

For one-dimensional exponential families with the canonical link function, it is known that the log-likelihood of the GLM is globally concave in (see, for example, Reference 1). Hence, the MLE can be found using methods such as gradient descent or coordinate descent. When non-canonical links are used, the GLM’s log-likelihood is no longer guaranteed to be concave in . However, in some situations we can stilll show that the log-likelihood is concave in . In this post, we show that the log-likelihood for probit regression is concave in .

In the probit regression model, , where is the cumulative distribution function (CDF) of the standard normal distribution. The responses are binary with . The likelihood function is

and the log-likelihood function is

To show that is concave in , we make two reductions:

1. Since the sum of concave functions is concave, it is enough to show that is concave.
2. Since composition with an affine function preserves concavity, it is enough to show that is concave in . (Here, .)

From here, we can show that is concave by showing that its second derivative is negative: for all . Since can only take on the values of 0 and 1, we can consider those cases separately.

Let denote the probability density function of the standard normal distribution. Recall that . When ,

since for all (see this previous post for a proof).

When ,

To show concavity of , it remains to show that for all . But this is true: see this previous post for a proof.

References:

1. Chapter 9: Generalized Linear Models.

A lower bound on x Phi(x) + phi(x)

I recently came across this lower bound that I thought was interesting and not immediately obvious. Let and be the probability density function (PDF) and the cumulative distribution function (CDF) of the standard normal distribution, i.e.

Then the following inequality holds:

The plot below illustrates this inequality with the function on the LHS and RHS being in red and black respectively:

x <- seq(-30, 30) / 10
plot(x, x * pnorm(x) + dnorm(x), type = 'l', col = "red",
     ylab = "f(x)")
lines(x, pmax(0, x), lty = 2)


The derivation of the inequality below is mine; there may be faster ways to do it. A key formula we will use in the proof is that for the derivative of : .
Let . First, we show that  for all . Since, for all ,

the function  is strictly increasing on the entire real line. Hence, using L’Hôpital’s rule,

Next, we show that  for all . Note that we only have to prove this for , since  for . Letting , we see that proving  for  is equivalent to proving  for . Since

for all , the function  is strictly decreasing for . Thus, for this range of  values,

as required. (We used L’Hôpital’s rule again in the sequence of equalities above.)
					

tidytuesdayR package for easy downloading on tidy tuesday data

Note to self: the tidytuesdayR R package (available on CRAN) makes it easy to load Tidy Tuesday datasets into R. See the README at the package’s github repository for a short vignette on how to use the package.

What is Gower’s distance?

Let’s say we have two observations and , and we want a measure of how similar (or dissimilar) they are to each other. If each of the entries in these two vectors are quantitative (i.e. take on real-number values), then we can use distance measures as a dissimilarity measure (e.g. Euclidean distance, Manhattan distance). But what if the features that make up aren’t all numeric? What if some of them are categorical or binary?

Gower’s distance, introduced in Gower (1971) (Reference 1), is a general similarity measure that can be used in this setting. For each feature , we define a score . If and are close to each other along feature , then the score is close to 1. Conversely, if they are far apart along feature the score is close to 0.

How the score is computed depends on the type of feature (to be described shortly). A quantity is also computed: if and can be compared along feature , then . If and cannot be compared along feature (because, for example, of missing values), is set to zero. Gower’s distance is just the average of the (known) scores:

Now, let’s talk about how to compute the scores for each type of feature. Gower describes 3 different types of features:

1. Quantitative variables (numeric variables): , where is the range of the feature , either in the population or in the sample.
2. Qualitative variables (categorical variables): .
3. Dichotomous variables: This is something we don’t talk about so much nowadays. A dichotomous variable is one where the feature is either present or absent, and “whose absence in both of a pair of individuals is not taken as a match” (emphasis mine). For such variables, and are determined by the table below (taken from Reference 1).

Finally, why is this called a distance? If and have similarity , then we can think of the distance between them as . For this to be a valid distance, it must satisfy the triangle inequality, i.e. for any 3 observations , and , their similarities must satisfy

Reference 1 proves that if there are no missing values (i.e. for all , and ), then the distance defined above satisfies the triangle inequality.

(Credit: I learned of Gower’s distance from Max Kuhn‘s book Feature Engineering and Selection: A Practical Approach for Predictive Models (Section 8.5).)

References:

1. Gower, J. C. (1971). A general coefficient of similarity and some of Its properties.

The Box-Cox and Yeo-Johnson transformations for continuous variables

The Box-Cox and Yeo-Johnson transformations are two different ways to transform a continuous (numeric) variable so that the resulting variable looks more normally distributed. They are often used in feature engineering to reduce skew in the raw variables.

Box-Cox transformation

George Box and David Cox proposed the Box-Cox transformation in Reference 1. The transformation is really a family of transformations indexed by a parameter :

This is also known as a power transformation as the variable is raised to a particular power. (Note that under translation and scaling, the Box-Cox transformation is really the same as the transformation ; the translation and scaling above is set so that the transformation is continuous at , making it easier for theoretical analysis.)

The transformation above only works for . If can be negative but we always have for some , then we can use

The key question is how to choose the best value of . Box and Cox propose using maximum likelihood estimation to do so.

Yeo-Johnson transformation

Yeo and Johnson (2000) (Reference 2) note that the tweak above only works when is bounded from below, and also that standard asymptotic results of maximum likelihood theory may not apply. They propose the following transformation:

The motivation for this transformation is rooted in the concept of relative skewness introduced by van Zwet (1964) that I won’t go into here. Some nice properties:

• is concave in for and convex for .
• The constant shift of makes is such that the transformed value will always have the same sign as the original value.
• The constant shift of also allows to become the identity transformation when .
• The new transformations on the positive line are equivalent to the Box-Cox transformation for (after accounting for the constant shift), so the Yeo-Johnson transformation can be viewed as a generalization of the Box-Cox transformation.

As with the Box-Cox transformation, the value of is chosen via maximum likelihood estimation.

(Credit: I first heard of the Box-Cox transformation in Art Owen‘s STATS 305A class, and I learned of the Yeo-Johnson transformation from Max Kuhn‘s book Feature Engineering and Selection: A Practical Approach for Predictive Models.)

References:

1. Box, G. E. P., and Cox, D. R. (1964). An analysis of transformations.
2. Yeo, I. K., and Johnson, R. A. (2000). A new family of power transformations to improve normality or symmetry.
3. var Zwet, W. R. (1964). Convex transformations of random variables.

What is the jackknife+?

This post describes the jackknife+, a new method developed by Barber et al. (2021) (Reference 1) for constructing prediction intervals with provable coverage guarantees.

The data and the goal

Assume we have i.i.d. training data and a new test point drawn from the same distribution. We want to fit a regression model using the training data, then produce a prediction interval for the test point, i.e. an interval around the prediction that is likely to contain the true response . Concretely, we want a prediction interval (as a function of the training points) such that

where is a pre-specified level of desired coverage. The probability is taken with respect to a new test point as well as with respect to the training data.

The jackknife and the jackknife+

Before introducing the methods, let me introduce some notation that will help us express the methods concisely.

• For , let denote the regression model fit on the training data with the th observation removed.
• For , let denote the leave-one-out residual for observation .
• For any values , let and denote the and quantiles of the empirical distribution of these values respectively. (Note: The paper has a more precise, technical definition of these quantiles.)

The original jackknife prediction interval is defined as

The newly proposed jackknife+ prediction interval is defined as

You might have to squint a little to see the difference between the two intervals. The jackknife centers the prediction interval on the predicted value fitted on the full training data, while the jackknife+ uses leave-one-out-points instead. The figure below is a visual example of how the construction of the two types of prediction intervals might differ.

The beauty of this little change is that it comes with provable coverage guarantees, without any assumptions! Well, OK, extremely minor assumptions. We assume that the method used to obtain the regression model is invariant to the ordering of the data, and that the sample size . No assumptions on the data generating distribution and on the regression method otherwise. Here’s the guarantee (Theorem 1):

Theorem: The jackknife+ prediction interval satisfies

The alert reader will notice that the RHS is and not as we would like. The authors address this in 3 different ways:

1. is a lower bound on the coverage. Empirically it looks like coverage is often very close to the target level .
2. They propose a version of the jackknife+, which they call jackknife-minimax, that has an (see Section 2.2 of Reference 1 for details). However, the authors were of the opinion that this method is generally too conservative.
3. While there are pathological examples that cause the jackknife+ to give coverage very close to (Theorem 2), this does not happen often. In Section 5 of the paper, the authors introduce the notion of “in-sample stability” and “out-of-sample stability”, and give coverage results that are much closer to under some stability conditions (Theorem 5).

Extension to cross-validation (CV+)

The computation of the jackknife and jackknife+ intervals can be inefficient because we need to compute for . This amounts to refitting the model times. The authors develop an extension of the jackknife+ that works in the -fold cross-validation (CV) setting. Here, the model refitting only happens times.

Assume that the training data is splt into disjoint subsets each of size . Define the following quantities:

• For , let be the regression model fit on the training data with removed.
• For , let , where is the subset that contains observation .

The CV+ prediction interval is defined as

(Note that the jackknife+ prediction interval is simply the CV+ prediction interval with .) The coverage guarantee for CV+ is as follows (Theorem 5):

Theorem: The CV+ prediction interval satisfies

Comparison with other methods

Reference 1 compares the jackknife+ against other methods in terms of theoretical guarantees and computational efficiency. I’ve reproduced the tables here, with the newly proposed methods highlighted in green.

Note: All figures in this post are taken from Reference 1.

References:

1. Barber, R. F. et al. (2021). Predictive inference with the jackknife+.

What are spaghetti plots and lasagna plots?

I recently heard of lasagna plots in a talk and I hadn’t heard of them before, so I looked it up. It turns out that the lasagna plot and its close cousin, the spaghetti plot, are plots that you’ve probably come across before. They are just technical terms that are used in certain domains with certain connotations.

According to Reference 1,

[the spaghetti plot] plots a subject’s values for the repeated outcome measure (vertical axis) versus time (horizontal axis) and connects the dots chronologically, using lines of uniform color and intensity for each subject.

Sounds like a line plot to me! A line plot with one line per subject and with time as the x-axis. The figure below demonstrates what a spaghetti plot is and what could go wrong. The left panel shows the time series for a single individual (one strand of spaghetti). Unfortunately, the entire dataset consists of 118 subjects, so plotting all the strands of spaghetti results in a hot mess.

Swihart et al. (2010) (Reference 1) introduce the lasagna plot as a potential remedy for the spaghetti plot’s shortcomings. The lasagna plot is a heatmap with some special connotations. Instead of one line per subject in the spaghetti plot, each subject gets one horizontal slice. Instead of depicting values using the y-axis, they are presented through the color of each part of the slice. Here is an example of a lasagna plot with 4 subjects and 6 timepoints:

Here is the lasagna plot for the same data which the spaghetti plot made a mess of. Since there are 118 subjects, each horizontal slice is very thin.

Certainly more palatable than the spaghetti plots, with colors to make it look like a real lasagna too! The panel on the left shows the plot with subjects labeled according to their position in the dataset. The right panel shows the plot with subjects ordered according to disease status, and within each disease status, ordered by data of recording. Some trends are obvious now: the top half tend to have lower values, and there seems to be a systematic pattern of missing data for some subjects.

Reference 1 goes into more detail about the possible ways to sort the rows and columns of the lasagna plot, as well as pointers on how to choose a suitable color palette.

Note: All figures were taken from Reference 1.

References:

1. Swihart, B. J. (2010). Lasagna plots: A saucy alternative to spaghetti plots.