Regression Modeling With Proportion Data (Part 1)

Setup

library("conflicted")
suppressPackageStartupMessages(library("tidyverse"))
conflict_prefer("filter", "dplyr")
#> [conflicted] Will prefer dplyr::filter over any other package
library("betareg")
library("broom")
library("kableExtra")
theme_set(theme_light(base_size = 12))
theme_update(panel.grid.minor = element_blank())

# Function to transform attendance rates for beta regression because values of 0
# and/or 1 may occur
transform_perc <- function(percentage_vec) {
    # See Cribari-Neto & Zeileis (2010)
    (percentage_vec * (length(percentage_vec) - 1) + 0.5)/length(percentage_vec)
}

Selected Variables

Dependent variable:

• Beta regression: Attendance rate; values were transformed to the interval (0, 1) using transform_perc()
• Quasi-binomial regression: Attendance rate in the interval [0, 1]
• Linear regression: Attendance (i.e., count)
• In all cases, entries where the attendance was larger than the capacity were replaced with the maximum capacity.

Independent variables:

• Home team: Different teams have different average attendance rates and/or stadium capacities. For the beta model, home team was also used to predict the dispersion parameter phi.
• Away team: Generally, more successful away teams like Kiel attract more spectators.
• Year/season: Used to capture potential differences across the three seasons.
• Weekday: Matches on weekends probably attract more spectators.
• All variables are categorical variables and appropriate coding schemes (i.e., dummy or effect coding where chosen).

Illustrative Plot of Estimates

The model estimates can also be illustrated graphically. The parameters mean (mu) and precision (phi) from above correspond to a beta distribution, which is specific to each home team in the fitted model. To calculate the densities of the beta distribution using dbeta(), I first have to transform the estimates from the model.

# Function to calculate a and b parameter of beta distribution from mu and phi
# used in betareg with link function logit for mu and log for phi
calc_ab <- function(mu, phi) {
    mu <- plogis(mu)
    phi <- exp(phi)
    b <- phi - mu*phi
    a <- phi - b
    return(cbind(a = a, b = b)[, 1:2])
}

# Calculate predicted beta distributions for Flensburg and Lemgo
preds1 <- 
    # Estimates of mu and phi from above
    calc_ab(2.08 + c(1.42, -.42), 2.4 + c(.49, .03)) %>%
    asplit(1) %>%
    setNames(c("SG Flensburg-Handewitt", "TBV Lemgo Lippe")) %>%
    # Calculate densities via dbeta()
    lapply(function(y) data.frame(x = seq(0, 1, .001),
                                  y = dbeta(seq(0, 1, .001), y[1], y[2]))) %>%
    enframe("home_team") %>% 
    unnest(value)

# Print some values
slice(preds1, seq(701, by = 100, length.out = 4))
#> # A tibble: 4 x 3
#>   home_team                  x        y
#>   <chr>                  <dbl>    <dbl>
#> 1 SG Flensburg-Handewitt   0.7   0.0133
#> 2 SG Flensburg-Handewitt   0.8   0.145 
#> 3 SG Flensburg-Handewitt   0.9   1.40  
#> 4 SG Flensburg-Handewitt   1   Inf

Below, I plotted the model-implied beta distributions for the two teams, which are illustrated using the green lines. Furthermore, the plot shows the predicted means (green point on the x-axis), the raw densities calculated by geom_density() across the 34 matches per team, as well as the observed attendance rate for each match (placed below the densities by geom_rug()).

In summary, the plot nicely illustrates the estimates from above, namely, very high and homogeneous attendance rates for Flensburg and lower and more varied attendance rates for Lemgo. This aligns quite well with the observed data. Note also that I picked these two teams more or less at random, but estimates are of course available for all teams in the data set.

# Plot observed data using densities and rugs
p1 <- df21 %>%
    filter(home_team %in% c("SG Flensburg-Handewitt", "TBV Lemgo Lippe")) %>%
    ggplot(aes(x = perc, fill = home_team)) +
    geom_density(trim = TRUE) +
    geom_rug() +
    facet_wrap(~ home_team) +
    coord_cartesian(xlim = c(.5, 1), ylim = c(0, 12)) +
    scale_fill_viridis_d() +
    labs(y = "Density", x = "Attendance Rate", fill = NULL,
         title = "Beta Regression: Observed and Fitted Attendance Rates") +
    theme(legend.position = "none")

# Add lines corresponding to estimated distributions
p1 <- p1 + 
    geom_line(data = preds1, aes(x = x, y = y),
              inherit.aes = FALSE, size = 1,
              color = rgb(0, 155, 149, maxColorValue = 255))

