Generalized Mixed Effects Models?

Introduction

In the NFL, teams can score points in three ways:

1. Safety: Two points

• Occurs less than 1% of possessions

2. Field Goal: Three points

• Occurs about 40% of teams’ offensive possessions

3. Touchdown: Six points

• Occurs about 20% of teams’ offensive possessions

A field goal attempt occurs when a team attempts to kick the ball through goal placed at both ends of the playing field. If the ball passes between the uprights, the team scores 3 points. Otherwise, the other team’s offense starts where the kicking team was on the field. If a team misses a field goal, not only do they lose out on three points, but the opponent has a higher probability of scoring points!

There are two types of stadium roofs that teams play in: open or closed. One common comment about stadiums with closed roofs (called domes) is teams have an easier time scoring since the weather doesn’t have any influence on the result of a play. We want to examine the data to determine if there is any legitimacy of these types of comments for field goal attempts.

Getting the data

The code chunk below will load all NFL plays from the last 5 years and clean it to keep only field goal attempts made by players that attempted at least 25 field goals.

We’ll start just by loading in the NFL plays data for the last 6 seasons and keeping just the field goal attempts

# Loading the most recent 6 seasons
pbp <- load_pbp((year(Sys.Date()) - 5):year(Sys.Date()))

FieldGoals <- 
  pbp  |>  
  
  # Only keeping attempts without penalty during the regular season
  filter(
    field_goal_attempt == 1, 
    penalty==0, 
    season_type=="REG"
  ) |> 
  
  mutate(
    # Changing roof_type to just indicate if the kick took place in or outdoors
    roof_type = if_else(condition = roof %in% c("outdoors", "open"), 
                        true = "outdoors",
                        false = "indoors"),
    
    # If the field goal was successful or not   
    fg_result = if_else(condition = field_goal_result == "made", 
                        true = 1, 
                        false = 0)
  ) |> 
  
  # Changing the name to something simpler  
  rename(
    kicker = kicker_player_name, 
    distance = kick_distance, 
    kicker_id = kicker_player_id
  ) |> 
  
  # Picking the variables we need  
  dplyr::select(
    team = posteam, kicker, kicker_id, distance, 
    roof_type, fg_result, season
  )

FieldGoals

## # A tibble: 5,862 × 7
##    team  kicker     kicker_id  distance roof_type fg_result season
##    <chr> <chr>      <chr>         <dbl> <chr>         <dbl>  <int>
##  1 SF    R.Gould    00-0023252       52 outdoors          1   2020
##  2 ARI   Z.Gonzalez 00-0033862       52 outdoors          0   2020
##  3 SF    R.Gould    00-0023252       24 outdoors          1   2020
##  4 ARI   Z.Gonzalez 00-0033862       56 outdoors          1   2020
##  5 ARI   Z.Gonzalez 00-0033862       49 outdoors          0   2020
##  6 DET   M.Prater   00-0023853       27 indoors           1   2020
##  7 CHI   C.Santos   00-0031203       35 indoors           1   2020
##  8 DET   M.Prater   00-0023853       32 indoors           1   2020
##  9 CHI   C.Santos   00-0031203       28 indoors           1   2020
## 10 DET   M.Prater   00-0023853       44 indoors           1   2020
## # ℹ 5,852 more rows

# Removing the pbp data set since it is large and not used again
rm(pbp)

The seven variables are:

1. team: The team attempting the field goal

2. kicker: The player who attempted the field goa

3. distance: How far away the uprights are from the player attempting the field goal

4. roof_type: One of two different values:

• outdoors: a stadium that doesn’t have a roof or a stadium with the roof open
• indoors: a stadium with a closed roof

5. fg_result: If the field goal attempt was successful

• fg_result = 1 -> the attempt was made
• fg_result = 0 -> the attempt was missed

6. season: The season the FGA took place during

The code below will count how many field goals each player attempted and only keep the field goals attempted by kickers with at least 25 attempts:

kicker_counts <- 
  FieldGoals |> 
  # Counting the number attempts for each player
  count(kicker) |> 
  # Only keeping players with at least 25 kicks
  filter(n >= 25) 



