I need your help!

I want your feedback to make the book better for you and other readers. If you find typos, errors, or places where the text may be improved, please let me know. The best ways to provide feedback are by GitHub or hypothes.is annotations.

You can leave a comment at the bottom of the page/chapter, or open an issue or submit a pull request on GitHub: https://github.com/isaactpetersen/Fantasy-Football-Analytics-Textbook

Hypothesis Alternatively, you can leave an annotation using hypothes.is. To add an annotation, select some text and then click the symbol on the pop-up menu. To see the annotations of others, click the symbol in the upper right-hand corner of the page.

25  Time Series Analysis

This chapter provides an overview of time series analysis.

25.1 Getting Started

25.1.1 Load Packages

Code
library("petersenlab")
library("xts")
library("zoo")
library("forecast")
library("brms")
library("rstan")
library("plotly")
library("tidyverse")

25.1.2 Load Data

Code
load(file = "./data/player_stats_weekly.RData")
load(file = "./data/player_stats_seasonal.RData")

We created the player_stats_weekly.RData and player_stats_seasonal.RData objects in Section 4.4.3. The following code loads the Bayesian model object that was fit in Section 12.4.5.

Code
load(url("https://osf.io/download/q6rjf/")) # Bayesian model object

25.2 Overview of Time Series Analysis

Time series analysis is useful when trying to generate forecasts from longitudinal data. That is, time series analysis seeks to evaluate change over time to predict future values.

There are many different types of time series analyses. For simplicity, in this chapter, we use autoregressive integrated moving average (ARIMA) and exponential smoothing models to demonstrate one approach to time series analysis. We also leverage Bayesian mixed models to generate forecasts of future performance and plots of individuals model-implied performance by age and position.

There are four (potential) key components of time series data:

  • seasonality
  • trend
  • cycle
  • irregularity

Seasonality corresponds to fluctuations that occur at fixed and known periods (e.g., time of day, day of week, or month of year) (Hyndman & Athanasopoulos, 2021). A trend represents a long-term change (i.e., increase or decrease) in the series. A cycle reflects one or more rises and falls that are not of a fixed frequency, with durations typically lasting at least two years (Hyndman & Athanasopoulos, 2021). Irregularity is random variation or noise.

25.3 Autoregressive Integrated Moving Average (ARIMA) Models

Hyndman & Athanasopoulos (2021) provide a nice overview of ARIMA models. As noted by Hyndman & Athanasopoulos (2021), ARIMA models aim to describe how a variable is correlated with itself over time (autocorrelation)—i.e., how earlier levels of a variable are correlated with later levels of the same variable. ARIMA models perform best when there is a clear pattern where later values are influenced by earlier values. ARIMA models incorporate autoregression effects, moving average effects, and differencing.

ARIMA models can have various numbers of terms and model complexity. They are specified in the following form: \(\text{ARIMA}(p,d,q)\), where:

  • \(p =\) the number of autoregressive terms
  • \(d =\) the number of differences between consecutive scores (to make the time series stationary by reducing trends and seasonality)
  • \(q =\) the number of moving average terms

ARIMA models assume that the data are stationary (i.e., there are no long-term trends), are non-seasonal (i.e., there is no consistency of the timing of the peaks or troughs in the line), and that earlier values influence later values. This may not strongly be the case in fantasy football, so ARIMA models may not be particularly useful in forecasting fantasy football performance. Other approaches, such as exponential smoothing, may be useful for data that show longer-term trends and seasonality (Hyndman & Athanasopoulos, 2021). Nevertheless, ARIMA models are widely used in forecasting financial markets and economic indicators. Thus, it is a useful technique to learn.

Adapted from: https://rc2e.com/timeseriesanalysis (Long & Teetor, 2019; archived at https://perma.cc/U5P6-2VWC).

25.3.1 Create the Time Series Objects

Code
weeklyFantasyPoints_tomBrady <- player_stats_weekly |> 
  filter(
    player_id == "00-0019596" | player_display_name == "Tom Brady")

weeklyFantasyPoints_peytonManning <- player_stats_weekly |> 
  filter(
    player_id == "00-0010346" | player_display_name == "Peyton Manning")

ts_tomBrady <- xts::xts(
  x = weeklyFantasyPoints_tomBrady["fantasyPoints"],
  order.by = weeklyFantasyPoints_tomBrady$gameday)

ts_peytonManning <- xts::xts(
  x = weeklyFantasyPoints_peytonManning["fantasyPoints"],
  order.by = weeklyFantasyPoints_peytonManning$gameday)

