time-series-analysis.qmd

# Time Series Analysis {#sec-timeSeries}

## Getting Started {#sec-timeSeriesGettingStarted}

### Load Packages {#sec-timeSeriesLoadPackages}

```{r}
library("petersenlab")
library("xts")
library("zoo")
library("forecast")
library("brms")
library("rstan")
library("plotly")
library("tidyverse")
```

### Load Data {#sec-timeSeriesLoadData}

```{r}
#| eval: false
#| include: false

load(file = file.path(path, "/OneDrive - University of Iowa/Teaching/Courses/Fantasy Football/Data/player_stats_weekly.RData", fsep = ""))
load(file = file.path(path, "/OneDrive - University of Iowa/Teaching/Courses/Fantasy Football/Data/player_stats_seasonal.RData", fsep = ""))
load(file = file.path(path, "/OneDrive - University of Iowa/Teaching/Courses/Fantasy Football/Data/bayesianMixedModelFit.RData", fsep = ""))
```

```{r}
load(file = "./data/player_stats_weekly.RData")
load(file = "./data/player_stats_seasonal.RData")
load(url("https://osf.io/download/q6rjf/"))
```

## Overview of Time Series Analysis {#sec-timeSeriesOverview}

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 many different types of time series analyses.
For simplicity, in this chapter, we use [autoregressive integrated moving average](#sec-timeSeriesARIMA) (ARIMA) models to demonstrate one approach to time series analysis.
We also leverage [Bayesian mixed models](#sec-fantasyPointsByAgeExperienceModelsBayesian) to generate forecasts of future performance and plots of individuals model-implied performance by age and position.

## Autoregressive Integrated Moving Average (ARIMA) Models {#sec-timeSeriesARIMA}

@Hyndman2021 provide a nice overview of ARIMA models.
As noted by @Hyndman2021, 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 [@Hyndman2021].
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> (archived at <https://perma.cc/U5P6-2VWC>).

### Create the Time Series Objects {#sec-createTimeSeriesObjects}

```{r}
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
ts_peytonManning

ts_combined <- merge(
  ts_tomBrady,
  ts_peytonManning
)

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

### Plot the Time Series {#sec-plotTimeSeries}

```{r}
#| label: fig-timeSeriesHistoricalTomBrady
#| fig-cap: "Tom Brady's Historical Fantasy Points by Game."
#| fig-alt: "Tom Brady's Historical Fantasy Points by Game."

plot(
  ts_tomBrady,
  main = "Tom Brady's Fantasy Points by Game")
```

```{r}
#| label: fig-timeSeriesHistoricalBradyManning
#| fig-cap: "Historical Fantasy Points by Game for Tom Brady and Peyone Manning."
#| fig-alt: "Historical Fantasy Points by Game for Tom Brady and Peyone Manning."

plot(
  ts_combined,
  legend,
  legend.loc = "topright",
  main = "Fantasy Points by Game")
```

### Rolling Mean/Median {#sec-rollingMeanMedian}

```{r}
zoo::rollmean(
  x = ts_tomBrady,
  k = 5)

zoo::rollmedian(
  x = ts_tomBrady,
  k = 5)
```

### Autocorrelation {#sec-autoCorrelation}

The autocorrelation function (ACF) plot depicts the autocorrelation of scores as a function of the length of the lag.
Significant autocorrelation is detected when the autocorrelation exceeds the dashed blue lines, as is depicted in @fig-timeSeriesACFtomBrady.

```{r}
#| label: fig-timeSeriesACFtomBrady
#| fig-cap: "Autocorrelation Function (ACF) Plot of Tom Brady's Historical Fantasy Points by Game."
#| fig-alt: "Autocorrelation Function (ACF) Plot of Tom Brady's Historical Fantasy Points by Game."

acf(ts_tomBrady)
```

```{r}
Box.test(ts_tomBrady)
```

### Fit an Autoregressive Integrated Moving Average Model {#sec-arimaModel}

```{r}
forecast::auto.arima(ts_tomBrady)
forecast::auto.arima(ts_peytonManning)
```

```{r}
#| label: fig-arimaTomBrady
#| fig-cap: "Model Summary of Autoregressive Integrated Moving Average Model fit to Tom Brady's Historical Performance by Game."
#| fig-alt: "Model Summary of Autoregressive Integrated Moving Average Model fit to Tom Brady's Historical Performance by Game."

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

summary(arima_tomBrady)
confint(arima_tomBrady)
forecast::checkresiduals(arima_tomBrady)
```

```{r}
#| label: fig-arimaModifiedTomBrady
#| fig-cap: "Model Summary of modified Autoregressive Integrated Moving Average Model fit to Tom Brady's Historical Performance by Game."
#| fig-alt: "Model Summary of modified Autoregressive Integrated Moving Average Model fit to Tom Brady's Historical Performance by Game."

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)
confint(arima_tomBrady_removeNonSigTerms)
forecast::checkresiduals(arima_tomBrady_removeNonSigTerms)
```

```{r}
#| label: fig-arimaPeytonManning
#| fig-cap: "Model Summary of Autoregressive Integrated Moving Average Model fit to Peyton Manning's Historical Performance by Game."
#| fig-alt: "Model Summary of Autoregressive Integrated Moving Average Model fit to Peyton Manning's Historical Performance by Game."

arima_peytonManning <- arima(
  ts_peytonManning,
  order = c(1, 0, 1))

summary(arima_peytonManning)
confint(arima_peytonManning)
forecast::checkresiduals(arima_peytonManning)
```

### Generate the Model Forecasts {#sec-arimaForecasts}

```{r}
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
forecast_peytonManning
```

### Plot the Model Forecasts {#sec-plotARIMAforecast}

```{r}
#| label: fig-timeSeriesProjectionTomBrady
#| fig-cap: "Tom Brady's Historical and Projected Fantasy Points by Game."
#| fig-alt: "Tom Brady's Historical and Projected Fantasy Points by Game."

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()
```

```{r}
#| label: fig-timeSeriesProjectionPeytonManning
#| fig-cap: "Peyton Manning's Historical and Projected Fantasy Points by Game."
#| fig-alt: "Peyton Manning's Historical and Projected Fantasy Points by Game."

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()
```

## Bayesian Mixed Models {#sec-timeSeriesBayesian}

The Bayesian longitudinal mixed models were estimated in @sec-fantasyPointsByAgeExperienceModelsBayesian.

### Prepare New Data Object {#sec-timeSeriesBayesianNewData}

```{r}
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)
```

```{r}
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) 
  )