kicker_counts |> 
  arrange(-n)

## # A tibble: 49 × 2
##    kicker          n
##    <chr>       <int>
##  1 D.Carlson     204
##  2 C.Boswell     185
##  3 Y.Koo         181
##  4 J.Myers       180
##  5 M.Gay         177
##  6 B.McManus     176
##  7 J.Tucker      176
##  8 K.Fairbairn   175
##  9 N.Folk        175
## 10 J.Slye        172
## # ℹ 39 more rows

There are a total of 51 kickers that attempted at least 25 field goals during the most recent six seasons. We’ll keep only the field goals attempted by these kickers to remove any players that were not accurate enough to be an NFL quality field goal kicker.

fg_df <- 
  FieldGoals |> 
  filter(kicker %in% kicker_counts$kicker) 

tibble(fg_df)

## # A tibble: 5,535 × 7
##    team  kicker     kicker_id  distance roof_type fg_result season
##    <chr> <chr>      <chr>         <dbl> <chr>         <dbl>  <int>
##  1 SF    R.Gould    00-0023252       52 outdoors          1   2020
##  2 ARI   Z.Gonzalez 00-0033862       52 outdoors          0   2020
##  3 SF    R.Gould    00-0023252       24 outdoors          1   2020
##  4 ARI   Z.Gonzalez 00-0033862       56 outdoors          1   2020
##  5 ARI   Z.Gonzalez 00-0033862       49 outdoors          0   2020
##  6 DET   M.Prater   00-0023853       27 indoors           1   2020
##  7 CHI   C.Santos   00-0031203       35 indoors           1   2020
##  8 DET   M.Prater   00-0023853       32 indoors           1   2020
##  9 CHI   C.Santos   00-0031203       28 indoors           1   2020
## 10 DET   M.Prater   00-0023853       44 indoors           1   2020
## # ℹ 5,525 more rows

Our cleaned data has a total of 5,535 kick attempts! They will be used to try to answer if it actually is easier to make a field goal indoors.

Logistic Regression

We can estimate the probability a field goal is successful based on distance and location (indoor/outdoor) using logistic regression.

The logistic regression model is:

\[\log\left(\frac{\pi(x)}{1-\pi(x)}\right) = \beta_0 + \beta_1x_1 + \beta^I_2 + \beta_{12}^Ix_1\]

where

• \(\beta_0\) is the log odds of a successful attempt if it is 0 yards out and in an outdoor stadium

• \(\beta_1\) is the effect of distance on the log odds of success for each 1 additional yard out

• \(\beta^I_2\) is the difference in the log odds of a successful attempt being indoors vs outdoors

• \(\beta_{12}^I\) is the interaction term between distance and roof type.

Summarized Data

# Summarizing the data first:
fg_sum <- 
  fg_df |> 
  # Count for each combo of distance, roof_type, and result
  count(distance, roof_type, fg_result) |> 
  
  # Calculating the total number of attempts per distance
  mutate(
    .by = distance,
    dist_n = sum(n)
  ) |> 
  
  # Calculating the success rate by distance and roof type
  mutate(
    .by = c(distance, roof_type),
    dist_field_n = sum(n),
    success_rate = n/dist_field_n
  ) |> 
  
  # Only keeping rows for the success rate and distance above 20
  filter(
    fg_result == 1,
    dist_n > 20
  ) |> 
  
  dplyr::select(-fg_result, -dist_n, -n) |> 
  
  rename(attempts = dist_field_n)

fg_sum

## # A tibble: 80 × 4
##    distance roof_type attempts success_rate
##       <dbl> <chr>        <int>        <dbl>
##  1       20 indoors         20        1    
##  2       20 outdoors        41        0.976
##  3       21 indoors         23        1    
##  4       21 outdoors        63        1    
##  5       22 indoors         35        1    
##  6       22 outdoors        73        0.973
##  7       23 indoors         50        1    
##  8       23 outdoors        85        1    
##  9       24 indoors         40        0.975
## 10       24 outdoors        79        1    
## # ℹ 70 more rows

Fixed Effect Model