ts_tomBrady
           fantasyPoints
2000-11-23          0.24
2001-09-23          2.74
2001-09-30          6.92
2001-10-07          4.34
2001-10-14         22.56
2001-10-21         19.88
2001-10-28         10.02
2001-11-04         22.00
2001-11-11          8.18
2001-11-18          9.00
       ...              
2022-10-27         17.10
2022-11-06         15.20
2022-11-13         17.02
2022-11-27         18.04
2022-12-05         17.14
2022-12-11         10.12
2022-12-18         20.58
2022-12-25         11.34
2023-01-01         37.68
2023-01-08          7.36
Code
ts_peytonManning
           fantasyPoints
1999-09-12         15.06
1999-09-19         17.22
1999-09-26         29.56
1999-10-10         20.66
1999-10-17         10.10
1999-10-24         17.86
1999-10-31         18.52
1999-11-07         20.60
1999-11-14         15.18
1999-11-21         22.80
       ...              
2015-09-13          4.90
2015-09-17         20.24
2015-09-27         18.86
2015-10-04          8.32
2015-10-11          6.64
2015-10-18          9.60
2015-11-01         11.60
2015-11-08         15.24
2015-11-15         -6.60
2016-01-03          2.56
Code
ts_combined <- merge(
  ts_tomBrady,
  ts_peytonManning
)

names(ts_combined) <- c("Tom Brady","Peyton Manning")

25.3.2 Plot the Time Series

Code
plot(
  ts_tomBrady,
  main = "Tom Brady's Fantasy Points by Game")
Tom Brady's Historical Fantasy Points by Game.
Figure 25.1: Tom Brady’s Historical Fantasy Points by Game.
Code
plot(
  ts_combined,
  legend,
  legend.loc = "topright",
  main = "Fantasy Points by Game")
Historical Fantasy Points by Game for Tom Brady and Peyton Manning.
Figure 25.2: Historical Fantasy Points by Game for Tom Brady and Peyton Manning.

25.3.3 Rolling Mean/Median

Code
zoo::rollmean(
  x = ts_tomBrady,
  k = 5)
           fantasyPoints
2001-09-30         7.360
2001-10-07        11.288
2001-10-14        12.744
2001-10-21        15.760
2001-10-28        16.528
2001-11-04        13.816
2001-11-11        15.384
2001-11-18        15.084
2001-11-25        11.568
2001-12-02        11.888
       ...              
2022-10-16        18.332
2022-10-23        15.892
2022-10-27        15.348
2022-11-06        16.212
2022-11-13        16.900
2022-11-27        15.504
2022-12-05        16.580
2022-12-11        15.444
2022-12-18        19.372
2022-12-25        17.416
Code
zoo::rollmedian(
  x = ts_tomBrady,
  k = 5)
           fantasyPoints
2001-09-30          4.34
2001-10-07          6.92
2001-10-14         10.02
2001-10-21         19.88
2001-10-28         19.88
2001-11-04         10.02
2001-11-11         10.02
2001-11-18          9.00
2001-11-25          8.52
2001-12-02          9.00
       ...              
2022-10-16         17.10
2022-10-23         15.20
2022-10-27         15.20
2022-11-06         17.02
2022-11-13         17.10
2022-11-27         17.02
2022-12-05         17.14
2022-12-11         17.14
2022-12-18         17.14
2022-12-25         11.34

25.3.4 Autocorrelation

The autocorrelation function (ACF) plot depicts the autocorrelation of scores as a function of the length of the lag. We can generate an ACF plot using the stats::acf() function. Significant autocorrelation is detected when the autocorrelation exceeds the dashed blue lines, as is depicted in Figure 25.3.

Code
acf(ts_tomBrady)
Autocorrelation Function (ACF) Plot of Tom Brady's Historical Fantasy Points by Game.
Figure 25.3: Autocorrelation Function (ACF) Plot of Tom Brady’s Historical Fantasy Points by Game.
Code
Box.test(ts_tomBrady)

    Box-Pierce test

data:  ts_tomBrady
X-squared = 4.7241, df = 1, p-value = 0.02974

25.3.5 Fit an Autoregressive Integrated Moving Average Model

We can fit an ARIMA model using the forecast::auto.arima() function of the forecast package (Hyndman et al., 2024; Hyndman & Khandakar, 2008):

Code
forecast::auto.arima(ts_tomBrady)
Series: ts_tomBrady 
ARIMA(2,1,2) 

