Reverse Engineering Fantasy Points: Frequentist Method with Generalised Linear Models

Setting up the Formula

Formulas in R are set up with the following syntax: <response> ~ <variable1> + <variable2> + .... The response in this case is fp and variables (or predictors) are the stats.

In this example we won’t be investigating higher order interactions between terms.

formula <- "fp ~ kicks + handballs + marks + tackles + inside_50_s + bounces + frees_for + 
frees_against + one_percenters + hit_outs + goals + behinds + goal_assists"

Choosing a (Distribution) Family

Next step in setting up a glm is to specify the distribution family, and a good way to choose one is to find out what domain the response exists in. This requires some knowledge of the data set and using a density plot can aid us in making a decision. Common distributions include:

• Gaussian (normal): Continuous unbounded (positive and negative real numbers)
• Gamma: Continuous positive
• Poisson: Integer positive
• Logistic: Continuous between 0 and 1 (or exclusively 0 and 1)

library(ggplot2)
ggplot(afltables_2019_fp, aes(x = fp)) +
  geom_density() +
  geom_rug() +
  ggtitle("Distribution of Fantasy Points") + ylab("") + xlab("FP")

It’s important to note that FP can be negative, therefore FP is continuous and unbounded. With this knowledge, it makes sense to choose the normal (Gaussian) distribution.

Selecting a Link

We also need to select a link, and have 3 options with the Gaussian family:

1. Identity: An increase in a stat corresponds to an linear increase in FP
2. Log: An increase in a stat corresponds to an proportional increase in FP
3. Inverse: An increase in a stat corresponds to an inverse increase in FP

We are going to be using the identity link, mainly because it’s the easiest to interpret and falls in line with standard linear regressions.

Running the Model + Interpretation

Putting it all together, and we can check the model summary.

glm1 <- glm(
  formula = formula,
  data = afltables_2019_fp,
  family = gaussian(link = "identity")
)

summary(glm1)

## 
## Call:
## glm(formula = formula, family = gaussian(link = "identity"), 
##     data = afltables_2019_fp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -23.8761   -2.0752    0.0152    2.0318   25.9347  
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     2.0892003  0.1099437  19.002   <2e-16 ***
## kicks           3.0123931  0.0100210 300.608   <2e-16 ***
## handballs       1.9856901  0.0084920 233.830   <2e-16 ***
## marks           2.9898337  0.0164047 182.254   <2e-16 ***
## tackles         4.0360824  0.0161849 249.373   <2e-16 ***
## inside_50_s    -0.0506948  0.0217571  -2.330   0.0198 *  
## bounces         0.0411818  0.0462881   0.890   0.3737    
## frees_for       0.9838578  0.0349754  28.130   <2e-16 ***
## frees_against  -3.0309144  0.0337349 -89.845   <2e-16 ***
## one_percenters  0.0000577  0.0147666   0.004   0.9969    
## hit_outs        0.9948819  0.0049551 200.778   <2e-16 ***
## goals           6.0652085  0.0381539 158.967   <2e-16 ***
## behinds         0.9645699  0.0491625  19.620   <2e-16 ***
## goal_assists   -0.0034078  0.0537769  -0.063   0.9495    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 9.806355)
## 
##     Null deviance: 6561783  on 9107  degrees of freedom
## Residual deviance:   89179  on 9094  degrees of freedom
## AIC: 46657
## 
## Number of Fisher Scoring iterations: 2

The model has been run and the results are being displayed above. Let’s examine it one step at a time.

Call shows us what model has been run. Deviance Residuals give a quick insight into the distribution of the residuals, which should be normally distributed around 0. Looks good so far, median \(\approx\) 0, and the quartiles and min/max are approximately symmetrical. We’ll be generating residuals plots to better understand it later.

Coefficients contains the information regarding estimates for the intercept and each stat. The estimate column (when using the “identity” link) shows the mean effect of each term in the model. For example, an increase by one (unit of) tackle increases FP by an estimate of 4.03.

Std. Error is the error term associated with the term’s estimate. Low standard error means the model in more certain that the mean estimate. The Pr(>|t|) is the p-value, a calculation based off the t-value and the degrees of freedom in the model.