# Fitting a fixed effect logistic regression model with an interaction term
fg_fixed <- 
  glm(formula = success_rate ~ distance * roof_type,
      # Swapping the order of indoor and outdoor so outdoor is the baseline
      data = fg_sum |> mutate(roof_type = fct_rev(roof_type)),
      weights = attempts,
      family = binomial)


# Checking the model estimates using the tidy function and rounding them to 3 digits
library(broom)
tidy(fg_fixed) |> 
  mutate(across(.cols = where(is.numeric),
                .fns = round,
                digits = 3)) |> 
  gt::gt()

term	estimate	std.error	statistic	p.value
(Intercept)	6.205	0.274	22.614	0.000
distance	-0.105	0.006	-17.663	0.000
roof_typeindoors	-0.094	0.490	-0.192	0.848
distance:roof_typeindoors	0.008	0.010	0.743	0.457

When examining the necessity of each model term, we should always start by examining the interaction term (we should always maintain a hierarchical model). From the table above, the p-value for the interaction term (distance:roof_typeindoors) is 0.457, indicating the interaction term is not statistically significant and the effect of distance on the probability a field goal attempt is successful the same for an indoor or outdoor stadium.

Let’s reparameterize and refit the model without the interaction term:

\[\log\left(\frac{\pi(x)}{1-\pi(x)}\right) = \beta_0 + \beta_1x_1 + \beta^I_2\]

# Fitting a fixed effect logistic regression model without an interaction term
fg_fixed2 <- 
  glm(formula = success_rate ~ distance + roof_type,
      data = fg_sum |> mutate(roof_type = fct_rev(roof_type)),
      family = binomial,
      weights = attempts)


# Checking the model estimates using the tidy function and rounding them to 3 digits
glm_estimates <- 
  tidy(fg_fixed2, 
       #exponentiate = T, # convert from log odds to odds
       conf.int = T) |>  # calculate a 95% confidence interval
  # Rounding the numeric columns to 3 decimal places
  mutate(
    across(
      .cols = where(is.numeric),
      .fns = round,
      digits = 3
    )
  ) 


glm_estimates |> 
  gt::gt()

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	6.093	0.227	26.859	0.000	5.656	6.546
distance	-0.102	0.005	-20.975	0.000	-0.112	-0.093
roof_typeindoors	0.264	0.088	3.001	0.003	0.093	0.438

Now both the effect of distance and roof_type are significant!

The odds a kick is successful when attempted indoors is about 30% higher if it is attempted at the same distance.

Interpreting the confidence interval effect, we are 95% confident that the odds an indoor kick is successful is 10% to 55% higher compared to an outdoor kick when attempted at the same distance.

The conclusions we can arrive from our model are only accurate if the model fits the data reasonably well. Let’s look at the fit statistics:

fit_stats_logit_fixed <- 
  glance(fg_fixed2) |> 
  mutate(p.val = pchisq(q = deviance, df = df.residual, lower.tail = F)) |> 
  select(nobs, deviance, df.residual, p.val)

gt::gt(fit_stats_logit_fixed)

nobs	deviance	df.residual	p.val
80	94.8	77	0.08231

The p-value is for a hypothesis test with the null hypothesis stating the model fits the data and the alternative hypothesis is that the model does not fit the data. With a p-value of 0.0823, there is some evidence that the model doesn’t fit the data, but not strong evidence. We’ll continue on the assumption that the model is a (some-what) reliable generalization about the relationship between field goal success rate, distance, and roof type.

Let’s graph the results of the logistic regression:

# Adding the estimated probability of success at 59 yards away
pred_df <- 
  fg_sum |> 
  mutate(
    est_success_rate = predict.glm(object = fg_fixed2, 
                                   newdata = fg_sum, 
                                   type = "response")
  )