Coefficients:
         ar1     ar2      ma1     ma2
      0.5777  0.0710  -1.4986  0.5109
s.e.  0.2627  0.0658   0.2590  0.2508

sigma^2 = 62.73:  log likelihood = -1164.34
AIC=2338.68   AICc=2338.86   BIC=2357.73
Code
forecast::auto.arima(ts_peytonManning)
Series: ts_peytonManning 
ARIMA(1,0,1) with non-zero mean 

Coefficients:
         ar1      ma1     mean
      0.8659  -0.7314  18.1494
s.e.  0.0970   0.1281   0.9568

sigma^2 = 58:  log likelihood = -860.86
AIC=1729.72   AICc=1729.89   BIC=1743.81

We can generate a plot of an ARIMA model using the stats::arima() function.

Code
arima_tomBrady <- arima(
  ts_tomBrady,
  order = c(5, 1, 4))

summary(arima_tomBrady)

Call:
arima(x = ts_tomBrady, order = c(5, 1, 4))

Coefficients:
         ar1      ar2      ar3      ar4     ar5      ma1     ma2     ma3
      0.6946  -0.0921  -0.4413  -0.1324  0.2165  -1.6067  0.7292  0.4654
s.e.  0.1531   0.2410   0.1556   0.0713  0.0686   0.1513  0.3245  0.3163
          ma4
      -0.5617
s.e.   0.1284

sigma^2 estimated as 59.47:  log likelihood = -1157.56,  aic = 2335.13

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE      MASE
Training set 0.6849811 7.700254 6.250196 -21.00147 50.53162 0.7358044
                    ACF1
Training set 0.002983661
Code
confint(arima_tomBrady)
          2.5 %       97.5 %
ar1  0.39449821  0.994738991
ar2 -0.56445714  0.380292004
ar3 -0.74622928 -0.136340954
ar4 -0.27206259  0.007285576
ar5  0.08202734  0.350885782
ma1 -1.90316251 -1.310212551
ma2  0.09317143  1.365135540
ma3 -0.15444207  1.085273172
ma4 -0.81325437 -0.310090232
Code
forecast::checkresiduals(arima_tomBrady)

    Ljung-Box test

data:  Residuals from ARIMA(5,1,4)
Q* = 3.7477, df = 3, p-value = 0.29

Model df: 9.   Total lags used: 12
Model Summary of Autoregressive Integrated Moving Average Model fit to Tom Brady's Historical Performance by Game.
Figure 25.4: Model Summary of Autoregressive Integrated Moving Average Model fit to Tom Brady’s Historical Performance by Game.
Code
arima_tomBrady_removeNonSigTerms <- arima(
  ts_tomBrady,
  order = c(5, 1, 4),
  fixed = c(NA, NA, 0, NA, NA, NA, NA, NA, NA))

summary(arima_tomBrady_removeNonSigTerms)

Call:
arima(x = ts_tomBrady, order = c(5, 1, 4), fixed = c(NA, NA, 0, NA, NA, NA, 
    NA, NA, NA))

Coefficients:
          ar1      ar2  ar3     ar4     ar5     ma1      ma2      ma3      ma4
      -0.9126  -0.7976    0  0.1011  0.1551  0.0220  -0.0196  -0.6575  -0.2559
s.e.   0.3146   0.1756    0  0.1019  0.0936  0.3219   0.2174   0.1514   0.1086

sigma^2 estimated as 59.34:  log likelihood = -1157.17,  aic = 2332.34

Training set error measures:
                    ME     RMSE      MAE       MPE     MAPE      MASE
Training set 0.6609404 7.691756 6.166178 -20.94507 50.26368 0.7259134
                     ACF1
Training set -0.006499149
Code
confint(arima_tomBrady_removeNonSigTerms)
          2.5 %      97.5 %
ar1 -1.52919038 -0.29605302
ar2 -1.14182198 -0.45338814
ar3          NA          NA
ar4 -0.09850619  0.30079015
ar5 -0.02830776  0.33860502
ma1 -0.60894381  0.65298527
ma2 -0.44563802  0.40636324
ma3 -0.95418084 -0.36079739
ma4 -0.46874055 -0.04309774
Code
forecast::checkresiduals(arima_tomBrady_removeNonSigTerms)

    Ljung-Box test

data:  Residuals from ARIMA(5,1,4)
Q* = 3.448, df = 3, p-value = 0.3276