# Add predicted means
p1 <- p1 + 
    geom_point(data = data.frame(y = 0,
                                 x = plogis(2.01 + c(1.49, -.35)),
                                 home_team = c("SG Flensburg-Handewitt", "TBV Lemgo Lippe")),
               aes(x = x, y = y), show.legend = FALSE, size = 3,
               color = rgb(0, 155, 149, maxColorValue = 255))
p1

Residuals

Finally, I want to have a quick look at the residuals (which are on the logit scale). I decided to plot them chronologically, that is, against match day. As you see, residuals tend to get higher across the course of a season indicating that attendance rates near the end of the season are slightly higher than predicted. Thus, it seems useful to add variables to the model that can capture this effect (see below).

m1_res <- augment(m1, data = df21)

ggplot(m1_res, aes(x = round_number, y = .resid)) +
        geom_point() +
        stat_smooth(method = "loess", formula = "y ~ x") +
        facet_wrap(~ year) +
        labs(y = "Residuals", x = "Match Day",
             title = "Residuals as a Function of Match Day")

Models Considered

In the following, I will compare different models:

• The models will either model the variable matchday using a linear trend. Or the trend over time will be captured by the categorical variable month (thus allowing different attendance rates in each month). More information about the coding for variable month is provided in part 2.
• Furthermore, I will compare different different link functions or forms:
  • Linear model (lm)
  • Quasi-binomial model (binomial)
  • Beta regression (beta)
  • Beta regression with dispersion parameter and logit link (beta+ logit)
  • Beta regression with dispersion parameter and loglog link (beta+ loglog)

The performance of the models will be evaluated relative to the training data set from above (season 2016/17 and 2017/18) and to a holdout or cross-validation data set (season 2018/19). Furthermore, I will compare the models relative to simply predicting the average attendance rate of the home team.

res11 <- list()

### Average ###

# This "model" simply predicts for each match the yearly average attendance rate
# of the home team
res11$average_no <- 
  glm(perc ~ home_team*year, data = df21,
        family = quasibinomial(logit), weights = capacity)

### beta+dispersion (logit) ###

res11$betaDsplogit_lin <-
    betareg(y ~ home_team + away_team + weekday + year + half2 + 
                matchday | home_team, data = df21)
res11$betaDsplogit_mon <-
    betareg(y ~ home_team + away_team + weekday + year + half2 + 
                Oct + Nov + Dec + Mar + Apr + May | home_team, data = df21)

### beta+dispersion (loglog) ###

res11$betaDsploglog_lin <-
    update(res11$betaDsplogit_lin, link = "loglog")
res11$betaDsploglog_mon <-
    update(res11$betaDsplogit_mon, link = "loglog")

### beta (logit) ###

res11$beta_lin <-
    betareg(y ~ home_team + away_team + weekday + year + half2 + 
                matchday, data = df21)
res11$beta_mon <-
    betareg(y ~ home_team + away_team + weekday + year + half2 + 
                Oct + Nov + Dec + Mar + Apr + May, data = df21)

### lm - linear model ###

res11$lm_lin <-
    lm(attendance ~ home_team + away_team + weekday + year + half2 + 
           matchday, data = df21)
res11$lm_mon <-
    lm(attendance ~ home_team + away_team + weekday + year + half2 + 
           Oct + Nov + Dec + Mar + Apr + May, data = df21)

### quasi-binomial ###

res11$quasibinomial_lin <-
    glm(perc ~ home_team + away_team + weekday + year + half2 + 
            matchday, data = df21,
        family = quasibinomial(logit), weights = capacity)
res11$quasibinomial_mon <-
    glm(perc ~ home_team + away_team + weekday + year + half2 + 
            Oct + Nov + Dec + Mar + Apr + May, data = df21,
        family = quasibinomial(logit), weights = capacity)

res12 <- enframe(res11)

res12
#> # A tibble: 11 x 2
#>    name              value    
#>    <chr>             <list>   
#>  1 average_no        <glm>    
#>  2 betaDsplogit_lin  <betareg>
#>  3 betaDsplogit_mon  <betareg>
#>  4 betaDsploglog_lin <betareg>
#>  5 betaDsploglog_mon <betareg>
#>  6 beta_lin          <betareg>
#>  7 beta_mon          <betareg>
#>  8 lm_lin            <lm>     
#>  9 lm_mon            <lm>     
#> 10 quasibinomial_lin <glm>    
#> 11 quasibinomial_mon <glm>

Model Performance

For all models stored in res12, I will make predictions for the training data df21 and for the holdout data df22. I will then compare these predictions with the actual observations using the following three metrics:

• Mean absolute error (MAE)
• Root mean squared error (RMSE)
• R squared (R2)