# Creating the graph
ggplot(
  data = pred_df,
  mapping = aes(
    x = distance,
    
    color = roof_type
  )
) +
  
  geom_point(
    mapping = aes(y = success_rate)
  ) + 
  
  # Adding the line of estimated probabilities
  geom_line(
    mapping = aes(y = est_success_rate)
  ) +
  
  # Adding indoor and outdoor and the end of their respective line
  geom_text(
    data = pred_df |> filter(distance == max(distance)),
    mapping = aes(
      y = est_success_rate,
      label = roof_type
    ),
    show.legend = F,
    nudge_x = 2.3
  ) +
  
  theme_bw() + 
  
  labs(
    title = "Observed and Estimated Field Goal Success",
    x = "Kicking Distance",
    y = "Success Rate",
    color = NULL
  ) + 
  
  scale_y_continuous(labels = scales::percent) + 
  
  scale_color_manual(values = c("#624a2e", "#567d46"),
                     labels = c("Indoors", "Outdoors")) +
  
  scale_x_continuous(breaks = 2:5*10,
                     labels = paste(2:5*10, "yds")) +
  
  theme(
    plot.title = element_text(hjust = 0.5),
    legend.position = c(0.1,0.2)
  )

Logistic regression model: \(\log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1x_1 + \beta^I_2\)

Mixed Effects Model

One big assumption most models have is that the rows in our data are independent once we account for the explanatory variables (distance and location in our data).

However, there is still one source of association between the rows (each field goal attempt): the players attempting the field goals!

We can account for the kickers attempting the field goals by including a random effect in the model:

\[\log\left(\frac{\pi(x)}{1-\pi(x)}\right) = \beta_0 + \beta_1x_1 + \beta^I_2 + \beta_{12}^Ix_1 + u_i\]

where \(u_i\) is the random effect of a specific kicker (kicker \(i\)) on the log odds a field goal attempt is successful.

The difference between a fixed effect (like \(\beta_1\)) and a random effect (like \(u_i\)) is that we assume a probability distribution for the random effect. Typically, we assume the random effects are Normally distributed:

\[u_i \sim N(0, \sigma_u)\]

The code chunk below will estimate the mixed-effects logistic regression model. It’s called a mixed-effects model because it has both fixed effects and random effect(s).

# Fitting the mixed effects model using glmer in the lme4 package
library(lme4)
fg_mem <- 
  glmer(
    fg_result ~ distance * roof_type + (1|kicker), # (1|kicker) -> random effect
    data = fg_df |> mutate(roof_type = fct_rev(roof_type)), 
    nAGQ = 50,  # How many quadrature points to use when estimating u_i
    family = binomial, 
    control = glmerControl(optimizer = "bobyqa")
  )

# Viewing the regression estimate table
summary(fg_mem)$coefficients |> 
  data.frame() |> 
  round(digits = 3) |> 
  rownames_to_column(var = "Term") |> 
  mutate(Estimate = exp(Estimate)) |> 
  rename(
    SE = Std..Error,
    `p-value` = Pr...z..,
    `test stat` = z.value
  ) |> 
  gt::gt()

Term	Estimate	SE	test stat	p-value
(Intercept)	585.2271	0.271	23.471	0.000
distance	0.8967	0.006	-18.927	0.000
roof_typeindoors	0.8098	0.476	-0.443	0.657
distance:roof_typeindoors	1.0111	0.010	1.085	0.278

Like with the fixed effects model, the interaction term doesn’t appear to be needed since the p-value is large. It will be removed and the model reran:

fg_mem2 <- 
  glmer(
    fg_result ~ distance + roof_type + (1|kicker), 
    data = fg_df |> mutate(roof_type = fct_rev(roof_type)), 
    nAGQ = 50,
    family = binomial, 
    control = glmerControl(optimizer = "bobyqa")
  )

summary(fg_mem2)$coefficients |> 
  data.frame() |> 
  round(digits = 3) |> 
  rownames_to_column(var = "Term") |> 
  mutate(Estimate = exp(Estimate)) |> 
  rename(
    SE = Std..Error,
    `p-value` = Pr...z..,
    `test stat` = z.value
  ) |> 
  gt::gt()

Term	Estimate	SE	test stat	p-value
(Intercept)	499.6960	0.225	27.584	0.000
distance	0.8994	0.005	-22.248	0.000
roof_typeindoors	1.3458	0.091	3.273	0.001

Like with our fixed effects model, we have strong evidence that the location has an effect on the probability a field goal attempt is successful after accounting for the effect of distance and the person attempting the field goal.