Model df: 9.   Total lags used: 12
Model Summary of modified Autoregressive Integrated Moving Average Model fit to Tom Brady's Historical Performance by Game.
Figure 25.5: Model Summary of modified Autoregressive Integrated Moving Average Model fit to Tom Brady’s Historical Performance by Game.
Code
arima_peytonManning <- arima(
  ts_peytonManning,
  order = c(1, 0, 1))

summary(arima_peytonManning)

Call:
arima(x = ts_peytonManning, order = c(1, 0, 1))

Coefficients:
         ar1      ma1  intercept
      0.8659  -0.7314    18.1494
s.e.  0.0970   0.1281     0.9568

sigma^2 estimated as 57.31:  log likelihood = -860.86,  aic = 1729.72

Training set error measures:
                      ME    RMSE     MAE     MPE     MAPE      MASE       ACF1
Training set 0.001963284 7.57018 5.93687 -101.59 126.2628 0.7576497 0.02314746
Code
confint(arima_peytonManning)
               2.5 %     97.5 %
ar1        0.6758287  1.0558937
ma1       -0.9824854 -0.4802354
intercept 16.2740870 20.0247285
Code
forecast::checkresiduals(arima_peytonManning)

    Ljung-Box test

data:  Residuals from ARIMA(1,0,1) with non-zero mean
Q* = 9.6606, df = 8, p-value = 0.2897

Model df: 2.   Total lags used: 10
Model Summary of Autoregressive Integrated Moving Average Model fit to Peyton Manning's Historical Performance by Game.
Figure 25.6: Model Summary of Autoregressive Integrated Moving Average Model fit to Peyton Manning’s Historical Performance by Game.

25.3.6 Generate the Model Forecasts

We can generate model forecasts from the ARIMA model using the forecast::forecast() function.

Code
forecast_tomBrady <- forecast::forecast(
  arima_tomBrady,
  level = c(80, 95)) # 80% and 95% confidence intervals

forecast_peytonManning <- forecast::forecast(
  arima_peytonManning,
  level = c(80, 95)) # 80% and 95% confidence intervals

forecast_tomBrady
    Point Forecast     Lo 80    Hi 80      Lo 95    Hi 95
336       17.53222  7.649182 27.41525  2.4174197 32.64701
337       22.03333 12.112159 31.95449  6.8602096 37.20644
338       14.24893  4.286668 24.21119 -0.9870347 29.48489
339       21.74409 11.583707 31.90447  6.2051255 37.28305
340       17.77125  7.610050 27.93246  2.2310342 33.31147
341       19.36256  9.074481 29.65064  3.6283007 35.09682
342       19.53110  9.226845 29.83535  3.7721023 35.29009
343       18.57753  8.261160 28.89390  2.8000019 34.35506
344       19.34576  9.025269 29.66626  3.5619297 35.12959
345       18.82221  8.501399 29.14301  3.0378949 34.60652
Code
forecast_peytonManning
    Point Forecast    Lo 80    Hi 80       Lo 95    Hi 95
251       12.03918 2.337602 21.74076 -2.79810196 26.87646
252       12.85880 3.069861 22.64773 -2.11208824 27.82968
253       13.56847 3.714551 23.42240 -1.50180107 28.63875
254       14.18295 4.280590 24.08532 -0.96140568 29.32732
255       14.71501 4.776482 24.65354 -0.48465672 29.91468
256       15.17570 5.210142 25.14125 -0.06530308 30.41669
257       15.57459 5.588819 25.56035  0.30267291 30.84650
258       15.91997 5.919075 25.92086  0.62492056 31.21502
259       16.21902 6.206802 26.23125  0.90665236 31.53140
260       16.47796 6.457258 26.49867  1.15261775 31.80331

25.3.7 Plot the Model Forecasts

We can plot the model forecasts using the forecast::autoplot() function.

Code
forecast::autoplot(forecast_tomBrady) + 
  labs(
    x = "Game Number",
    y = "Fantasy Points",
    title = "Tom Brady's Historical and Projected Fantasy Points by Game",
    subtitle = "(if he were to have continued playing additional seasons)"
  ) +
  theme_classic()
Tom Brady's Historical and Projected Fantasy Points by Game Based on Autoregressive Integrated Moving Average (ARIMA) Model.
Figure 25.7: Tom Brady’s Historical and Projected Fantasy Points by Game Based on Autoregressive Integrated Moving Average (ARIMA) Model.
Code
forecast::autoplot(forecast_peytonManning) + 
  labs(
    x = "Game Number",
    y = "Fantasy Points",
    title = "Peyton Manning's Historical and Projected Fantasy Points by Game",
    subtitle = "(if he were to have continued playing additional seasons)"
  ) +
  theme_classic()