```

### Generate Predictions {#sec-timeSeriesBayesianGeneratePredictions}

```{r}
player_stats_seasonal_offense_subsetCC$fantasyPoints_bayesian <- predict(
  bayesianMixedModelFit,
  newdata = player_stats_seasonal_offense_subsetCC
)[,"Estimate"]
```

### Table of Next Season Predictions {#sec-timeSeriesBayesianNextSeasonPredictions}

```{r}
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)

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)

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)

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)
```

### Plot of Individuals' Model-Implied Predictions {#sec-timeSeriesBayesianModelImpliedPredictions}

#### Quarterbacks {#sec-fantasyPointsIndividualByAgeTimeSeriesPlotsQB}

```{r}
#| label: fig-individualFantasyPointsBayesianByAgeQB
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Quarterbacks."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Quarterbacks."
#| code-fold: true

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()

ggplotly(
  plot_individualFantasyPointsByAgeQB,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
```

#### Running Backs {#sec-fantasyPointsIndividualByAgeTimeSeriesPlotsRB}

```{r}
#| label: fig-individualFantasyPointsBayesianByAgeRB
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Running Backs."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Running Backs."
#| code-fold: true

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()

ggplotly(
  plot_individualFantasyPointsByAgeRB,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
```

#### Wide Receivers {#sec-fantasyPointsIndividualByAgeTimeSeriesPlotsWR}

```{r}
#| label: fig-individualFantasyPointsBayesianByAgeWR
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Wide Receivers."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Wide Receivers."
#| code-fold: true

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()

ggplotly(
  plot_individualFantasyPointsByAgeWR,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
```

#### Tight Ends {#sec-fantasyPointsIndividualByAgeTimeSeriesPlotsTE}

```{r}
#| label: fig-individualFantasyPointsBayesianByAgeTE
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Tight Ends."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Bayesian Generalized Additive Model, for Tight Ends."
#| code-fold: true

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: Wide Receivers"
  ) +
  theme_classic()

ggplotly(
  plot_individualFantasyPointsByAgeTE,
  tooltip = c("age","fantasyPoints_bayesian","text","label")
)
```

## Conclusion {#sec-timeSeriesConclusion}

::: {.content-visible when-format="html"}

## Session Info {#sec-timeSeriesSessionInfo}

```{r}
sessionInfo()
```

:::