Shrinkage estimation of NFL field goal success probabilities

2.1Data collection: web scraping with R

Like most other work in the area, we assume kicks to be independent and model the dependence of success probability on explanatory variables using generalized linear models. The explanatory information we consider includes only the distance of the attempt and who the kicker is. Using the readLines() command and regular expressions functionality in the R statistical software package (R Core Team, 2014),, the distance of a field goal attempt is gleaned from the website with url. This website has complete play-by-play information going back as far as 1998. More specifically, the string “yard field goal" is matched and the kicker and distance occurring in the string before the match are recorded along with the outcome that follows the match. Every field goal attempted, even those on which a penalty occurred, but for which a result was given in the play-by-play transcript, are included in the dataset. The task loops over seasons 1998-2014, weeks within season (including the NFL playoffs), games within week, and finally plays within game.

The dataset contains 17, 104 attempts from more than K = 111 kickers. Kickers who made or missed all of their career attempts are excluded from consideration to avoid separability issues with maximum likelihood estimation. (Jaret Holmes played for three teams from 1999-2001, making all five of his field goal attempts and all four of his extra point attempts successfully. Danny Boyd was even better, making 5/5 field goals and 7/7 extra points for Jacksonville in 2002.)

Outcomes were dichotomized as “good” or “not good”, the latter event includes blocks and fumbled holds. Our reasoning is that the blame for failed attempts usually rightly falls on the kicker, but sometimes the snap, the hold or the blocking are at fault. The records online do not easily allow for determination of blame. An example of an entry that was scraped from the url (http://www.pro-football-reference.com/boxscores/199809060cin.htm). (a game from week 1 of the 1998-1999 season) appears in Fig. 1.

Fig.1

Play-by-play item containing relevant information on a field goal attempt.

Occasionally, comparisons between total kick frequencies computed from play-by-play data with kicker career totals will reveal slight discrepancies, where certain plays are absent from the pro-football-reference play-by-play information. For example, in a 1998 Week 1 game of the New England Patriots at the Denver Broncos, Adam Vinatieri of the New England Patriots had a 37 yard field goal attempt during the second quarter. The pro-football-reference play-by-play account of the event is blank (see Fig. 2). We resorted to the site to attempt to rectify this particular record. Though we could not find links on the site to play-by-play information for games before 2001 (a situation in which automated screen-scraping is then difficult), we were able to find some urls, with links provided by a search engine.

Fig.2

Play-by-play item for a blocked kick by Adam Vinatieri in a 1998 Week 1 game that appeared as a blank on pro-football-reference.com

2.2Choice of a link function

For notation, we let k = 1, …, K index the kickers (K = 111 in our application) and let j = 1, …, N (N = 17104) index all of the field goal attempts. For each attempt, we let binary outcome O_j = 1 if the kick is good, and O_j = 0 otherwise. The probability of success of attempt j is denoted π_j. The general form for the log-likelihood function of all of the probability parameters {π_j} for independent, dichotomous data with outcomes o_j is given by

Generalized linear models adopt a link function g, through which the probabilities π_j are transformed to an expression that is linear in regression coefficients introduced to model the effects of available explanatory variables (and reduce the dimension of the parameter space). We define d_j as the distance of the jth attempt and define indicator variables for kicker k as X_jk = I (jthattemptisbykicker k), for k = 1, …, K, with regression coefficients β₁, …, β_K. The probability parameters are now, after transformation by the link function g, linear in the regression coefficients,

We refer to β_k as a kicker-specific slope and can test the significance of estimated regression coefficients by comparison with a generalized linear model with only one slope, β, to accommodate all kickers,

In applications in which statistical regression models are developed for binary response variables, logistic regression seems to be the most popular choice (Pasteur and Cunningham-Rhoads, 2014; Morris and Rolph, 1981; Bilder and Loughin, 1998; Berry and Wood, 2004), perhaps because of its interpretability, particularly through odds ratios. However, several other link functions are easily explored using statistical software for generalized linear models. We introduce a superscript notation for probability functions associated with kickers that will simplify the shrinkage estimation exposition in Section 3. Let π^k (d) denote the success probability function for kicker k at any distance d. We consider three choices for the link function: logistic, probit, and CLL. Let Φ and Φ^-1 denote the cumulative distribution function and quantile function, respectively, of a random variable with standard normal distribution. Then the links and inverse links (probability functions) for these three candidates are given below:

The log-likelihood function of a specified link function uses the reparameterizations above. As an example the log-likelihood for the logistic link is given by

Table 2 summarizes the maximized log-likelihood function when using the logit, probit, and CLL link functions for estimating field goal success probability in our data set. Tests comparing the nested models within a column of the table are all highly significant on 110 degrees of freedom (comparing models with K = 111 different slopes to a model with a single slope), indicating significant differences among the estimated slopes. With such a large sample, there does appear to be enough data to detect potentially subtle differences among kickers.

Figure 3 presents the inverse link functions applied to the empirical success frequencies against distance. The lines overlaid on the plots are from maximum likelihood estimation of the single-slope model. Residuals computed as standardized differences between transformed success frequencies and the estimated link functions are plotted against distances in Fig. 4, providing an enhanced visualization of the fits. Link-specific lowess smoothers overlaying these plots suggest that the residuals from the CLL fit exhibit the least dependence on distance.

Fig.3

The fits provided by three inverse link functions. Plot (a) presents the three fits to the empirical frequency data: logit link fit solid/blue, probit link dashed/red, complementary log-log link dot-dashed/green. Plots (b), (c) and (d) present the fits after link transformation. All four plots have distance in yards as the horizontal axis.

Fig.4

Residuals from three link functions in generalized linear models (logit as blue squares, probit as red triangles and complementary log-log as green circles). Residuals are differences from link-transformed empirical success frequencies and estimated link functions and have been normalized to have unit variance. Plotted against distances with at least 100 attempts, with Lowess smoothers overlaid using solid lines.

In an analytical comparison of link functions, McCullagh et al. (1989) observes that for trials with intermediate success probability, 0.1 ≤ π ≤ .9 (the vast majority of field goal attempts), probit and logistic functions are nearly linearly related, so that discrimination between the fits of models using these two links is difficult. Additionally, for trials with success probability approaching one (which occur as field goal attempt distances decrease towards 17 yards), the CLL approaches one more slowly than either of the other two link functions. The CLL link can also differ (Dobson and Barnett, 2011) from logit and probit links for small values of π, which are often of greatest interest when assessing kickers. Close inspection of Fig. 3 and Fig. 4 seem to reflect these observations, with better agreement with empirical frequencies on short kicks than that provided by either logit or probit. Fitting these models using kicks at the extremes may be problematic in that misses at short distances are likely due to mishandled snaps and attempts at long distances may lead to selection bias since these attempts are probably only afforded to very good kickers. In light of these issues, log-likelihood ratios are also reported after restricting attention to kicks with distances between 25 and 50 yards. With this restriction, the CLL link achieves the maximum likelihood among the link functions for models with either a single slope or with kicker-specific slopes.

The statistic proposed by Hosmer and Lemeshow (2004) (HL) offers another way to assess goodness-of-fit in regression models for binary data with continuous explanatory variables. The HL statistic is a Pearson-style measure of goodness-of-fit whose sampling distribution under a correctly specified model can be approximated by a chi-square distribution. The HL procedure partitions all kicks into 10 groups of roughly equal size. (We accomplished this using the algorithm laid out in the documentation for PROC LOGISTIC within the SAS statistical software package (SAS Institute, 2011). The procedure then computes, for each group, observed (O_i) and expected (E_i) counts of successful kicks by averaging estimated probabilities within group, in order to compute the χ² statistic, X:

The summary statistics necessary for computation of the HL statistic are given in Table 3 for the logit and CLL links. The logistic regression model does not provide a good fit for the 10th group, where the estimated success probability is highest. The average estimated probability among these 1708 kicks is πˆ¯10=0.9692 for an estimated expected cell count of e₁₀ = 1708 × 0.9692 = 1655.34. However, the observed count is o₁₀ = 1676 successful kicks, resulting in a large contribution, (1676 - 1655.3) ²/{1655.3 (1 -0.9692)} =8.4, to the chi-square statistic.

Table 4 gives HL chi-square statistics andp-values on 8 degrees of freedom for each of the three link functions with either a single slope, or with kicker-specific slopes. The models with the CLL link function are the only ones that do not exhibit significant lack-of-fit. With such a large sample size (N = 17, 104 attempts), even subtle lack-of-fit is likely to be detected, as in the logistic (p = .0291,. 0657) and probit (p < .0001, p = .0001) models. Despite this significance, neither of the two provides a bad fit.

Lastly, a link function could be chosen on the basis of how well probabilities estimated using the link can predict successful and unsuccessful field goal attempts. A decision rule which predicts a successful field goal when the estimated success probability exceeds some threshold (πˆ>π0) has as its sensitivity and specificity the corresponding relative frequencies of successes among kicks which fall, respectively, above or below this threshold. A plot of sensitivity against one minus specificity is called the receiver operating characteristics (ROC) curve and the area under this curve (AUC) is a common way of assessing the capacity of a fitted model to predict successes and failures. The AUC for models with a single slope that are not kicker-specific will necessarily be the same for all three links, as the ordering of the estimated probabilities is determined only by the distance of the attempts. However, small differences in AUC exist among the three links when kicker-specific slopes are included in the model. This ordering is consistent with that implied by visual inspection of Fig. 3. The AUC values, all close to 0.75, are presented in Table 5.

Since the CLL link enjoys the smallest HL and largest AUC under cross-validation (as described in Section 3.1), computations under the other two link functions are excluded from presentation in the remainder of the paper.

3.1Proposed shrinkage estimators

For shrinkage estimators of the success probability at a given distance, π (d), we consider estimators that are weighted averages of the MLE and empirical frequency estimators defined above, πˆk(d) and πˆEF(d), respectively. These shrinkage estimators may be written in a general form using a weighting function ω:

The first shrinkage estimator we propose is simply the midpoint between the two components. We refer to this shrinkage estimate, obtained by setting ω = 1/2, as the midpoint estimator (denoted by superscript M):

One approach towards constructing a shrinkage estimator that weights model-based kicker-specific estimates more heavily is simply to weight components in inverse proportion to their standard errors. This will upweight the MLE, πˆk(d) for veteran kickers for whom more information is available. For the MLE, we rely upon the GLIMMIX procedure of SAS to produce asymptotic standard errors to accompany the maximum likelihood estimator, denoted by SEk=SE(πˆk(d)). For the empirical frequency component πˆEF(d), we use SEEF=SE(πˆEF(d))=πˆEF(d)(1-πˆEF(d))/Nd where N_d denotes the total number of kicks attempted at distance d in the entire dataset. If the empirical frequency is 0, such as at d = 65, we used SEEF=0.5·0.5/Nd. The resulting inverse variance (IV) weighted shrinkage estimate can then be written in the general weighted average form by taking ω = SE_EF/(SE_k + SE_EF):

Note that because sample sizes differ, this weighting scheme can lead to different weights at different distances.

Another approach is to choose a weighting function for the MLE that is increasing in n_k towards a maximum weight of one. Any cumulative distribution function (CDF) has this property. We choose the exponential distribution function, which for a suitably chosen rate parameter a, increases rapidly for small n with diminishing returns on large n_k:

This choice leads to the exponentially weighted shrinkage estimator,

The quantity a in the exponential weights estimator can be viewed as a tuning parameter, whose value is “fit” using the entire dataset and whose performance is evaluated by cross-validation (see Section 3.2). For the a parameter, mean values 100 < a < 700 were considered, reasoning that the kick counts fell in the interval 4 < n_k < 572. A search of values of a, to the nearest multiple of 10, achieves a Hosmer-Lemeshow (HL) statistic of HL = 9.32 at an optimized value of a = 670. Using the exponential CDF with mean a = 670 for a weight function, the corresponding weights on the model-based estimates for kickers with sample sizes of n_k = 10,  50,  100,  200,  400 would be 0.02,  0.07,  0.14,  0.26,  0.45, respectively.

Figure 5 presents a plot of the weight function, ω for these different shrinkage estimators. The two separate curves provided for the IV weighted estimator corresponds to kicks at d = 45 and d = 55 yards. The empirical frequencies at these two distances were πˆEF(45)=353/484=0.73 and πˆEF=26/52=0.50, respectively. For these two distances, plots against sample size suggest that the model-based standard errors may be roughly approximated by 0.64/n and 0.54/n, respectively. Because of the dependence of the empirical frequency on sample size, with little information to go by for kicks as long as d = 55 yards, the IV weighted estimator places more weight on the first component for a kick at d = 55 yards than one at d = 45 yards. All of these curves, with vertical axis given on the right of the plot, are superimposed over a frequency histogram showing how the N = 17, 104 attempts are distributed across the K = 111 kickers. The frequencies of kickers with different experience levels is given on the left vertical axis.

Fig.5

Weight functions for shrinkage estimation, “exp” plots the weights based on the exponential distribution and IV plots the weights proportional to the inverse of the variances of the components.

The estimators of success probability can be naturally shrunk towards central values by means of a Bayesian generalized linear model with CLL link utilizing empirical Bayes estimation (see e.g.,Carlin and Louis (1997)). In such a model, as fit using the bayesglm function of the arm package in R (Gelman et al., 2014), the regression coefficients are assumed to be independent, normally distributed random variables with mean and variance characterized by hyper-parameters that are specified from data, in our case:

The hyperparameters for the intercept are the default choices for the bayesglm function and those for the kicker-specific slope were chosen from the sample mean and variance of maximum likelihood estimates from the frequentist generalized linear model fit.

Fig. 7 illustrates the reduction in variance due to shrinkage. For attempt distances, d = 30,  40,  50,  60 yards, random samples of 10 kicks were selected. The estimated success probability estimates are plotted using different colors, with estimation technique on the horizontal axis. Not only is the reduction in variation apparent, but the model-based estimates at 60 yards are shifted downward. The explanation for this shift is lack-of-fit of the model for very long attempts. Shrinkage towards the empirical frequency tends to alleviate bias caused by this poor fit.

Fig.7

Model-based maximum-likelihood and shrinkage (midpoint) estimates of success probability, colored according to attempt distance (30, 40, 50, and 60 yards plotted with black, red,green and blue lines, respectively).

3.2Assessment using cross-validation

While HL is useful for quantifying how well several fitted models align with observed data, it may not tell the entire story as to how well the fitted model would generalize to observations made in the future, or out-of-sample data. Cross-validation (CV) can be used to quantify the degree to which a model fitted to observed data generalizes to data not used to fit the model or to compare the capacity of several models to generalize. Pasteur and Cunningham-Rhoads (2014) uses five-fold CV to select weather, pressure, and fatigue variables for a multiple logistic regression model. In five-fold CV, the dataset is partitioned into five subsets of equal size. The first subset is held out of the process where the model is fit, and then the fitted model is used to estimate the probabilities of the events that were held out. HL is then computed only for the held-out subset. This step is repeated four more times, holding out each subset exactly once. As can be seen from inspection of the first two numerical columns of Table 6, neither the model-based estimate (CLL) nor the empirical frequency (EF) generalize well under 5-fold CV, in the sense that there is disagreement between the number of observed and expected successful kicks when attempts are put into the 10 groups for the HL statistic. Recall that the midpoint estimated probability for any kick is simply the midpoint between CLL and EF. While this shrinkage estimate has a greater HL than its CLL and EF components when the entire dataset is used both to fit the model and to compute the goodness-of-fit statistic, it shows better agreement between observed and expected successful kicks under CV. The CV procedure was repeated four times, using four random partitions of the dataset into five subsets. Each time, the midpoint estimate enjoys a better fit for the held-out data than either of its two components, as quantified by the HL statistic. The shrinkage estimators based on the exponential weights, with a = 660,  670 or 680 also performed better than either component by itself. In computation of these estimates, the step of tuning the a “parameter” was not repeated during cross-validation. Rather, estimates based on values of a = 660,  670,  680 (which did the best on the entire dataset) were computed for every observation, and only those results are reported here. If this weighting scheme is tolerable, then it appears to perform comparably to the simple midpoint shrinkage estimator. Among the six estimators reported in Table 6, there is no dominant choice when considering performance under cross-validation. One clear observation, however, is that the midpoint and exponential shrinkage estimators are generalizing to out-of-sample data better than either the CLL or EF components by themselves.

4.1Ranking kickers

In addition to estimating the success probability for a given kicker from a given distance, the estimation methodology can be used to rank NFL kickers. Historically, comparisons among kickers in sports media have been made using the overall number (as on espn.go.com) or proportion of successful attempts (as on sports.yahoo.com or si.com) without regard to distance, as it not obvious how to incorporate distance and degree of difficulty based on summary statistics. Alternatively, Berry and Berry (1985) suggests averaging model-based kicker-specific estimates over the league-wide distribution of attempt distances. A histogram for these attempts is given in Fig. 8. The average distance over all 17, 104 attempts was 36.6 yards. Tables 7, 8, and 9 present a ranking of kickers under this summary statistic, using several different methods of estimation, along with the standard “overall” metric.

Fig.8

Frequency histogram of attempt distances.

For veteran kickers with a large number of kicks that are similar in distribution to that of the whole league, the effect of shrinkage is minimal. Adam Vinatieri has succeeded on 84% of 572 attempted kicks at an average distance of 35.8. His model-based estimate averaged over the distribution of all NFL kicks in 1998-2014 is 0.83 and the shrinkage estimate is only slightly different at 0.82. Conversely kickers like Patrick Murray and Greg Zuerlein have achieved very successful records so far in their young careers, hitting 20/24 (83%) and 73/89 (82%), at longer-than-average average distances of 41.1 yards and 40.0 yards respectively. Their model-based estimates averaged over the league-wide distribution of distances is more favorable (87% and 86%; note the bump up in rank from that based on the empirical frequency) but the shrinkage estimate brings these closer to more average kickers (84% for both kickers).

These tables suggest a ranking that places many current kickers among the best over the last 16 year period. Justin Tucker has been successful on 105 out of 116 attempts, with average distance d¯=38.1 for an overall mark of 91%. The kicker-specific, model-based estimate averaged over the league-wide distribution of attempt distances is 92%, but the midpoint technique shrinks this back to 87%. By any of the three methods, Tucker is ranked to be the most accurate among the K = 111 kickers underconsideration.

4.2Stadium effects

A topic of interest for many fans, writers and announcers is the effects of certain stadiums on the chances of making a field goal. The open end of Heinz Field receives particular attention for increased difficulty due to wind (Batista, 2002). Modeling fixed factorial effects for stadiums enables investigation of this issue. Though a likelihood ratio test finds evidence of stadium effects (χ² = 112, p < 0.0001, df = 50), their inclusion does not lead to an improved fit as assessed by Hosmer-Lemeshow, where there is a greater discrepancy between observed and expected counts than the model based on kickers and distances alone (HL = 8.3with stadium effects and HL = 1.9 without stadium effects, from Table 3). Nevertheless, including these effects does help to identify the more difficult and more favorable stadiums in which to kick, though again the estimation may suffer from selectionbias.

Success probabilities for each stadium can be estimated by back-transforming the marginal mean on the CLL scale. These marginal means average equally over all kickers and correspond to the average distance of d¯=36.6 yards. These estimates provide a ranking of stadium difficulty and are given in Table 11 along with empirical game and field goal attempt counts. We make no attempt to shrink these estimates, rather, we simply restrict our attention to stadiums for which more than 100 attempts are available in our dataset. Table 11 also gives average distances at which coaches have allowed their kickers to attempt field goals, game counts and empirical successrates.

The ranking in Table 11 seems to corroborate the anecdote that Heinz Field is a tough place to kick, as does the lower-than-average distance of field goal attempts. The new stadium in Denver, “Sports Authority Field at Mile High,” also has a reputation as a good place to kick and it does turn up high on this list. The estimates based on marginal means also tend to be lower than the empirical success rates. This observation may possibly be due to the fact that marginal means average equally over all kickers. Empirical success rates will tend to involve more attempts from high-frequency kickers, who tend to be better kickers.

4.3Big legs

The capacity to kick long-distance field goals is of particular interest to many football fans as well as NFL team management. One way to enrich the model to allow rankings among kickers to vary across distance is to expand the degree of the polynomial to quadratic. Using L and Q superscripts for the linear and quadratic kicker effects leads to the following representations of the generalized linearmodel:

For this investigation, attention is restricted to the K′ = 54 veteran kickers with at least 100 attempts in their careers, since small sample sizes present greater separability issues for fitting quadratic models. A comparison of the model with kicker-specific slopes nested within the quadratic model is not significant using this restricted dataset (χ² = 63.3, df = 54, p = 0.181). This non-significance could possibly be explained by selection bias, as kickers with many attempts will all tend to be good enough to remain in the league. While not significant, there is a preponderance of kickers for whom the estimated quadratic coefficients have small p-values, suggesting that there are some kickers for whom quadratic coefficients are non-zero. In fact the average p-value for these 54 tests of the form βkQ=0 is 0.16 providing anecdotal evidence of curvature. For one further bit of evidence, a comparison of nested models based on the logit link function is weakly significant, (χ² = 72.8, df = 54, p = 0.0450), though the fit is not as good as that provided bythe CLL.

Table 10 enables simultaneous inspection of kickers who have estimated success probabilities that rank highly (top 20) in at least one of the three distances, 35, 45 or 55 yards. Entries in the table are sorted by the success probability estimate at 55 yards, πˆk(55). Estimates at d = 35, 45 are also given along with ranks which enable identification of kickers who seem to have strengths at either medium or long distances (or both). While some kickers are accurate from any distance, for example current kickersBailey, Tucker, Walsh, and Barth, others have a decent shot at making long kicks but can be less accurate than is desired at short distances, for example Edinger and Hollis. Lastly are the kickers who rarely miss from intermediate distances, but are not estimated to have a big enough leg to have much more than a coin flip chance of success on long field goals, for example Gould, Kasay, and Stover. For comparison with observed data, the last three columns give empirical success frequencies after binning kicks into 10-yard intervals surrounding d = 35,  45,  55 yards. Clustering can be used to group kickers, based on their trivariate estimated success probability at the three distances, 35,45 and 55 yards, such that within groups, kickers are similar. The default hierarchical clustering algorithm implemented by SAS PROC CLUSTER suggests the existence of three clusters, and that cluster assignment is mostly determined by πˆ(55). The first 8 kickers in Table 9 form a “big leg” cluster, the last four kickers form the “accurate but no leg” cluster, and the remainder fall somewhere in between. The means of the estimated success probabilities for these three groups are plotted against yards in Fig. 6. This does suggest that there may be somewhat of a trade-off between a “big leg” and accuracy at closer ranges.

Fig.6

Three clusters of kickers using shrinkage estimates of success probabilities at three distances.

Table 13 helps quantify the variance-reducing property of shrinkage estimation by comparison of the coefficient of variation among kickers with maximum likelihood (no shrinkage) estimates. The table gives results for both linear and quadratic fits, with diminished coefficient of variation for shrinkage estimates, particularly for shorter attempts.

Another way to rank kickers is to consider big-legs and stadium effects and shrinkage estimation all together. Stadium effects can bias a ranking and by including them in the model, together with quadratic distance effects, one can attempt to rank kickers on an even footing. The model is then given by

4.4Non-stationarity

The analyses in this paper cover a wide time span. Kicker ability appears to be changing, especially in the last 10 years, with improving accuracy and experience among kickers (Clark et al., 2013). The frequencies of long field goal attempts are also increasing. Figure 9 shows the log of attempt frequencies against year, with separate plotting characters and colors for each distance, rounded to the nearest five yards. The log-scale is better for visualizing the increase in frequencies. In fact, the log of the frequency of kicks between 53 and 57 yards, denoted in the graph as “55”, and also of kicks between 48 and 52 yards, denoted as “50”, appear to be linearly increasing with season. The figure overlays lines using estimated coefficients from Poisson regressions of frequency on year, separate for each 5-yard distance interval, with the natural log-link function for the generalized linear model. The slopes are significant for distances of 50, 55, and 60 yards.

Fig.9

Log of attempt frequencies over season separated by distances rounded to the nearest five yards starting at 40 yards. Lines are linear regressions of log-transformed attempt frequencies on year.