Peyton Manning's Historical and Projected Fantasy Points by Game Based on Autoregressive Integrated Moving Average (ARIMA) Model.
Figure 25.8: Peyton Manning’s Historical and Projected Fantasy Points by Game Based on Autoregressive Integrated Moving Average (ARIMA) Model.

25.4 Exponential Smoothing

25.4.1 Fit the Exponential Smoothing Model

We fit exponential smoothing models using the forecast::ets() function of the forecast package (Hyndman et al., 2024; Hyndman & Khandakar, 2008):

Code
ets_tomBrady <- forecast::ets(ts_tomBrady)
ets_peytonManning <- forecast::ets(ts_peytonManning)

We can generate a plot of an ARIMA model using the stats::arima() function.

Code
summary(ets_tomBrady)
ETS(A,N,N) 

Call:
forecast::ets(y = ts_tomBrady)

  Smoothing parameters:
    alpha = 0.0427 

  Initial states:
    l = 13.5045 

  sigma:  7.9429

     AIC     AICc      BIC 
3340.152 3340.224 3351.594 

Training set error measures:
                    ME     RMSE     MAE       MPE     MAPE      MASE       ACF1
Training set 0.3720536 7.919124 6.37574 -42.07214 70.00621 0.7505841 0.04197058
Code
forecast::checkresiduals(ets_tomBrady)

    Ljung-Box test

data:  Residuals from ETS(A,N,N)
Q* = 19.555, df = 10, p-value = 0.03376

Model df: 0.   Total lags used: 10
Model Summary of Exponential Smoothing Model fit to Tom Brady's Historical Performance by Game.
Figure 25.9: Model Summary of Exponential Smoothing Model fit to Tom Brady’s Historical Performance by Game.
Code
summary(ets_peytonManning)
ETS(A,N,N) 

Call:
forecast::ets(y = ts_peytonManning)

  Smoothing parameters:
    alpha = 0.1374 

  Initial states:
    l = 18.0157 

  sigma:  7.7078

     AIC     AICc      BIC 
2405.471 2405.569 2416.036 

Training set error measures:
                     ME     RMSE      MAE       MPE     MAPE      MASE
Training set -0.2418402 7.676873 6.036916 -107.9422 132.1856 0.7704173
                   ACF1
Training set 0.04724071
Code
forecast::checkresiduals(ets_peytonManning)

    Ljung-Box test

data:  Residuals from ETS(A,N,N)
Q* = 10.9, df = 10, p-value = 0.3653

Model df: 0.   Total lags used: 10
Model Summary of Exponential Smoothing Model fit to Peyton Manning's Historical Performance by Game.
Figure 25.10: Model Summary of Exponential Smoothing Model fit to Peyton Manning’s Historical Performance by Game.

25.4.2 Generate the Model Forecasts

We can generate model forecasts from the ARIMA model using the forecast::forecast() function.

Code
forecastETS_tomBrady <- forecast::forecast(
  ets_tomBrady,
  level = c(80, 95)) # 80% and 95% confidence intervals

forecastETS_peytonManning <- forecast::forecast(
  ets_peytonManning,
  level = c(80, 95)) # 80% and 95% confidence intervals

forecastETS_tomBrady
    Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
336       18.82327 8.644072 29.00247 3.255530 34.39101
337       18.82327 8.634807 29.01173 3.241361 34.40518
338       18.82327 8.625551 29.02099 3.227206 34.41933
339       18.82327 8.616304 29.03023 3.213063 34.43348
340       18.82327 8.607065 29.03947 3.198933 34.44761
341       18.82327 8.597834 29.04870 3.184815 34.46172
342       18.82327 8.588611 29.05793 3.170711 34.47583
343       18.82327 8.579397 29.06714 3.156619 34.48992
344       18.82327 8.570191 29.07635 3.142540 34.50400
345       18.82327 8.560994 29.08554 3.128473 34.51807
Code
forecastETS_peytonManning
    Point Forecast      Lo 80    Hi 80     Lo 95    Hi 95