The p-value is a uniquely frequentist form of test set up to determine whether the estimated mean has a statistically significant difference to 0 in the form of a null hypothesis. It can be described by following equation:

\[ H_0:\mu = 0 \]

where \(\mu\) is the term (variable). Normally a p-value of 0.05 is tested against, with a p-value greater than 0.05 meaning we “Accept the null hypothesis” (that the estimate is equal to 0, therefore has no effect). If the p-value is less than 0.05, then we “Fail to reject the null hypothesis” (that the estimate is not equal to 0, and has an effect).

(Intercept) has an estimate of 2.09, and is statistically significant (2e-16 is less than 0.05). We can use the Std. Error (Standard Error (SE)) to calculate a confidence interval, which is much easier to interpret (we use the constant 1.96 due to using a 95%CI, a 99%CI would use 2.58):

\[ 95\% \texttt{ Confidence Interval}: \texttt{Estimate}\pm1.96*\texttt{SE}\\ \texttt{Intercept} = 2.09(95\%CI[1.87, 2.3]) \]

According to the model, we can be 95% sure that the intercept is between 1.87 and 2.3. The intercept is interpreted as the shift applied the entire dataset.

Looking at the other predictive variables, we can see that 10 of the 13 variables have p-values below the threshold of 0.05. You’ll notice that the variables with a p-value above the threshold have estimates close to (and overlapping with 95%CI) 0. For example, every bounce adds approximately 0 to FP, while kicks adds 3.012 (95%CI[2.99,3.03]) per kick to total FP.

par(mfrow=c(2,2))
plot(glm1)

The residuals vs fitted plot shows a nice linear constant distribution of residuals. If the variance increased (fanned out) we would have heteroscedasticity, and if there’s curvature then our model is most likely missing a squared (or nonlinear) term.

Outlier Detection

The QQ plot is quite nice, majority of the points fall on the dashed line, with the exception of a few outliers, which we will now remove. Using Cook’s Distance, we can see how individual observations effect the model and will remove the points that have the biggest effect. The bigger the value, the more influence that point has on the model.

par(mfrow=c(1,1))
plot(cooks.distance(glm1))

There appears to be some outliers in the dataset, we can remove about 30 data points (30/9000 or 0.3% of points). We are going to rerun the model with these outliers removed.

afltables_2019_fp_out_rem <- afltables_2019_fp %>% 
  mutate(
    id = row_number(),
    cooks_d = cooks.distance(glm1)
  ) %>% 
  arrange(-cooks_d) %>% # assign each obs its Cooks D
  slice(-(1:30)) %>%    # remove the top 30
  arrange(id)

# rerun the model without outliers
glm2 <- glm(
  formula = formula,
  data = afltables_2019_fp_out_rem, 
  family = gaussian(link = 'identity')
)

summary(glm2)

## 
## Call:
## glm(formula = formula, family = gaussian(link = "identity"), 
##     data = afltables_2019_fp_out_rem)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.3848   -2.0720    0.0121    2.0219   16.8099  
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     2.069185   0.106339  19.458   <2e-16 ***
## kicks           3.014910   0.009671 311.741   <2e-16 ***
## handballs       1.987187   0.008196 242.472   <2e-16 ***
## marks           2.983426   0.015866 188.044   <2e-16 ***
## tackles         4.033991   0.015660 257.604   <2e-16 ***
## inside_50_s    -0.054049   0.021004  -2.573   0.0101 *  
## bounces         0.055426   0.045887   1.208   0.2271    
## frees_for       1.000116   0.033831  29.562   <2e-16 ***
## frees_against  -3.040852   0.032581 -93.331   <2e-16 ***
## one_percenters  0.007380   0.014268   0.517   0.6050    
## hit_outs        0.996866   0.004812 207.166   <2e-16 ***
## goals           6.077538   0.036991 164.296   <2e-16 ***
## behinds         0.975389   0.047510  20.530   <2e-16 ***
## goal_assists   -0.009567   0.051917  -0.184   0.8538    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 9.103323)
## 
##     Null deviance: 6491380  on 9077  degrees of freedom
## Residual deviance:   82513  on 9064  degrees of freedom
## AIC: 45828
## 
## Number of Fisher Scoring iterations: 2