R2 should be high, and MAE and RMSE should be relatively low. These latter two measures are provided in the metric of the original attendance values, that is an MAE of 500 means that the predictions are on average off by a value of 500 spectators.

# Holdout data
df22 <- filter(df1, season_id == 55149) %>%
    filter(!home_team %in% c("SG BBM Bietigheim")) %>%
    filter(!away_team %in% c("SG BBM Bietigheim")) %>%
    mutate(year = "17/18") %>%
    filter(date < as.Date("2019-02-01"))

### Make predictions ###

res13 <- res12 %>%
    separate(name, into = c("mod", "matchday")) %>%
    mutate(matchday = factor(matchday)) %>%
    # Make "predictions" for data set used
    mutate(olddata = map(value, ~augment(.x, data = df21, type.predict = "resp"))) %>%
    # Make predictions for data from 2018/19 (cross validation)
    mutate(newdata = map(value, ~augment(.x, newdata = df22, type.predict = "resp"))) %>%
    # Calcuate predictions and residuals for actual attendance number rather
    # than attendance rate
    mutate(olddata = map(olddata, function(x = .x){
        if (x$.fitted[1] > 1) {
            x$.pred <- round(x$.fitted)
        } else {
            x$.pred <- round(x$.fitted*x$capacity)
        }
        x$.res <- x$attendance - x$.pred
        return(x)
    })) %>%
    mutate(newdata = map(newdata, function(x = .x){
        if (x$.fitted[1] > 1) {
            x$.pred <- round(x$.fitted)
        } else {
            x$.pred <- round(x$.fitted*x$capacity)
        }
        x$.res <- x$attendance - x$.pred
        return(x)
    }))

### Performance measures (yardstick package) ###

res14 <- res13 %>% 
    mutate(old = map(olddata, ~yardstick::metrics(.x, attendance, .pred))) %>%
    mutate(new = map(newdata, ~yardstick::metrics(.x, attendance, .pred))) %>%
    # Wrangle into long format
    select(-value, -olddata, -newdata) %>% 
    unnest(c(old, new), names_sep = "_") %>% 
    select(-old_.estimator, -new_.estimator) %>%
    pivot_longer(c(old_.metric, old_.estimate, new_.metric, new_.estimate),
                 names_sep = "_", names_to = c("data", ".value"))

tail(res14)
#> # A tibble: 6 x 5
#>   mod           matchday data  .metric .estimate
#>   <chr>         <fct>    <chr> <chr>       <dbl>
#> 1 quasibinomial mon      old   rmse      652.   
#> 2 quasibinomial mon      new   rmse      995.   
#> 3 quasibinomial mon      old   rsq         0.926
#> 4 quasibinomial mon      new   rsq         0.842
#> 5 quasibinomial mon      old   mae       419.   
#> 6 quasibinomial mon      new   mae       538.

### Polish labels etc ###

res15 <- res14 %>%
    mutate(matchday = fct_recode(matchday, Linear = "lin", Month = "mon"),
           data = fct_rev(fct_recode(data, "Holdout Data" = "new", "Training Data" = "old")),
           .metric = fct_recode(.metric,
                                RMSE = "rmse", R2 = "rsq", MAE = "mae"),
           mod = sub("Dsp(.*)", "+ (\\1)", mod),
           mod = factor(mod, 
                        levels = c("average", "lm", "quasibinomial", "beta", "beta+ (logit)",
                                   "beta+ (loglog)")))

As shown in the plot below, the quasi-binomial model shows the best performance in the training data set. The beta models perform slightly worse but outperforms the quasi-binomial with respect to R2 in the holdout data set. The linear model performs worst (with the exception for RMSE in the holdout data), even worse than predicting just averages without a real model. When comparing the different variants of the beta model, the enhanced models (beta+), which also model the precision in addition to the mean, perform slightly better, especially in the holdout data set. The difference between the logit and the loglog link is tiny and I would prefer the logit model here because it is more common and thus easier to communicate.

The time effect is slightly better captured by the categorical predictor months compared to a linear trend, more on that in part 2.

In summary, both the quasi-binomial and the beta model are able to capture the relevant aspects of the attendance rates of matches of the German Handball-Bundesliga. The beta model, however, has the advantage that it can provided prediction intervals if desired whereas the intervals of the quasi-binomial model are way too narrow with data of several thousand spectators (not shown herein).

ggplot(res15, aes(y = .estimate, x = mod, color = matchday)) +
    geom_point(position = position_dodge(width = .5), size = 3) +
    facet_wrap(data ~.metric, scales = "free_y") +
    scale_color_viridis_d(begin = .2) +
    labs(x = "Model", y = "", color = "Time\nTrend",
         title = "Comparison of Different Models for Proportion Data") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))