251       9.708678 -0.1692214 19.58658 -5.398265 24.81562
252       9.708678 -0.2620216 19.67938 -5.540191 24.95755
253       9.708678 -0.3539661 19.77132 -5.680808 25.09816
254       9.708678 -0.4450780 19.86243 -5.820151 25.23751
255       9.708678 -0.5353795 19.95273 -5.958256 25.37561
256       9.708678 -0.6248920 20.04225 -6.095153 25.51251
257       9.708678 -0.7136357 20.13099 -6.230875 25.64823
258       9.708678 -0.8016302 20.21899 -6.365451 25.78281
259       9.708678 -0.8888940 20.30625 -6.498909 25.91626
260       9.708678 -0.9754451 20.39280 -6.631278 26.04863

25.4.3 Plot the Model Forecasts

We can plot the model forecasts using the forecast::autoplot() function.

Code
forecast::autoplot(forecastETS_tomBrady) + 
  labs(
    x = "Game Number",
    y = "Fantasy Points",
    title = "Tom Brady's Historical and Projected Fantasy Points by Game",
    subtitle = "(if he were to have continued playing additional seasons)"
  ) +
  theme_classic()
Tom Brady's Historical and Projected Fantasy Points by Game Based on Exponential Smoothing.
Figure 25.11: Tom Brady’s Historical and Projected Fantasy Points by Game Based on Exponential Smoothing.
Code
forecast::autoplot(forecastETS_peytonManning) + 
  labs(
    x = "Game Number",
    y = "Fantasy Points",
    title = "Peyton Manning's Historical and Projected Fantasy Points by Game",
    subtitle = "(if he were to have continued playing additional seasons)"
  ) +
  theme_classic()
Peyton Manning's Historical and Projected Fantasy Points by Game Based on Exponential Smoothing.
Figure 25.12: Peyton Manning’s Historical and Projected Fantasy Points by Game Based on Exponential Smoothing.

25.5 Bayesian Mixed Models

The Bayesian longitudinal mixed models were estimated in Section 12.4.5.

25.5.1 Prepare New Data Object

Code
player_stats_seasonal_offense_subset <- player_stats_seasonal |> 
  dplyr::filter(position_group %in% c("QB","RB","WR","TE") | position %in% c("K"))

player_stats_seasonal_offense_subset$position[which(player_stats_seasonal_offense_subset$position == "HB")] <- "RB"

player_stats_seasonal_offense_subset$player_idFactor <- factor(player_stats_seasonal_offense_subset$player_id)
player_stats_seasonal_offense_subset$positionFactor <- factor(player_stats_seasonal_offense_subset$position)
Code
player_stats_seasonal_offense_subsetCC <- player_stats_seasonal_offense_subset |>
  filter(
    !is.na(player_idFactor),
    !is.na(fantasyPoints),
    !is.na(positionFactor),
    !is.na(ageCentered20),
    !is.na(ageCentered20Quadratic),
    !is.na(years_of_experience))

player_stats_seasonal_offense_subsetCC <- player_stats_seasonal_offense_subsetCC |> 
  filter(player_id %in% bayesianMixedModelFit$data$player_idFactor) |> 
  mutate(positionFactor = droplevels(positionFactor))

player_stats_seasonal_offense_subsetCC <- player_stats_seasonal_offense_subsetCC |>
  group_by(player_id) |> 
  group_modify(~ add_row(.x, season = max(player_stats_seasonal_offense_subsetCC$season) + 1)) |> 
  fill(player_display_name, player_idFactor, position, position_group, positionFactor, team, .direction = "downup") |> 
  ungroup()

player_stats_seasonal_offense_subsetCC <- player_stats_seasonal_offense_subsetCC |> 
  left_join(
    player_stats_seasonal_offense_subsetCC |> 
      filter(season == max(player_stats_seasonal_offense_subsetCC$season) - 1) |> 
      select(player_id, age_lastYear = age, years_of_experience_lastYear = years_of_experience),
    by = "player_id") |>
  mutate(
    age = if_else(season == max(player_stats_seasonal_offense_subsetCC$season), age_lastYear + 1, age), # increment age by 1
    ageCentered20 = age - 20,
    years_of_experience = if_else(season == max(player_stats_seasonal_offense_subsetCC$season), years_of_experience_lastYear + 1, years_of_experience)) # increment experience by 1

activePlayers <- unique(player_stats_seasonal_offense_subsetCC[c("player_id","season")]) |> 
  filter(season == max(player_stats_seasonal_offense_subsetCC$season) - 1) |> 
  select(player_id) |> 
  pull()

inactivePlayers <- player_stats_seasonal_offense_subsetCC$player_id[which(player_stats_seasonal_offense_subsetCC$player_id %ni% activePlayers)]