par(mfrow=c(2,2))
plot(glm2)

QQ plot tails are nicer.

par(mfrow=c(1,1))
plot(cooks.distance(glm2))

No individual points really stand out, no more outliers we can confidently remove.

Remove Non-informative Variables

The next step is to remove variables that have no predictive power on the model, and we can calculate this using the AIC metric. Akaike’s Information Criterion (AIC) is a penalisation term for models. Essentially, lower AIC is better, and we can use the following code to perform this automatically. The process may look a lot is going on, but I’ll do my best to break it down.

glm2_step <- MASS::stepAIC(glm2, direction = 'both')

## Start:  AIC=45828.27
## fp ~ kicks + handballs + marks + tackles + inside_50_s + bounces + 
##     frees_for + frees_against + one_percenters + hit_outs + goals + 
##     behinds + goal_assists
## 
##                  Df Deviance   AIC
## - goal_assists    1    82513 45826
## - one_percenters  1    82515 45827
## - bounces         1    82526 45828
## <none>                 82513 45828
## - inside_50_s     1    82573 45833
## - behinds         1    86350 46239
## - frees_for       1    90468 46662
## - frees_against   1   161809 51940
## - goals           1   328241 58361
## - marks           1   404412 60256
## - hit_outs        1   473206 61682
## - handballs       1   617719 64101
## - tackles         1   686607 65061
## - kicks           1   967196 68171
## 
## Step:  AIC=45826.3
## fp ~ kicks + handballs + marks + tackles + inside_50_s + bounces + 
##     frees_for + frees_against + one_percenters + hit_outs + goals + 
##     behinds
## 
##                  Df Deviance   AIC
## - one_percenters  1    82515 45825
## - bounces         1    82526 45826
## <none>                 82513 45826
## + goal_assists    1    82513 45828
## - inside_50_s     1    82579 45832
## - behinds         1    86353 46237
## - frees_for       1    90473 46660
## - frees_against   1   161814 51938
## - goals           1   330714 58427
## - marks           1   404981 60266
## - hit_outs        1   473494 61685
## - handballs       1   619111 64119
## - tackles         1   686793 65061
## - kicks           1   970731 68202
## 
## Step:  AIC=45824.58
## fp ~ kicks + handballs + marks + tackles + inside_50_s + bounces + 
##     frees_for + frees_against + hit_outs + goals + behinds
## 
##                  Df Deviance   AIC
## - bounces         1    82528 45824
## <none>                 82515 45825
## + one_percenters  1    82513 45826
## + goal_assists    1    82515 45827
## - inside_50_s     1    82588 45831
## - behinds         1    86395 46240
## - frees_for       1    90519 46663
## - frees_against   1   161887 51940
## - goals           1   338089 58626
## - marks           1   412895 60440
## - hit_outs        1   477291 61756
## - handballs       1   628707 64257
## - tackles         1   687665 65071
## - kicks           1   970732 68200
## 
## Step:  AIC=45824
## fp ~ kicks + handballs + marks + tackles + inside_50_s + frees_for + 
##     frees_against + hit_outs + goals + behinds
## 
##                  Df Deviance   AIC
## <none>                 82528 45824
## + bounces         1    82515 45825
## + one_percenters  1    82526 45826
## + goal_assists    1    82528 45826
## - inside_50_s     1    82598 45830
## - behinds         1    86403 46239
## - frees_for       1    90520 46661
## - frees_against   1   161972 51943
## - goals           1   338123 58624
## - marks           1   413917 60461
## - hit_outs        1   477681 61761
## - handballs       1   629399 64265
## - tackles         1   689561 65094
## - kicks           1   991022 68386

We start off with AIC = 45828.27, and the stepwise AIC calculates that removing goal_assists reduces the AIC, and the new model has an AIC = 45826.3, thus improving the model. This process is repeated until removing variables no longer decreases the AIC. Using direction = 'both' allows for the algorithm to check if re-adding previously removed variables increases the AIC. In this case, it doesn’t.

Let’s check out the summary for our new model.

summary(glm2_step)