If a kick is attempted at the same distance by the same player, the odds a kick is successful is about 35.7% higher if it is attempted in an indoor stadium compared to an outdoor stadium.

Let’s visualize the difference using the mixed effect model. We’ll start by creating a data frame that has each combo of distance, stadium type, kicker, and the success probability estimated by the mixed effects model.

# Creating a data frame that has a row for each combination of distance, roof_type, and kicker to use to estimate the probability a kick is successful:
pred_data <- 
  expand.grid(
    distance = seq(20, 55, by = 0.250),   # Distance by quarter yard increments
    roof_type = unique(fg_df$roof_type),  # Two types of roof
    kicker = unique(fg_df$kicker)         # kickers included in the model
  )        


# We can use the predict function to calculate the estimated probabilities 
fg_predict <- 
  data.frame(
    pred_data, 
    est_made = predict(fg_mem, 
                       newdata = pred_data, 
                       type = "response")
  )

The graph below shows the median, middle 50%, and middle 90% of model probabilities that a field goal attempt is successful.

# Calculating the 5 number summary of success for each distance and roof_type combo
avg_probs_dist <- 
  fg_predict |>
  summarize(
    .by = c(roof_type, distance),
    avg_prob_made = mean(est_made),
    min_prob = min(est_made),
    Q1_prob = quantile(est_made, prob = 0.25),
    Q3_prob = quantile(est_made, prob = 0.75),
    max_prob = max(est_made),
    quantile95 = quantile(est_made, prob = 0.95),
    quantile05 = quantile(est_made, prob = 0.05)
  )