player_stats_seasonal_offense_subsetCC <- player_stats_seasonal_offense_subsetCC |> 
  filter(player_id %in% activePlayers | (player_id %in% inactivePlayers & season < max(player_stats_seasonal_offense_subsetCC$season) - 1)) |> 
  mutate(
    player_idFactor = droplevels(player_idFactor) 
  )

25.5.2 Generate Predictions

Code
player_stats_seasonal_offense_subsetCC$fantasyPoints_bayesian <- predict(
  bayesianMixedModelFit,
  newdata = player_stats_seasonal_offense_subsetCC
)[,"Estimate"]

25.5.3 Table of Next Season Predictions

Code
player_stats_seasonal_offense_subsetCC |> 
  filter(season == max(player_stats_seasonal_offense_subsetCC$season), position == "QB") |>
  arrange(-fantasyPoints_bayesian) |> 
  select(player_display_name, fantasyPoints_bayesian)
Code
player_stats_seasonal_offense_subsetCC |> 
  filter(season == max(player_stats_seasonal_offense_subsetCC$season), position == "RB") |>
  arrange(-fantasyPoints_bayesian) |> 
  select(player_display_name, fantasyPoints_bayesian)
Code
player_stats_seasonal_offense_subsetCC |> 
  filter(season == max(player_stats_seasonal_offense_subsetCC$season), position == "WR") |>
  arrange(-fantasyPoints_bayesian) |> 
  select(player_display_name, fantasyPoints_bayesian)
Code
player_stats_seasonal_offense_subsetCC |> 
  filter(season == max(player_stats_seasonal_offense_subsetCC$season), position == "TE") |>
  arrange(-fantasyPoints_bayesian) |> 
  select(player_display_name, fantasyPoints_bayesian)

25.5.4 Plot of Individuals’ Model-Implied Predictions

25.5.4.1 Quarterbacks