## 
## Call:
## glm(formula = fp ~ kicks + handballs + marks + tackles + inside_50_s + 
##     frees_for + frees_against + hit_outs + goals + behinds, family = gaussian(link = "identity"), 
##     data = afltables_2019_fp_out_rem)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.3897   -2.0712    0.0067    2.0277   16.7968  
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.093698   0.096962  21.593  < 2e-16 ***
## kicks          3.016683   0.009549 315.930  < 2e-16 ***
## handballs      1.986821   0.008106 245.117  < 2e-16 ***
## marks          2.983426   0.015636 190.809  < 2e-16 ***
## tackles        4.032478   0.015615 258.248  < 2e-16 ***
## inside_50_s   -0.055593   0.020097  -2.766  0.00568 ** 
## frees_for      0.999574   0.033735  29.630  < 2e-16 ***
## frees_against -3.041033   0.032551 -93.424  < 2e-16 ***
## hit_outs       0.996921   0.004785 208.359  < 2e-16 ***
## goals          6.073662   0.036245 167.574  < 2e-16 ***
## behinds        0.971072   0.047064  20.633  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 9.102053)
## 
##     Null deviance: 6491380  on 9077  degrees of freedom
## Residual deviance:   82528  on 9067  degrees of freedom
## AIC: 45824
## 
## Number of Fisher Scoring iterations: 2

The model suggests that inside_50_s is a statistically significant variable, but on closer inspection it seems something weird is going on. Each inside 50 adds -0.056 (95%CI[-0.095,-0.017]) to the FP score, which begs the question: does FP penalise for inside 50’s? Using our intuition it seems much more likely that the true effect is 0, and can reasonably argued for it to be excluded from the model.

formula_3 <- "fp ~ kicks + handballs + marks + tackles + frees_for + 
frees_against + hit_outs + goals + behinds"

glm3 <- glm(
  formula = formula_3,
  data = afltables_2019_fp_out_rem, 
  family = gaussian(link = 'identity')
)

summary(glm3)

## 
## Call:
## glm(formula = formula_3, family = gaussian(link = "identity"), 
##     data = afltables_2019_fp_out_rem)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.4164   -2.0786    0.0108    2.0108   16.8300  
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.107037   0.096878   21.75   <2e-16 ***
## kicks          3.004733   0.008519  352.72   <2e-16 ***
## handballs      1.982782   0.007976  248.59   <2e-16 ***
## marks          2.991265   0.015382  194.46   <2e-16 ***
## tackles        4.027332   0.015509  259.67   <2e-16 ***
## frees_for      1.000197   0.033746   29.64   <2e-16 ***
## frees_against -3.043979   0.032545  -93.53   <2e-16 ***
## hit_outs       0.996308   0.004781  208.38   <2e-16 ***
## goals          6.062463   0.036031  168.26   <2e-16 ***
## behinds        0.953944   0.046672   20.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 9.10873)
## 
##     Null deviance: 6491380  on 9077  degrees of freedom
## Residual deviance:   82598  on 9068  degrees of freedom
## AIC: 45830
## 
## Number of Fisher Scoring iterations: 2

Interestingly, the AIC did increase but the residual deviance decreased.

Presenting Results

Next step is to convey the estimates and confidence intervals.

library(broom)
glm3_coef <- tidy(glm3)
glm3_confint <- glm3_coef %>% 
  mutate(
    lb = estimate-1.96*std.error,
    ub = estimate+1.96*std.error,
  )

library(ggplot2)
ggplot(glm3_confint, aes(x=estimate, y=reorder(term, desc(term)))) + 
  geom_errorbar(aes(xmin=lb, xmax=ub), width=.3) +
  geom_point() +
  ggtitle("The Effect of AFL Statistics on Total Fantasy Points") + 
  xlab("Estimate (95% CI)") + ylab("Term")

We are 95% confident that the true mean of the term is between these confidence intervals.

glm3_confint_clean <- 
  glm3_confint %>% 
  mutate(
    estimate = round(estimate, 2),
    lb = round(lb, 2),
    ub = round(ub, 2)
  ) %>% 
  select(Term = term, Estimate = estimate, `Lower Bound` = lb, `Upper Bound` = ub)

DT::datatable(glm3_confint_clean, rownames = F, width = 500,
              caption = "The Effect of AFL Statistics on Total Fantasy Points (95%CI)", 
              options = list(dom = 't'))