# Now we'll plot the results:
ggplot(
  data = avg_probs_dist,
  mapping = aes(
    x = distance, 
    y = avg_prob_made, 
    fill = roof_type
  )
) + 
  
  # Adding the line for the mean estimated probability
  geom_line(
    mapping = aes(color = roof_type),
    linewidth = 1
  ) + 
  
  # Adding a confidence region for Q1 and Q3 
  geom_ribbon(
    mapping = aes(
      ymin = Q1_prob, 
      ymax = Q3_prob
    ),   
    alpha = 0.5, 
    linetype = "dashed"
  ) +
  
  annotate(
    geom = "text",
    label = "Shaded regions indicate the estimated\nmiddle 50% and 90% ranges",
    x = 26,
    y = 0.8
  ) +
  
  geom_ribbon(
    mapping = aes(
      ymin = quantile95,
      ymax = quantile05
    ),
    alpha = 0.25,
    linetype = "dotted"
  ) +
  
  #facet_wrap(~roof_type) + 
  
  # Changing the fill and color pallete to match the previous color choices
  scale_color_manual(
    values = c("#624a2e", "#567d46"), 
    aesthetics = c("colour","fill")
  ) + 
  
  # Changing the y-axis to percentages
  scale_y_continuous(labels = scales::percent) +
  
  # Changing the x-axis to include yards
  scale_x_continuous(
    breaks = 2:6*10,
    labels = paste(2:6*10, "yds")
  ) + 
  
  # Adding context
  labs(
    title = "Probability of a Successful Field Goal by Distance for 
             <span style='color:#624a2e;'>Indoor</span> and 
             <span style='color:#567d46;'>Outdoor</span> Stadiums",
    subtitle = "Seasons 2019-2024, Min 25 Attempts per Kicker, Between 19-58 Yards",
    x = "Field Goal Distance",
    y = "Model Probability of Success",
    caption = "data: nflfastR"
  ) + 
  
  # Moving the legend and centering the subtitle
  theme(
    legend.position = 'none',
    plot.title = ggtext::element_markdown(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

Mixed effect regression: \(\log\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1x_1 + \beta^I_2 + u_i\)

Indoor kicks are more likely to be made, and the further the distance, the larger the difference in the success rate! The difference is most notable for attempts 45 yards and further.

Comparing Kickers

While it doesn’t help answer our original question “Is it easier to attempt a field goal indoors”, we can use the random effects generated from our mixed effects model to answer “Who is the best kicker over the last 6 seasons, accounting for distance and stadium type?”

If the random effect is positive, the kicker is better than league average and if the random effect is negative, the kicker is worse than league average. The larger in absolute value, the better (or worse) the player is.

# Let's look at the effectiveness of each player
kicker_eff <- 
  ranef(fg_mem) |> 
  data.frame() |> 
  select(kicker = grp, condval)

kicker_eff <- 
  left_join(
    x = kicker_eff,
    y = kicker_counts,
    by = "kicker"
  )

# The top and bottom 5 kickers
bind_rows(
  .id = "type",
  "best"  = kicker_eff |> slice_max(condval, n = 5),
  "worst" = kicker_eff |> filter(n > 50) |> slice_min(condval, n = 5)
) |> 
  arrange(-condval) |> 
  mutate(`odds of success` = exp(condval) |> round(2)) |> 
  select(kicker, attempts = n, `odds of success`) |> 
  gt::gt()

kicker	attempts	odds of success
N.Folk	175	1.45
B.Aubrey	102	1.35
C.Boswell	185	1.32
C.Dicker	123	1.32
E.Pineiro	119	1.28
B.Grupe	91	0.86
C.Ryland	75	0.86
G.Joseph	121	0.84
M.Badgley	95	0.80
J.Moody	70	0.78

According to the model, Chris Boswell has been the best kicker in the league over the last six years, with the odds of a successful field goal attempt 33% above the league average. Brandon Aubrey, Justin Tucker, and Nick Folk are tied for second place, all having 27% odds of a successful field goal attempt above league average.

The worst kicker in the league (with a minimum of 50 attempts) over the last 6 seasons is Ryan Succop. The odds his field goal attempt is successful is 20% below the league average. The worst 5 players are as bad as the best five players are good. This makes sense since NFL kickers aren’t given many chances, and if they start performing below the league average, they won’t have a job for very long!

Let’s visualize kickers’ random effects for players who have at least 100 attempts:

# Loading needed packages
pacman::p_load(nflplotR)


gg_kicker_effects <- 
  kicker_eff |> 
  # Only keeping rows with at least 100 attempts
  filter(n >= 100) |> 
  
  left_join(y = FieldGoals |> select(kicker, kicker_id) |> distinct(),
            by = "kicker") |> 
  
  # blue = above average, red = below average and reordering the kickers
  mutate(
    bar_col = if_else(condval > 0, "#013369", "#D50A0A"),
    kicker = fct_reorder(kicker, condval)
  ) |> 
  
  # Bar chart for kicker effects
  ggplot(
    mapping = aes(
      x = fct_rev(kicker),
      y = abs(condval),
      fill = bar_col
    )
  ) + 
  
  # Adding a line for the bars to sit on at 0
  geom_hline(
    linewidth = 1,
    yintercept = 0
  ) +
  
  geom_col(
    color = "black"
  ) +
  
  annotate(
    geom = "text",
    color = "#013369",
    label = "More likely to make the attempt",
    x = "C.Santos",
    y = .3,
    fontface = "bold"
  ) + 
  
  annotate(
    geom = "text",
    color = "#D50A0A",
    label = "More likely to miss the attempt",
    x = "W.Lutz",
    y = 0.3,
    fontface = "bold"
  ) +
  
  scale_fill_identity() + 
  
  scale_y_continuous(
    limits = c(0, 0.35),
    breaks = (0):5/10,
    expand = c(0,0),
    labels = round(exp(0:5/10), 1)
  ) + 
  
  labs(
    x = NULL,
    y = "Odd of success", 
    title = "Odds of success for NFL Kickers",
    subtitle = "2019 - 2024 regular season; Minumum 100 attempts"
  ) + 
  
  theme(
    axis.text.x = element_text(
      angle = 90,
      hjust = 0.95,
      vjust = 0.15
    )
  ) +
  geom_nfl_headshots(
    mapping = aes(player_gsis = kicker_id),
    width = 0.08,
    height = 0.08,
    vjust = 0
  )  



gg_kicker_effects

Odds of success for a kicker accounting for distance and location

There are more players with a positive effect on the odds a field goal is successful. A possible explanation is if a kicker has a below average chance of success, they are cut from the team and are unlikely to get more than 100 attempts and wouldn’t be included in the graph.