Code
plot_individualFantasyPointsByAgeQB <- ggplot(
  data = player_stats_seasonal_offense_subsetCC |> filter(position == "QB"),
  mapping = aes(
    x = round(age, 2),
    y = round(fantasyPoints_bayesian, 2),
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    se = FALSE,
    linewidth = 0.5,
    color = "black") +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    size = 1,
    color = "transparent" # make points invisible but keep tooltips
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Quarterbacks"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsByAgeQB,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
Figure 25.13: Plot of Individuals’ Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Quarterbacks.

25.5.4.2 Running Backs

Code
plot_individualFantasyPointsByAgeRB <- ggplot(
  data = player_stats_seasonal_offense_subsetCC |> filter(position == "RB"),
  mapping = aes(
    x = age,
    y = fantasyPoints_bayesian,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    se = FALSE,
    linewidth = 0.5,
    color = "black") +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    size = 1,
    color = "transparent" # make points invisible but keep tooltips
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Running Backs"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsByAgeRB,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
Figure 25.14: Plot of Individuals’ Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Running Backs.

25.5.4.3 Wide Receivers

Code
plot_individualFantasyPointsByAgeWR <- ggplot(
  data = player_stats_seasonal_offense_subsetCC |> filter(position == "WR"),
  mapping = aes(
    x = age,
    y = fantasyPoints_bayesian,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    se = FALSE,
    linewidth = 0.5,
    color = "black") +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    size = 1,
    color = "transparent" # make points invisible but keep tooltips
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Wide Receivers"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsByAgeWR,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
Figure 25.15: Plot of Individuals’ Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Wide Receivers.

25.5.4.4 Tight Ends

Code
plot_individualFantasyPointsByAgeTE <- ggplot(
  data = player_stats_seasonal_offense_subsetCC |> filter(position == "TE"),
  mapping = aes(
    x = age,
    y = fantasyPoints_bayesian,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    se = FALSE,
    linewidth = 0.5,
    color = "black") +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_bayesian,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    size = 1,
    color = "transparent" # make points invisible but keep tooltips
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Tight Ends"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsByAgeTE,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
Figure 25.16: Plot of Individuals’ Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Tight Ends.

25.6 Conclusion

In sum, time series analysis seeks to evaluate change over time to predict future values. There are many different types of time series analyses. We demonstrated use of autoregressive integrated moving average (ARIMA) models to predict future fantasy points. ARIMA models aim to describe how earlier levels of a variable are correlated with later levels of the same variable. ARIMA models perform best when there is a clear pattern where later values are influenced by earlier values. ARIMA models assume that the data are stationary (i.e., there are no long-term trends), are non-seasonal (i.e., there is no consistency of the timing of the peaks or troughs in the line), and that earlier values influence later values. This may not strongly be the case in fantasy football, so ARIMA models may not be particularly useful in forecasting fantasy football performance. We also used Bayesian mixed models to generate forecasts of future performance and plots of individuals model-implied performance by age and position.

25.7 Session Info

Code
sessionInfo()
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
 [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
 [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
[10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   

time zone: UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lubridate_1.9.5    forcats_1.0.1      stringr_1.6.0      dplyr_1.2.1       
 [5] purrr_1.2.2        readr_2.2.0        tidyr_1.3.2        tibble_3.3.1      
 [9] tidyverse_2.0.0    plotly_4.12.1      ggplot2_4.0.3      rstan_2.32.7      
[13] StanHeaders_2.39.1 brms_2.23.0        Rcpp_1.1.2         forecast_9.0.2    
[17] xts_0.14.2         zoo_1.9-0          petersenlab_1.2.3 

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3   tensorA_0.36.2.1     rstudioapi_0.19.0   
  [4] jsonlite_2.0.0       magrittr_2.0.5       TH.data_1.1-5       
  [7] estimability_2.0.0   farver_2.1.2         nloptr_2.2.1        
 [10] rmarkdown_2.32       vctrs_0.7.3          minqa_1.2.8         
 [13] base64enc_0.1-6      htmltools_0.5.9      distributional_0.8.1
 [16] curl_8.0.0           Formula_1.2-6        htmlwidgets_1.6.4   
 [19] plyr_1.8.9           sandwich_3.1-3       emmeans_2.0.4       
 [22] lifecycle_1.0.5      pkgconfig_2.0.3      Matrix_1.7-5        
 [25] R6_2.6.1             fastmap_1.2.0        rbibutils_2.4.1     
 [28] digest_0.6.39        colorspace_2.1-3     ps_1.9.3            
 [31] crosstalk_1.2.2      Hmisc_5.3-0          labeling_0.4.3      
 [34] timechange_0.4.0     httr_1.4.9           abind_1.4-8         
 [37] mgcv_1.9-4           compiler_4.6.1       withr_3.0.3         
 [40] htmlTable_2.5.0      S7_0.2.2             backports_1.5.1     
 [43] inline_0.3.21        DBI_1.3.0            psych_2.6.5         
 [46] QuickJSR_1.11.0      pkgbuild_1.4.8       MASS_7.3-65         
 [49] loo_2.10.1           tools_4.6.1          pbivnorm_0.6.0      
 [52] foreign_0.8-91       otel_0.2.0           nnet_7.3-20         
 [55] glue_1.8.1           quadprog_1.5-8       nlme_3.1-169        
 [58] grid_4.6.1           cmdstanr_0.9.0.9002  checkmate_2.3.4     
 [61] cluster_2.1.8.2      reshape2_1.4.5       generics_0.1.4      
 [64] gtable_0.3.6         tzdb_0.5.0           data.table_1.18.6.1 
 [67] hms_1.1.4            pillar_1.11.1        posterior_1.7.0     
 [70] mitools_2.7          splines_4.6.1        lattice_0.22-9      
 [73] survival_3.8-6       tidyselect_1.2.1     mix_1.0-13          
 [76] knitr_1.52           reformulas_0.4.4     gridExtra_2.3.1     
 [79] V8_8.2.0             urca_1.3-4           stats4_4.6.1        
 [82] xfun_0.60            bridgesampling_1.2-1 timeDate_4052.112   
 [85] matrixStats_1.5.0    stringi_1.8.9        yaml_2.3.12         
 [88] boot_1.3-32          evaluate_1.0.5       codetools_0.2-20    
 [91] cli_3.6.6            RcppParallel_6.2.1   rpart_4.1.27        
 [94] xtable_1.8-8         Rdpack_2.6.6         processx_3.9.0      
 [97] lavaan_0.7-2         coda_0.19-4.1        parallel_4.6.1      
[100] rstantools_2.7.1     fracdiff_1.5-4       bayesplot_1.16.0    
[103] Brobdingnag_1.2-9    lme4_2.0-6           viridisLite_0.4.3   
[106] mvtnorm_1.4-2        scales_1.4.0         rlang_1.3.0         
[109] multcomp_1.4-32      mnormt_2.1.2

Feedback

Please consider providing feedback about this textbook, so that I can make it as helpful as possible. You can provide feedback at the following link: https://forms.gle/LsnVKwqmS1VuxWD18

Email Notification

The online version of this book will remain open access. If you want to know when the print version of the book is for sale, enter your email below so I can let you know.