mixed-models.qmd

{{< include _chunk-timing.qmd >}}

# Mixed Models {#sec-mixedModels}

This chapter provides an overview of mixed models.

## Getting Started {#sec-mixedModelsGettingStarted}

### Load Packages {#sec-mixedModelsLoadPackages}

```{r}
library("lme4")
library("lmerTest")
library("effectsize")
library("MuMIn")
library("emmeans")
library("sjstats")
library("mgcv")
library("AICcmodavg")
library("bbmle")
library("rstan")
library("brms")
library("cmdstanr") # todo: install.packages("cmdstanr", repos = c("https://stan-dev.r-universe.dev", getOption("repos"))); cmdstanr::check_cmdstan_toolchain(); cmdstanr::install_cmdstan()
library("fitdistrplus")
library("performance")
library("parallelly")
library("broom.mixed")
library("tidybayes")
library("plotly")
library("viridis")
library("tidyverse")
```

### Specify Package Options {#sec-mixedModelsPackageOptions}

```{r}
emm_options(lmerTest.limit = 100000)
emm_options(pbkrtest.limit = 100000)
```

### Load Data {#sec-mixedModelsLoadData}

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

load(file = "./data/nfl_depthCharts.RData")
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 = ""))
load(file = file.path(path, "/OneDrive - University of Iowa/Teaching/Courses/Fantasy Football/Data/modelPredictions_bayesian.RData", fsep = ""))
```

```{r}
load(file = "./data/nfl_depthCharts.RData")
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 @sec-calculatePlayerAge.

## Overview of Mixed Models {#sec-mixedModelsOverview}

Mixed models are simply regression models that account for the nested or hierarchical structure of the data.
Mixed models can be be used to address the [assumption in multiple regression](#sec-assumptionsRegression) that the residuals are independent (i.e., uncorrelated).
When data come from the same group, that can lead to residuals being correlated, so accounting for the grouping-level structure can address the assumption.

The modeling framework goes by many terms, including mixed models, mixed-effects models, multilevel models, hierarchical linear models.
They are sometimes called multilevel models and hierarchical linear models, whose name emphasizes the hierarchical structure of the data because the data are nonindependent.
When observations (i.e., data points) are collected from multiple lower-level units (e.g., people) in an upper-level unit (e.g., married couple, family, classroom, school, neighborhood, team), the data from the lower-level units are considered "nested" within the upper-level unit.
In this context, *nested data* refers to multiple observations from the same upper-level unit.
For instance, longitudinal data are nested within the same participant.
Students are nested within classrooms, which are nested within schools.
Players are nested within teams.

When data are nested, the data from the lower-level unit are likely to be correlated, to some degree, because they come from the same upper-level unit.
For example, multiple players may come from the same team, and the players' performance on that team is likely interrelated because they share common experiences and influence one another.
Thus, data from multiple players on a given team are considered nested within that team.
Longitudinal data can also be considered nested data, in which time points are nested within the person (i.e., the same player provides an observation across multiple time points).
As we will discuss, it is important to account for levels of nesting when the observations are nonindependent.
Mixed models provide a framework for accounting for levels of nesting, in order to account for why data from the same upper-level unit are likely to be intercorrelated.

These models are also sometimes called mixed models or mixed-effects models because the models can include a mix of fixed and random effects.
Fixed effects are effects that are constant across individuals (i.e., upper-level units).
Random effects are effects that vary across individuals (i.e., upper-level units).
For instance, consider a longitudinal study of fantasy performance as a function of age.
If we have longitudinal data for multiple players, the time points are nested within players.
Examining the association between age as a fixed effect in relation to fantasy performance would examine the association between a player's age and their fantasy performance while holding the association between age and fantasy performance constant across all players.
That is, it would assume that all players show the same trajectory such as increase from ages 20 to 24 and then decrease.
Examining the association between age as a random effect in relation to fantasy performance would examine the association between a player's age and their fantasy performance while allowing the association between age and fantasy performance to vary across players.
That is, it would allow the possibility that some players improve with age, whereas other players decline in performance with age.

When including random effects of a variable (e.g., age) in a mixed model, it is also important to include fixed effects of that variable in the model.
This is because random effects have a mean of zero.
Fixed effects allow the mean to differ from zero.
Thus, inclusion of random effects without the corresponding fixed effect can lead to bias in estimation of the association between the [predictor variables](#sec-correlationalStudy) and the [outcome variable](#sec-correlationalStudy).

A visualization of mixed models, including the distinction between fixed and random intercepts and slopes, is provided by @Freeman2017: <http://mfviz.com/hierarchical-models> (archived at <https://perma.cc/H2GZ-P5RW>).

Mixed modeling is a more flexible modeling framework than [analysis of variance](#sec-anova) (ANOVA).
Unlike [ANOVA](#sec-anova), mixed modeling can accommodate both continuous and categorical [predictor variables](#sec-correlationalStudy) that vary within units [@Brauer2018].
[ANOVA](#sec-anova), cannot handle continuous [predictor variables](#sec-correlationalStudy) that vary within units [@Brauer2018].
In addition, [ANOVA](#sec-anova) yields biased inferences (i.e., a higher Type I error rate) when the same participants are exposed to multiple items or stimuli [i.e., responses by the same respondent tend to be correlated across items, and responses to the same item tend to be correlated across respondents; @Brauer2018].
Additionally, mixed models better handle missing data compared to [ANOVA](#sec-anova) [@Brauer2018].

### Ecological Fallacy {#sec-ecologicalFallacy}

As described in @sec-fallaciesEcological, the ecological fallacy is the error of drawing inferences about an individual from group-level data.
A type of ecological fallacy is [Simpson's paradox](#sec-simpsonsParadox).

### Simpson's Paradox {#sec-simpsonsParadox}

Simpson's paradox occurs when the association between the [predictor variable](#sec-correlationalStudy) and [outcome variable](#sec-correlationalStudy) for the subgroups differ from the association when the subgroups are combined.
Examples of Simpson's paradox are depicted in Figures [-@fig-simpsonParadox] and [-@fig-simpsonsParadox].
In the example below, there is a positive association between the [predictor variable](#sec-correlationalStudy) and [outcome variable](#sec-correlationalStudy) for each group.
However, when the groups are combined, there is a *negative* association between [predictor variable](#sec-correlationalStudy) and [outcome variable](#sec-correlationalStudy).
That is, the sign of an association can differ at different levels of analysis (e.g., group level versus person level).

::: {#fig-simpsonsParadox}
![](images/simpsonsParadox.png){fig-alt="An Example of Simpson's Paradox. In this example, there is a positive association between `x` and `y` within each group (i.e., red or blue), but there is a negative association between `x` and `y` when the groups are combined. (Figure retrieved from <https://en.wikipedia.org/wiki/File:Simpson%27s_paradox_continuous.svg>)"}

An Example of Simpson's Paradox. In this example, there is a positive association between `x` and `y` within each group (i.e., red or blue), but there is a negative association between `x` and `y` when the groups are combined. (Figure retrieved from <https://en.wikipedia.org/wiki/File:Simpson%27s_paradox_continuous.svg>)
:::

Consider that we observe a between-person association between a [predictor](#sec-correlationalStudy), for example, how much sports drink a player drinks and their performance, such that players who drink more sports drink before games tend to perform better during games.
However, based on this association, if we draw the inference that sports drink consumption leads to better performance, this could be a faulty inference.
It is possible that, at the *within*-person level, there is no association or even a negative association between sports drink consumption and performance.
For helping to approximate causality, a much stronger test than relying on the between-person association would be to examine the association within the individual.
That is, examining the association within the individual would examine: when a player consumes more sports drink, whether they perform better than when the same player consumes less sports drink.
We describe within-person analyses to approximate causal inferences in @sec-causalInferenceWithinSubject.

In short, we often need to account for the groups or "levels" within which people are nested.
When multiple observations are nested within the same upper-level unit, they are often correlated.
That is, some variance is attributable to the lower-level unit (e.g., the player) and some variance is attributable to the upper-level unit (e.g., the team).
It is important to evaluate how much variance is attributable to the upper-level unit.
A way of evaluating how much variance is attributable to an upper-level unit is with the intraclass correlation coefficient (ICC).
If substantial variance is attributable to the upper-level unit, it is important to account for the nested structure of the data.
Mixed models are a useful approach to account for data with a nested structure so you can avoid committing the ecological fallacy.

### Data Structure {#sec-mixedModelDataStructure}

For fitting mixed models, the data should be in [long form](#sec-longWideForm)—each player (or unit of interest) should have multiple rows based on the grouping structure.
The cluster variable is the variable or variables that represent the upper-level unit(s), within which the data are nested.
When estimating mixed models, it is important to know the grouping structure of the data because the data structure influences the model syntax.
Examples of data from nonindependent units include nested data, multiple membership data, and cross-classified data, as depicted in @fig-nestedDataStructures.

::: {#fig-nestedDataStructures}
![](images/nestedData.png){fig-alt="Various Data Structures Involving Data from Nonindependent Units, Including Nested data, Multiple Membership Data, and Cross-Classified Data."}

Various Data Structures Involving Data from Nonindependent Units, Including Nested data, Multiple Membership Data, and Cross-Classified Data.
:::

One example of nonindependent data is if your data are *nested*.
When data are collected from multiple lower-level units (e.g., people) in an upper-level unit (e.g., married couple, family, classroom, school, neighborhood), and each lower-level unit belongs to one and only one higher-level unit, the data from the lower-level units are considered nested within the upper-level unit.
For example, multiple participants in your sample may come from the same classroom.
Data from multiple people can be nested within the same family, classroom, school, etc.
Longitudinal data can also be considered nested data, in which time points are nested within the participant (i.e., the same participant provides an observation across multiple time points).
And if you have multiple informants of a child's behavior (e.g., parent-, teacher-, friend-, and self-report), the ratings could also be nested within the participant (i.e., there are ratings from multiple informants of the same child).

Another form of nonindependent data is when data involve *multiple membership*.
Data involve multiple membership when the lower-level units belong to more than one upper-level unit.
As an example of multiple membership, children may have more than one teacher, and therefore, in a modeling sense, children "belong" to more than one teacher cluster.

Another form of nonindependent data is when data are *cross-classified* (also called crossed).
Data are cross-classified when the lower-level units are classified by two or more upper-level units, but the upper-level units are not hierarchical or nested within one another.
For instance, children may be nested within the crossing of schools and neighborhoods.
That is, children are nested within schools, and children are also nested within neighborhoods; however, children attending the same school may not necessarily be from the same neighborhood, and children from the same neighborhood may not necessarily attend the same school.
The data would still have a two-level hierarchy, but the data would be nested in specific school-neighborhood combinations.
That is, children are nested within the cross-classification of schools and neighborhoods.
An example of cross-classified data in football is that players are nested within teams and within agents.
Agents are not nested within teams, and teams are not nested within agents.
That is, teams are nested within the cross-classification of teams and agents.

Applied to longitudinal fantasy football data, most simply, players' longitudinal performance data would be modeled as nested data with a two-level model: timepoints nested within players, with the player's ID as the cluster variable.
However, you could get more complicated.
For instance, you could estimate a three-level model of timepoints nested within players nested within teams.
This would help account for players' productivity benefitting or taking a hit from being on a given team, and it would account for players changing teams.
Alternatively, consistent with @Usami2018, you could estimate a cross-classified model of players' individual observations/measurement occasions belonging to two upper-level units: players and time points (e.g., week 1, week 2, 2025 season, etc.).
This would help account for systematic effects that are specific to a given timepoint (e.g., week or season).

## Handling of Missing Data {#sec-mixedModelMissingness}

Mixed models handle missing data better than [ANOVA](#sec-anova) and [multiple regression](#sec-multipleRegression) [@Brauer2018].
[ANOVA](#sec-anova) and [multiple regression](#sec-multipleRegression) use listwise deletion by default and, thus, exclude participants who have missing values on any of the variables in the model (i.e., [independent/dependent](#sec-experiment) or [predictor/outcome](#sec-correlationalStudy)).
By contrast, with the data in [long form](#sec-longWideForm), mixed models make use of all rows that have values on all of the [predictor variables](#sec-correlationalStudy).
That is, mixed models exclude only those rows that have a missing value on any of the [predictor variables](#sec-correlationalStudy).
Unlike [ANOVA](#sec-anova) and [multiple regression](#sec-multipleRegression), mixed models retain rows with missing values on the [outcome variables](#sec-correlationalStudy).

### Comparing Model Fit {#sec-mixedModelCompareModels}

When trying to determine what [predictors](#sec-correlationalStudy) to include and how complex of a model to fit, it can be helpful to compare a variety a models to determine how well they fit the data (i.e., how much variance they account for).
To compare non-nested models, you can use fit criteria such as Akaike information criterion (AIC), which penalizes model complexity.
If two models are nested, you can compare their model fit using a likelihood ratio test.
Two models are considered nested if one model includes all of the terms of the original model along with additional terms.
The simpler model is considered "nested within" the more complex model.
In the likelihood ratio test, a significant difference indicates that the more complex model is significantly better fitting than the simpler model.
If the likelihood ratio test is non-significant, it indicates that the more complex model is not significantly better fitting, and thus you should prefer the more parsimonious model.

## Fantasy Points Per Season by Position, Age, and Experience {#sec-fantasyPointsByAgeExperience}

In this chapter, we use mixed models to demonstrate how to model trajectories (growth curves) of players' fantasy points as a function of age.
When estimating growth curve models, time points are nested within individuals (i.e., players).
Mixed models provide a flexible framework for estimating growth curve models because they allow unequal time intervals between measurement occasions, time intervals that differ between individuals, and different numbers of observations for each individual [@Brauer2018].

One of the challenges when modeling players' fantasy points as a function of age is that fantasy points are a function of both ability and opportunity.
With age, the player's ability and opportunity may both decline, but it is difficult to disentangle the extent to which a player's decline is related to declines in ability versus opportunity.
As players get older, they may be supplanted by younger, more talented players.
With age, despite still having latent (unobserved) ability, players will get fewer and, eventually, no opportunity.
And we do not know the counterfactual—how many fantasy points they would have scored had they been given the full opportunity each season.
Thus, players may go from scoring 100+ points in a season to all of a sudden scoring way fewer points, with little in the way of intermediate steps.
Thus, this should inform how you interpret the resulting trajectories.
For instance, the steep declines can make it look like the model-implied fantasy points go well below zero at later ages, even though it is obviously not likely for players to obtain such a magnitude of negative fantasy points.
The negative model-implied values of fantasy points at later ages are likely an artifact of the decreasing opportunity that players tend to get as they age.
The crossover point, where the positive values becomes negative might indicate, for instance, that players at that age tend to have retired and thus have no opportunity.
Regardless, the form/shape/steepness of the decline is still potentially meaningful even if the precise negative number at a given age is not.

```{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}
#| eval: false
#| include: false

for(i in unique(nfl_depthCharts$season)){
  print(i)
  nfl_depthCharts %>%
    filter(season == i) %>% 
    select(week) %>% 
    arrange(week) %>% 
    unique() %>% 
    pull() %>% 
    print()
}
```

```{r}
seasons17week <- 2001:2020
seasons18week <- 2021:max(nfl_depthCharts$season, na.rm = TRUE)

endOfSeasonDepthCharts <- nfl_depthCharts %>% 
  filter((season %in% seasons17week & week == 18) | (season %in% seasons18week & week == 19)) # get end-of-season depth charts

qb1s <- endOfSeasonDepthCharts %>% 
  filter(position == "QB", depth_team == 1)

fb1s <- endOfSeasonDepthCharts %>% 
  filter(position == "FB", depth_team == 1)

k1s <- endOfSeasonDepthCharts %>% 
  filter(position == "K", depth_team == 1)

rb1s <- endOfSeasonDepthCharts %>% 
  filter(position == "RB", depth_team == 1)

wr1s <- endOfSeasonDepthCharts %>% 
  filter(position == "WR", depth_team == 1)

te1s <- endOfSeasonDepthCharts %>% 
  filter(position == "TE", depth_team == 1)

player_stats_seasonal_offense_subsetDepth <- player_stats_seasonal_offense_subset %>% 
  filter(player_id %in% c(
    qb1s$gsis_id,
    fb1s$gsis_id,
    k1s$gsis_id,
    rb1s$gsis_id,
    wr1s$gsis_id,
    te1s$gsis_id
    ))
```

Create a `newdata` object for generating the plots of model-implied fantasy points by age and position:

```{r}
pointsPerSeason_positionAge_newData <- expand.grid(
  positionFactor = factor(c("FB","QB","RB","TE","WR")), #,"K"
  age = seq(from = 20, to = 40, length.out = 10000)
)

pointsPerSeason_positionAge_newData$ageCentered20 <- pointsPerSeason_positionAge_newData$age - 20
pointsPerSeason_positionAge_newData$ageCentered20Quadratic <- pointsPerSeason_positionAge_newData$ageCentered20 ^ 2
pointsPerSeason_positionAge_newData$years_of_experience <- floor(pointsPerSeason_positionAge_newData$age - 22) # assuming that most players start at age 22 (i.e., rookie year) and thus have 1 year of experience at age 23
pointsPerSeason_positionAge_newData$years_of_experience[which(pointsPerSeason_positionAge_newData$years_of_experience < 0)] <- 0
```

Create an object with complete cases for generating the plots of individuals' model-implied fantasy points by age and 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))
```

### Scatterplots of Fantasy Points by Age and Position {#sec-fantasyPointsByAgeScatterplot}

[Scatterplots](#sec-scatterplot) are a helpful tool for quickly examining the association between two variables.
However, [scatterplots](#sec-scatterplot)—as well as [correlation](#sec-correlation) and [multiple regression](#sec-multipleRegression)—can hide meaningful associations that differ across units of analysis.

#### Quarterbacks {#sec-fantasyPointsByAgeScatterplotQB}

A [scatterplot](#sec-scatterplot) of Quarterbacks' fantasy points by age is in @fig-fantasyPointsByAgeScatterplotQB.

```{r}
#| label: fig-fantasyPointsByAgeScatterplotQB
#| fig-cap: "Scatterplot of Fantasy Points by Age for Quarterbacks."
#| fig-alt: "Scatterplot of Fantasy Points by Age for Quarterbacks."
#| code-fold: true

plot_scatterplotFantasyPointsByAgeQB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "QB") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasyPoints),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Quarterbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_scatterplotFantasyPointsByAgeQB,
  tooltip = c("age","fantasyPoints","text","label"))
```

Based on the [scatterplot](#sec-scatterplot) (and the bivariate association below), Quarterbacks' fantasy points appear to (slightly) increase with age.

```{r}
cor.test(
  formula = ~ age + fantasyPoints,
  data = player_stats_seasonal_offense_subset %>% filter(position == "QB")
)
```

#### Fullbacks {#sec-fantasyPointsByAgeScatterplotFB}

A [scatterplot](#sec-scatterplot) of Fullbacks' fantasy points by age is in @fig-fantasyPointsByAgeScatterplotFB.

```{r}
#| label: fig-fantasyPointsByAgeScatterplotFB
#| fig-cap: "Scatterplot of Fantasy Points by Age for Fullbacks."
#| fig-alt: "Scatterplot of Fantasy Points by Age for Fullbacks."
#| code-fold: true

plot_scatterplotFantasyPointsByAgeFB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "FB") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasyPoints),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Fullbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_scatterplotFantasyPointsByAgeFB,
  tooltip = c("age","fantasyPoints","text","label"))
```

Based on the [scatterplot](#sec-scatterplot), Fullbacks' fantasy points appear to be relatively stable across ages.

#### Running Backs {#sec-fantasyPointsByAgeScatterplotRB}

A [scatterplot](#sec-scatterplot) of Running Backs' fantasy points by age is in @fig-fantasyPointsByAgeScatterplotRB.

```{r}
#| label: fig-fantasyPointsByAgeScatterplotRB
#| fig-cap: "Scatterplot of Fantasy Points by Age for Running Backs."
#| fig-alt: "Scatterplot of Fantasy Points by Age for Running Backs."
#| code-fold: true

plot_scatterplotFantasyPointsByAgeRB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "RB") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasyPoints),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Running Backs"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_scatterplotFantasyPointsByAgeRB,
  tooltip = c("age","fantasyPoints","text","label"))
```

Based on the [scatterplot](#sec-scatterplot), Running Backs' fantasy points appear to be relatively stable across ages.

#### Wide Receivers {#sec-fantasyPointsByAgeScatterplotWR}

A [scatterplot](#sec-scatterplot) of Wide Receivers' fantasy points by age is in @fig-fantasyPointsByAgeScatterplotWR.

```{r}
#| label: fig-fantasyPointsByAgeScatterplotWR
#| fig-cap: "Scatterplot of Fantasy Points by Age for Wide Receivers."
#| fig-alt: "Scatterplot of Fantasy Points by Age for Wide Receivers."
#| code-fold: true

plot_scatterplotFantasyPointsByAgeWR <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "WR") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasyPoints),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Wide Receivers"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_scatterplotFantasyPointsByAgeWR,
  tooltip = c("age","fantasyPoints","text","label"))
```

Based on the [scatterplot](#sec-scatterplot), Wide Receivers' fantasy points appear to be relatively stable across ages.

#### Tight Ends {#sec-fantasyPointsByAgeScatterplotTE}

A [scatterplot](#sec-scatterplot) of Tight Ends' fantasy points by age is in @fig-fantasyPointsByAgeScatterplotTE.

```{r}
#| label: fig-fantasyPointsByAgeScatterplotTE
#| fig-cap: "Scatterplot of Fantasy Points by Age for Tight Ends."
#| fig-alt: "Scatterplot of Fantasy Points by Age for Tight Ends."
#| code-fold: true

plot_scatterplotFantasyPointsByAgeTE <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "TE") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_point(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  geom_smooth(
    mapping = aes(
    x = age,
    y = fantasyPoints),
    inherit.aes = FALSE
  ) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Tight Ends"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_scatterplotFantasyPointsByAgeTE,
  tooltip = c("age","fantasyPoints","text","label"))
```

Based on the [scatterplot](#sec-scatterplot) (and the bivariate association below), Tight Ends' fantasy points appear to increase with age.

```{r}
cor.test(
  formula = ~ age + fantasyPoints,
  data = player_stats_seasonal_offense_subset %>% filter(position == "TE")
)
```

### Plots of Raw Trajectories of Fantasy Points By Age and Player {#sec-fantasyPointsByAgeExperienceRawData}

[Scatterplots](#sec-scatterplot) can be helpful for quickly visualizing the association between two variables.
However, as mentioned earlier, [scatterplots](#sec-scatterplot) can hide the association between variables at different units of analysis.
For instance, consider that we are trying to predict how a player will perform based on their age.
We are interested not only in what the association is between age and fantasy points between players (i.e., a between-person association).
We are also interested in what the association is between age and fantasy points *within* a given player (and within each player; i.e., a within-individual association).
Arguably, the within-individual association between age and fantasy points is more relevant to the prediction of performance than the association between age and fantasy points between players.
Assuming that the between-player association between age and fantasy points is the same as the within-player association when it is not is an example of the [ecological fallacy](#sec-fallaciesEcological).

Below, we depict players' raw trajectories of fantasy points as a function of age.
These are known as spaghetti plots.
By examining the trajectory for each player, we can get a better understanding of how performance changes (within an individual) as a function of age.

#### Quarterbacks {#sec-fantasyPointsByAgeExperienceRawDataQB}

A plot of Quarterbacks' raw fantasy points data by age is in @fig-rawFantasyPointsByAgeQB.

```{r}
#| label: fig-rawFantasyPointsByAgeQB
#| fig-cap: "Plot of Raw Trajectories of Fantasy Points by Age for Quarterbacks."
#| fig-alt: "Plot of Raw Trajectories of Fantasy Points by Age for Quarterbacks."
#| code-fold: true

plot_rawFantasyPointsByAgeQB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "QB") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Quarterbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_rawFantasyPointsByAgeQB,
  tooltip = c("age","fantasyPoints","text","label"))
```

#### Fullbacks {#sec-fantasyPointsByAgeExperienceRawDataFB}

A plot of Fullbacks' raw fantasy points data by age is in @fig-rawFantasyPointsByAgeFB.

```{r}
#| label: fig-rawFantasyPointsByAgeFB
#| fig-cap: "Plot of Raw Trajectories of Fantasy Points by Age for Fullbacks."
#| fig-alt: "Plot of Raw Trajectories of Fantasy Points by Age for Fullbacks."
#| code-fold: true

plot_rawFantasyPointsByAgeFB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "FB") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Fullbacks"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_rawFantasyPointsByAgeFB,
  tooltip = c("age","fantasyPoints","text","label"))
```

#### Running Backs {#sec-fantasyPointsByAgeExperienceRawDataRB}

A plot of Running Backs' raw fantasy points data by age is in @fig-rawFantasyPointsByAgeRB.

```{r}
#| label: fig-rawFantasyPointsByAgeRB
#| fig-cap: "Plot of Raw Trajectories of Fantasy Points by Age for Running Backs."
#| fig-alt: "Plot of Raw Trajectories of Fantasy Points by Age for Running Backs."
#| code-fold: true

plot_rawFantasyPointsByAgeRB <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "RB") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Running Backs"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_rawFantasyPointsByAgeRB,
  tooltip = c("age","fantasyPoints","text","label"))
```

#### Wide Receivers {#sec-fantasyPointsByAgeExperienceRawDataWR}

A plot of Wide Receivers' raw fantasy points data by age is in @fig-rawFantasyPointsByAgeWR.

```{r}
#| label: fig-rawFantasyPointsByAgeWR
#| fig-cap: "Plot of Raw Trajectories of Fantasy Points by Age for Wide Receivers."
#| fig-alt: "Plot of Raw Trajectories of Fantasy Points by Age for Wide Receivers."
#| code-fold: true

plot_rawFantasyPointsByAgeWR <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "WR") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Wide Receivers"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_rawFantasyPointsByAgeWR,
  tooltip = c("age","fantasyPoints","text","label"))
```

#### Tight Ends {#sec-fantasyPointsByAgeExperienceRawDataTE}

A plot of Tight Ends' raw fantasy points data by age is in @fig-rawFantasyPointsByAgeTE.

```{r}
#| label: fig-rawFantasyPointsByAgeTE
#| fig-cap: "Plot of Raw Trajectories of Fantasy Points by Age for Tight Ends."
#| fig-alt: "Plot of Raw Trajectories of Fantasy Points by Age for Tight Ends."
#| code-fold: true

plot_rawFantasyPointsByAgeTE <- ggplot(
  data = player_stats_seasonal_offense_subset %>% 
    filter(position == "TE") %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints = round(fantasyPoints, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints,
    color = player_id)) +
  geom_line(
    aes(
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
  )) +
  scale_color_viridis(discrete = TRUE) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age: Tight Ends"
  ) +
  theme_classic() +
  theme(legend.position = "none")

plotly::ggplotly(
  plot_rawFantasyPointsByAgeTE,
  tooltip = c("age","fantasyPoints","text","label"))
```

### Linear Regression Models {#sec-fantasyPointsByAgeExperienceModelsRegression}

#### Null Model {#sec-fantasyPointsByAgeExperienceModelsRegressionNull}

We can estimate model fit using the `MuMIn` package [@R-MuMIn]:

```{r}
pointsPerSeason_nullModel <- lm(
  fantasyPoints ~ 1,
  data = player_stats_seasonal_offense_subset,
  na.action = "na.exclude"
)

summary(pointsPerSeason_nullModel)
print(effectsize::standardize_parameters(pointsPerSeason_nullModel, method = "refit"), digits = 2)

summary(pointsPerSeason_nullModel)$r.squared
AIC(pointsPerSeason_nullModel)
MuMIn::AICc(pointsPerSeason_nullModel)
```

A plot of the model-implied trajectories of fantasy points by age from the null model is in @fig-fantasyPointsNull.

```{r}
#| label: fig-fantasyPointsNull
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Null Model."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Null Model."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_nullModel <- predict(
  object = pointsPerSeason_nullModel,
  newdata = pointsPerSeason_positionAge_newData
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_nullModel
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Null Model"
  ) +
  theme_classic()
```

#### Linear Model {#sec-fantasyPointsByAgeExperienceModelsRegressionLinear}

```{r}
pointsPerSeason_linearRegression <- lm(
  fantasyPoints ~ positionFactor + ageCentered20 + positionFactor:ageCentered20,
  data = player_stats_seasonal_offense_subset,
  na.action = "na.exclude"
)

summary(pointsPerSeason_linearRegression)
print(effectsize::standardize_parameters(pointsPerSeason_linearRegression, method = "refit"), digits = 2)

summary(pointsPerSeason_linearRegression)$r.squared
AIC(pointsPerSeason_linearRegression)
MuMIn::AICc(pointsPerSeason_linearRegression)
```

A plot of the model-implied trajectories of fantasy points by age from the linear regression model is in @fig-fantasyPointsLinearRegression.

```{r}
#| label: fig-fantasyPointsLinearRegression
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Linear Regression Model."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Linear Regression Model."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_linearRegression <- predict(
  object = pointsPerSeason_linearRegression,
  newdata = pointsPerSeason_positionAge_newData
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_linearRegression,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Linear Regression Model",
    color = "Position"
  ) +
  theme_classic()
```

#### Quadratic Model {#sec-fantasyPointsByAgeExperienceModelsRegressionQuadratic}

```{r}
pointsPerSeason_quadraticRegression <- lm(
  fantasyPoints ~ positionFactor + ageCentered20 + ageCentered20Quadratic + positionFactor:ageCentered20 + positionFactor:ageCentered20Quadratic,
  data = player_stats_seasonal_offense_subset,
  na.action = "na.exclude"
)

summary(pointsPerSeason_quadraticRegression)
print(effectsize::standardize_parameters(pointsPerSeason_linearRegression, method = "refit"), digits = 2)

summary(pointsPerSeason_quadraticRegression)$r.squared
AIC(pointsPerSeason_quadraticRegression)
MuMIn::AICc(pointsPerSeason_quadraticRegression)
```

A plot of the model-implied trajectories of fantasy points by age from the regression model with a quadratic term for age is in @fig-fantasyPointsQuadraticRegression.

```{r}
#| label: fig-fantasyPointsQuadraticRegression
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Quadratic Regression Model."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Quadratic Regression Model."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_quadraticRegression <- predict(
  object = pointsPerSeason_quadraticRegression,
  newdata = pointsPerSeason_positionAge_newData
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_quadraticRegression,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Quadratic Regression Model",
    color = "Position"
  ) +
  theme_classic()
```

#### Compare Models {#sec-fantasyPointsByAgeExperienceModelsRegressionCompare}

We compare nested models using an $F$-test with the `stats::anova()` function.
We compare non-nested models using the `bbmle` [@R-bbmle] (`bbmle::AICtab()`) and `MuMIn` [@R-MuMIn] (`MuMIn::AICc()`) packages.

```{r}
anova(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression
)

AIC(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression
  )

lmModels <- list(
  "nullModel" = pointsPerSeason_nullModel,
  "linearRegression" = pointsPerSeason_linearRegression,
  "quadraticRegression" = pointsPerSeason_quadraticRegression
)

bbmle::AICtab(lmModels)

MuMIn::AICc(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression
)

summary(pointsPerSeason_nullModel)$r.squared
summary(pointsPerSeason_linearRegression)$r.squared
summary(pointsPerSeason_quadraticRegression)$r.squared

deviance(pointsPerSeason_nullModel)
deviance(pointsPerSeason_linearRegression)
deviance(pointsPerSeason_quadraticRegression)

logLik(pointsPerSeason_nullModel)
logLik(pointsPerSeason_linearRegression)
logLik(pointsPerSeason_quadraticRegression)
```

### Mixed Models {#sec-fantasyPointsByAgeExperienceModels}

By accounting for which player each observation comes from using mixed models, we can examine the association between age and fantasy points in a more meaningful way, without violating the [assumption in multiple regression](#sec-assumptionsRegression) that the observations are independent (i.e., that the residuals are uncorrelated).

#### Random Intercepts Model {#sec-fantasyPointsByAgeExperienceModelsIntercepts}

We estimate the multilevel models using the `lmerTest::lmer()` function of the `lmerTest` package [@R-lmerTest], which is an extension of the `lme4` package [@R-lme4; @Bates2015_packages]:

```{r}
pointsPerSeason_randomIntercepts <- lmerTest::lmer(
  fantasyPoints ~ 1 + (1 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_randomIntercepts)
print(effectsize::standardize_parameters(pointsPerSeason_randomIntercepts, method = "refit"), digits = 2)

performance::r2(pointsPerSeason_randomIntercepts)
performance::icc(pointsPerSeason_randomIntercepts)
AIC(pointsPerSeason_randomIntercepts)
AICc(pointsPerSeason_randomIntercepts)
```

A plot of the model-implied trajectories of fantasy points by age from the mixed model with random intercepts is in @fig-fantasyPointsRandomIntercepts.

```{r}
#| label: fig-fantasyPointsRandomIntercepts
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Random Intercepts Mixed Model."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Random Intercepts Mixed Model."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_randomIntercepts <- predict(
  object = pointsPerSeason_randomIntercepts,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_randomIntercepts
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age",
    subtitle = "Random Intercepts Model"
  ) +
  theme_classic()
```

A plot of individuals' model-implied trajectories of fantasy points by age from the mixed model with random intercepts is in @fig-individualFantasyPointsRandomIntercepts.

```{r}
#| label: fig-individualFantasyPointsRandomIntercepts
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Mixed Model with Random Intercepts. Overlaid with the Model-Implied Trajectory."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Mixed Model with Random Intercepts. Overlaid with the Model-Implied Trajectory."
#| code-fold: true

player_stats_seasonal_offense_subsetCC$fantasyPoints_randomIntercepts <- predict(
  object = pointsPerSeason_randomIntercepts,
  newdata = player_stats_seasonal_offense_subsetCC
)

plot_individualFantasyPointsRandomIntercepts <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints_randomIntercepts = round(fantasyPoints_randomIntercepts, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints_randomIntercepts,
    group = player_id)) +
  geom_line(
    aes(
      x = age,
      y = fantasyPoints_randomIntercepts,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    linewidth = 0.5,
    color = "black") +
  geom_line(
    mapping = aes(
      x = age,
      y = fantasyPoints_randomIntercepts
    ),
    data = pointsPerSeason_positionAge_newData %>% 
      mutate(
        age = round(age, 2),
        fantasyPoints = round(fantasyPoints_randomIntercepts, 2)
        ),
    inherit.aes = FALSE,
    se = TRUE,
    color = "#3366FF",
    linewidth = 2
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position: Random Intercepts Model",
    #color = "Position"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsRandomIntercepts,
  tooltip = c("age","fantasyPoints_randomIntercepts","text","label")
)
```

#### Random Intercepts Model with Position as Fixed-Effect Predictor {#sec-fantasyPointsByAgeExperienceModelsPosition}

To obtain standardized regression coefficients, we use the `effectsize::standardize_parameters()` function of the `effectsize` package [@R-effectsize; @BenShachar2020_packages].
We use the `emmeans::emmeans()` function of the `emmeans` package [@R-emmeans] to determine the model-implied number of points per season by position.

```{r}
pointsPerSeason_position <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + (1 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_position)
print(effectsize::standardize_parameters(pointsPerSeason_position, method = "refit"), digits = 2)

performance::r2(pointsPerSeason_position)
emmeans::emmeans(pointsPerSeason_position, "positionFactor")
performance::icc(pointsPerSeason_position)
AIC(pointsPerSeason_position)
```

A plot of the model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts and a fixed effect of position is in @fig-fantasyPointsPosition.

```{r}
#| label: fig-fantasyPointsPosition
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Random Intercepts Mixed Model With Position as a Fixed-Effect Predictor."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age in Random Intercepts Mixed Model With Position as a Fixed-Effect Predictor."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_position <- predict(
  object = pointsPerSeason_position,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_position,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Random Intercepts Model With Position as Fixed-Effect Predictor",
    color = "Position"
  ) +
  theme_classic()
```

A plot of individuals' model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts and a fixed effect of position is in @fig-individualFantasyPointsPosition.

```{r}
#| label: fig-individualFantasyPointsPosition
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Mixed Model With Random Intercepts and a Fixed-Effect of Position. Overlaid with the Model-Implied Trajectory by Position."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Mixed Model With Random Intercepts and a Fixed-Effect of Position. Overlaid with the Model-Implied Trajectory by Position."
#| code-fold: true

player_stats_seasonal_offense_subsetCC$fantasyPoints_position <- predict(
  object = pointsPerSeason_position,
  newdata = player_stats_seasonal_offense_subsetCC
)

plot_individualFantasyPointsPosition <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints_position = round(fantasyPoints_position, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints_position,
    color = positionFactor,
    group = player_id)) +
  geom_line(
    aes(
      x = age,
      y = fantasyPoints_position,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    linewidth = 0.5) +
  geom_line(
    mapping = aes(
      x = age,
      y = fantasyPoints_position,
      color = positionFactor
    ),
    data = pointsPerSeason_positionAge_newData %>% 
      mutate(
        age = round(age, 2),
        fantasyPoints_position = round(fantasyPoints_position, 2)
        ),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position:\nRandom Intercepts Model With Position As Predictor",
    color = "Position"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsPosition,
  tooltip = c("age","fantasyPoints_position","text","label")
)
```

#### Identify the Best-Fitting Functional Form of Age {#sec-fantasyPointsByAgeExperienceModelsAge}

##### Linear Models {#sec-fantasyPointsByAgeExperienceModelsAgeLinear}

###### Random Intercepts, Fixed Linear Slopes {#sec-fantasyPointsByAgeExperienceModelsAgeLinearFixed}

```{r}
pointsPerSeason_positionAgeFixedLinearSlopes <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + ageCentered20 + (1 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_positionAgeFixedLinearSlopes)
print(effectsize::standardize_parameters(pointsPerSeason_positionAgeFixedLinearSlopes, method = "refit"), digits = 2)

performance::r2(pointsPerSeason_positionAgeFixedLinearSlopes)
emmeans::emmeans(pointsPerSeason_positionAgeFixedLinearSlopes, "positionFactor")
emmeans::emmeans(pointsPerSeason_positionAgeFixedLinearSlopes, "ageCentered20")
performance::icc(pointsPerSeason_positionAgeFixedLinearSlopes)
AIC(pointsPerSeason_positionAgeFixedLinearSlopes)
```

A plot of the model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts and fixed linear slopes is in @fig-fantasyPointsFixedLinearSlopes.

```{r}
#| label: fig-fantasyPointsFixedLinearSlopes
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts and Fixed Linear Slopes."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts and Fixed Linear Slopes."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_fixedLinearSlopes <- predict(
  object = pointsPerSeason_positionAgeFixedLinearSlopes,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_fixedLinearSlopes,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Mixed Model with Random Intercepts and Fixed Linear Slopes",
    color = "Position"
  ) +
  theme_classic()
```

A plot of individuals model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts and fixed linear slopes is in @fig-individualFantasyPointsFixedLinearSlopes.

```{r}
#| label: fig-individualFantasyPointsFixedLinearSlopes
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts and Fixed Slopes. Overlaid with the Model-Implied Trajectory by Position."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts and Fixed Slopes. Overlaid with the Model-Implied Trajectory by Position."
#| code-fold: true

player_stats_seasonal_offense_subsetCC$fantasyPoints_fixedLinearSlopes <- predict(
  object = pointsPerSeason_positionAgeFixedLinearSlopes,
  newdata = player_stats_seasonal_offense_subsetCC
)

plot_individualFantasyPointsFixedLinearSlopes <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints_fixedLinearSlopes = round(fantasyPoints_fixedLinearSlopes, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints_fixedLinearSlopes,
    color = positionFactor,
    group = player_id)) +
  geom_line(
    aes(
      x = age,
      y = fantasyPoints_fixedLinearSlopes,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    linewidth = 0.5) +
  geom_line(
    mapping = aes(
      x = age,
      y = fantasyPoints_fixedLinearSlopes,
      color = positionFactor
    ),
    data = pointsPerSeason_positionAge_newData,
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position:\nModel With Random Intercepts and Fixed Slopes",
    color = "Position"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsFixedLinearSlopes,
  tooltip = c("age","fantasyPoints_fixedLinearSlopes","text","label")
)
```

###### Random Intercepts, Random Linear Slopes {#sec-fantasyPointsByAgeExperienceModelsAgeLinearRandom}

```{r}
pointsPerSeason_positionAgeRandomLinearSlopes <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + ageCentered20 + (1 + ageCentered20 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_positionAgeRandomLinearSlopes)
print(effectsize::standardize_parameters(pointsPerSeason_positionAgeRandomLinearSlopes, method = "refit"), digits = 2)

performance::r2(pointsPerSeason_positionAgeRandomLinearSlopes)
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearSlopes, "positionFactor")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearSlopes, "ageCentered20")
performance::icc(pointsPerSeason_positionAgeRandomLinearSlopes)
AIC(pointsPerSeason_positionAgeRandomLinearSlopes)
```

A plot of the model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts and random linear slopes is in @fig-fantasyPointsRandomLinearSlopes.

```{r}
#| label: fig-fantasyPointsRandomLinearSlopes
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts and Random Linear Slopes."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts and Random Linear Slopes."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_randomLinearSlopes <- predict(
  object = pointsPerSeason_positionAgeRandomLinearSlopes,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_randomLinearSlopes,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Mixed Model with Random Intercepts and Random Linear Slopes",
    color = "Position"
  ) +
  theme_classic()
```

A plot of individuals' model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts and random linear slopes is in @fig-fantasyPointsRandomLinearSlopes.

```{r}
#| label: fig-individualFantasyPointsRandomLinearSlopes
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts and Random Linear Slopes. Overlaid with the Model-Implied Trajectory by Position."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts and Random Linear Slopes. Overlaid with the Model-Implied Trajectory by Position."
#| code-fold: true

player_stats_seasonal_offense_subsetCC$fantasyPoints_randomLinearSlopes <- predict(
  object = pointsPerSeason_positionAgeRandomLinearSlopes,
  newdata = player_stats_seasonal_offense_subsetCC
)

plot_individualFantasyPointsRandomLinearSlopes <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints_randomLinearSlopes = round(fantasyPoints_randomLinearSlopes, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints_randomLinearSlopes,
    color = positionFactor,
    group = player_id)) +
  geom_line(
    aes(
      x = age,
      y = fantasyPoints_randomLinearSlopes,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    linewidth = 0.5) +
  geom_line(
    mapping = aes(
      x = age,
      y = fantasyPoints_randomLinearSlopes,
      color = positionFactor
    ),
    data = pointsPerSeason_positionAge_newData %>% 
      mutate(
        age = round(age, 2),
        fantasyPoints_randomLinearSlopes = round(fantasyPoints_randomLinearSlopes, 2)
        ),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position:\nModel With Random Intercepts and Random Linear Slopes",
    color = "Position"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsRandomLinearSlopes,
  tooltip = c("age","fantasyPoints_randomLinearSlopes","text","label")
)
```

##### Quadratic Models {#sec-fantasyPointsByAgeExperienceModelsAgeQuadratic}

###### Random Intercepts, Random Linear Slopes, Fixed Quadratic Slopes {#sec-fantasyPointsByAgeExperienceModelsAgeQuadraticFixed}

```{r}
pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + ageCentered20 + ageCentered20Quadratic + (1 + ageCentered20 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
print(effectsize::standardize_parameters(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes, method = "refit"), digits = 2)

performance::r2(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes, "positionFactor")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes, "ageCentered20")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes, "ageCentered20Quadratic")
performance::icc(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
AIC(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
```

A plot of the model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts, random linear slopes, and fixed quadratic slopes is in @fig-fantasyPointsRandomLinearFixedQuadracticSlopes.

```{r}
#| label: fig-fantasyPointsRandomLinearFixedQuadracticSlopes
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_randomLinearFixedQuadraticSlopes <- predict(
  object = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_randomLinearFixedQuadraticSlopes,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Mixed Model with Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes",
    color = "Position"
  ) +
  theme_classic()
```

A plot of individuals' model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts, random linear slopes, and fixed quadratic slopes is in @fig-individualFantasyPointsRandomLinearFixedQuadracticSlopes.

```{r}
#| label: fig-individualFantasyPointsRandomLinearFixedQuadracticSlopes
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes. Overlaid with the Model-Implied Trajectory by Position."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes. Overlaid with the Model-Implied Trajectory by Position."
#| code-fold: true

player_stats_seasonal_offense_subsetCC$fantasyPoints_randomLinearFixedQuadraticSlopes <- predict(
  object = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  newdata = player_stats_seasonal_offense_subsetCC
)

plot_individualFantasyPointsRandomLinearFixedQuadracticSlopes <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints_randomLinearFixedQuadraticSlopes = round(fantasyPoints_randomLinearFixedQuadraticSlopes, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints_randomLinearFixedQuadraticSlopes,
    color = positionFactor,
    group = player_id)) +
  geom_line(
    aes(
      x = age,
      y = fantasyPoints_randomLinearFixedQuadraticSlopes,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    se = FALSE,
    linewidth = 0.5) +
  geom_line(
    mapping = aes(
      x = age,
      y = fantasyPoints_randomLinearFixedQuadraticSlopes,
      color = positionFactor
    ),
    data = pointsPerSeason_positionAge_newData %>% 
      mutate(
        age = round(age, 2),
        fantasyPoints_randomLinearFixedQuadraticSlopes = round(fantasyPoints_randomLinearFixedQuadraticSlopes, 2)
        ),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position:\nModel With Random Intercepts, Random Linear Slopes, and\nFixed Quadratic Slopes",
    color = "Position"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsRandomLinearFixedQuadracticSlopes,
  tooltip = c("age","fantasyPoints_randomLinearFixedQuadraticSlopes","text","label")
)
```

###### Random Intercepts, Random Linear Slopes, Random Quadratic Slopes {#sec-fantasyPointsByAgeExperienceModelsAgeQuadraticRandom}

```{r}
pointsPerSeason_positionAgeRandomLinearRandomQuadraticSlopes <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + ageCentered20 + ageCentered20Quadratic + (1 + ageCentered20 + ageCentered20Quadratic | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)
```

###### Random Intercepts, Random Linear Slopes, Fixed Quadratic Slopes in Interaction With Position {#sec-fantasyPointsByAgeExperienceModelsAgeQuadraticFixedInteraction}

```{r}
pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + ageCentered20 + ageCentered20Quadratic + positionFactor:ageCentered20 + positionFactor:ageCentered20Quadratic + (1 + ageCentered20 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
print(effectsize::standardize_parameters(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction, method = "refit"), digits = 2)

performance::r2(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction, "positionFactor")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction, "ageCentered20")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes, "ageCentered20Quadratic")
performance::icc(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
AIC(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
```

A plot of the model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts, random linear slopes, and fixed quadratic slopes in interaction with position is in @fig-fantasyPointsRandomLinearFixedQuadraticSlopesInteraction.

```{r}
#| label: fig-fantasyPointsRandomLinearFixedQuadraticSlopesInteraction
#| fig-cap: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes in Interaction With Position."
#| fig-alt: "Plot of Model-Implied Trajectories of Fantasy Points by Age and Position in Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes in Interaction With Position."
#| code-fold: true

pointsPerSeason_positionAge_newData$fantasyPoints_randomLinearFixedQuadraticSlopesInteraction <- predict(
  object = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_randomLinearFixedQuadraticSlopesInteraction,
    color = positionFactor
  )
) + 
  geom_line(linewidth = 2) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Mixed Model with Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes\nin Interaction With Position",
    color = "Position"
  ) +
  theme_classic()
```

A plot of individuals' model-implied trajectories of fantasy points by age and position from the mixed model with random intercepts, random linear slopes, and fixed quadratic slopes in interaction with position is in @fig-individualFantasyPointsRandomLinearFixedQuadraticSlopesInteraction.

```{r}
#| label: fig-individualFantasyPointsRandomLinearFixedQuadraticSlopesInteraction
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes in Interaction With Position. Overlaid with the Model-Implied Trajectory by Position."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age and Position, from a Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes in Interaction With Position. Overlaid with the Model-Implied Trajectory by Position."
#| code-fold: true

player_stats_seasonal_offense_subsetCC$fantasyPoints_randomLinearFixedQuadraticSlopesInteraction <- predict(
  object = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction,
  newdata = player_stats_seasonal_offense_subsetCC
)

plot_individualFantasyPointsRandomLinearFixedQuadraticSlopesInteraction <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% 
    mutate(
      age = round(age, 2),
      fantasyPoints_randomLinearFixedQuadraticSlopesInteraction = round(fantasyPoints_randomLinearFixedQuadraticSlopesInteraction, 2)
    ),
  mapping = aes(
    x = age,
    y = fantasyPoints_randomLinearFixedQuadraticSlopesInteraction,
    color = positionFactor,
    group = player_id)) +
  geom_line(
    aes(
      x = age,
      y = fantasyPoints_randomLinearFixedQuadraticSlopesInteraction,
      text = player_display_name, # add player name for mouse over tooltip
      label = season # add season for mouse over tooltip
    ),
    linewidth = 0.5) +
  geom_line(
    mapping = aes(
      x = age,
      y = fantasyPoints_randomLinearFixedQuadraticSlopesInteraction,
      color = positionFactor
    ),
    data = pointsPerSeason_positionAge_newData %>% 
      mutate(
        age = round(age, 2),
        fantasyPoints_randomLinearFixedQuadraticSlopesInteraction = round(fantasyPoints_randomLinearFixedQuadraticSlopesInteraction, 2)
        ),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position:\nModel With Random Intercepts, Random Linear Slopes, and\nFixed Quadratic Slopes in Interaction With Position",
    color = "Position"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsRandomLinearFixedQuadraticSlopesInteraction,
  tooltip = c("age","fantasyPoints_randomLinearFixedQuadraticSlopesInteraction","text","label")
)
```

###### Adding Fixed-Effect Predictor of Experience {#sec-fantasyPointsByAgeExperienceModelsAgeQuadraticFixedInteractionExperience}

```{r}
pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + ageCentered20 + ageCentered20Quadratic + positionFactor:ageCentered20 + positionFactor:ageCentered20Quadratic + years_of_experience + (1 + ageCentered20 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  REML = FALSE,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
print(effectsize::standardize_parameters(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, method = "refit"), digits = 2)

performance::r2(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "positionFactor")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "ageCentered20")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "ageCentered20Quadratic")
emmeans::emmeans(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "years_of_experience")
performance::icc(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
AIC(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
```

##### Compare Models {#sec-fantasyPointsByAgeExperienceModelsAgeCompare}

After fitting several models, we now must compare their fit to determine which model fits "best" while also considering parsimony.
Parsimonious models are more likely to be true and more likely to generalize to other samples, because more complex models are more likely to overfit the data.
Thus, more complex models will almost always fit better than simpler models.
Thus, we are not just interested in whether a more complex model fits better than the simpler model; we also care about whether the more complex model fits *significantly* better than the simpler model given its additional complexity.
For evaluating and comparing models, we examine the [likelihood ratio test](#sec-fantasyPointsByAgeExperienceModelsAgeCompareLRT), the [Akaike Information Criterion](#sec-fantasyPointsByAgeExperienceModelsAgeCompareAIC) (AIC), the [corrected AIC](#sec-fantasyPointsByAgeExperienceModelsAgeCompareAICc) (AICc), the [Bayesian Information Criterion](#sec-fantasyPointsByAgeExperienceModelsAgeCompareBIC) (BIC), [$R^2$](#sec-fantasyPointsByAgeExperienceModelsAgeCompareRsquared), [deviance](#sec-fantasyPointsByAgeExperienceModelsAgeCompareDeviance), and [log likelihood](#sec-fantasyPointsByAgeExperienceModelsAgeCompareLogLik).

The BIC penalizes model complexity more than the AIC does.
The BIC is preferable when there is a "true" model, and one intends to identify the true model.
The AIC is preferable when we are concerned more about predictive accuracy and when overfitting is less of a concern.
Because we are more concerned about predictive accuracy and we do not believe one of these models is the "true" model per se of age-related changes in fantasy performance, we will give more weight to AIC than BIC.

Below, we specify various groups of models for the model fit comparisons:

```{r}
lmVsMixedModel <- list(
  "nullModel" = pointsPerSeason_nullModel,
  "randomIntercepts" = pointsPerSeason_randomIntercepts
)

lmAndMixedModels <- list(
  "nullModel" = pointsPerSeason_nullModel,
  "linearRegression" = pointsPerSeason_linearRegression,
  "quadraticRegression" = pointsPerSeason_quadraticRegression,
  "randomIntercepts" = pointsPerSeason_randomIntercepts,
  "position" = pointsPerSeason_position,
  "fixedLinear" = pointsPerSeason_positionAgeFixedLinearSlopes,
  "randomLinear" = pointsPerSeason_positionAgeRandomLinearSlopes,
  "randomLinearFixedQuadratic" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  "randomLinearFixedQuadraticInteraction" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction
)

mixedModels <- list(
  "randomIntercepts" = pointsPerSeason_randomIntercepts,
  "position" = pointsPerSeason_position,
  "fixedLinear" = pointsPerSeason_positionAgeFixedLinearSlopes,
  "randomLinear" = pointsPerSeason_positionAgeRandomLinearSlopes,
  "randomLinearFixedQuadratic" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  "randomLinearFixedQuadraticInteraction" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction
)

mixedModels1 <- list(
  "randomIntercepts" = pointsPerSeason_randomIntercepts,
  "position" = pointsPerSeason_position
)

mixedModels2 <- list(
  "fixedLinear" = pointsPerSeason_positionAgeFixedLinearSlopes,
  "randomLinear" = pointsPerSeason_positionAgeRandomLinearSlopes,
  "randomLinearFixedQuadratic" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  "randomLinearFixedQuadraticInteraction" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction
)
```

###### Likelihood Ratio Test {#sec-fantasyPointsByAgeExperienceModelsAgeCompareLRT}

```{r}
anova(
  pointsPerSeason_randomIntercepts,
  pointsPerSeason_position,
  pointsPerSeason_positionAgeFixedLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction
)
```

###### Akaike Information Criterion (AIC) {#sec-fantasyPointsByAgeExperienceModelsAgeCompareAIC}

```{r}
AIC(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression,
  pointsPerSeason_randomIntercepts,
  pointsPerSeason_positionAgeFixedLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction
  )

bbmle::AICtab(lmAndMixedModels)
```

###### Corrected Akaike Information Criterion (AICc) {#sec-fantasyPointsByAgeExperienceModelsAgeCompareAICc}

We compute model fit using the using the `AICcmodavg` [@R-AICcmodavg], `bbmle` [@R-bbmle], and `MuMIn` [@R-MuMIn] packages

```{r}
#AICcmodavg::aictab(lmVsMixedModel) # throws error (can't mix lm with lmer)
bbmle::AICctab(lmVsMixedModel)

AICcmodavg::aictab(mixedModels) # throws error (can't mix lm with lmer)
#bbmle::AICctab(mixedModels) # throws error (different numbers of observations)

AICcmodavg::aictab(mixedModels1)
bbmle::AICctab(mixedModels1)

AICcmodavg::aictab(mixedModels2)
bbmle::AICctab(mixedModels2)

MuMIn::AICc(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression,
  pointsPerSeason_randomIntercepts,
  pointsPerSeason_positionAgeFixedLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction
)
```

###### Bayesian Information Criterion (BIC) {#sec-fantasyPointsByAgeExperienceModelsAgeCompareBIC}

```{r}
BIC(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression,
  pointsPerSeason_randomIntercepts,
  pointsPerSeason_positionAgeFixedLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction
  )

#AICcmodavg::bictab(lmAndMixedModels) # throws error (can't mix lm with lmer)
bbmle::BICtab(lmVsMixedModel)

AICcmodavg::bictab(mixedModels)
#bbmle::AICctab(mixedModels) # throws error (different numbers of observations)

AICcmodavg::bictab(mixedModels1)
bbmle::BICtab(mixedModels1)

AICcmodavg::bictab(mixedModels2)
bbmle::BICtab(mixedModels2)
```

###### $R^2$ {#sec-fantasyPointsByAgeExperienceModelsAgeCompareRsquared}

```{r}
summary(pointsPerSeason_nullModel)$r.squared
summary(pointsPerSeason_linearRegression)$r.squared
summary(pointsPerSeason_quadraticRegression)$r.squared
performance::r2(pointsPerSeason_randomIntercepts)
performance::r2(pointsPerSeason_positionAgeFixedLinearSlopes)
performance::r2(pointsPerSeason_positionAgeRandomLinearSlopes)
performance::r2(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
performance::r2(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
```

###### Deviance {#sec-fantasyPointsByAgeExperienceModelsAgeCompareDeviance}

```{r}
deviance(pointsPerSeason_nullModel)
deviance(pointsPerSeason_linearRegression)
deviance(pointsPerSeason_quadraticRegression)
deviance(pointsPerSeason_randomIntercepts)
deviance(pointsPerSeason_positionAgeFixedLinearSlopes)
deviance(pointsPerSeason_positionAgeRandomLinearSlopes)
deviance(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
deviance(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
```

###### Log Likelihood {#sec-fantasyPointsByAgeExperienceModelsAgeCompareLogLik}

```{r}
logLik(pointsPerSeason_nullModel)
logLik(pointsPerSeason_linearRegression)
logLik(pointsPerSeason_quadraticRegression)
logLik(pointsPerSeason_randomIntercepts)
logLik(pointsPerSeason_positionAgeFixedLinearSlopes)
logLik(pointsPerSeason_positionAgeRandomLinearSlopes)
logLik(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
logLik(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
```

#### Multicollinearity {#sec-mixedModelsMulticollinearity}

As with [multiple regression](#sec-multipleRegression), [multicollinearity](#sec-multipleRegressionMulticollinearity) is also a potential concern for mixed models.
Here are the variance inflation factors (VIFs) from an earlier model (VIF is described in @sec-multipleRegressionMulticollinearity):

```{r}
car::vif(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
```

#### Diagnostic Plots {#sec-fantasyPointsByAgeExperienceModelsDiagnosticPlots}

##### Regression Coefficients {#sec-fantasyPointsByAgeExperienceModelsRegressionCoefficients}

To visualize the regression coefficients, we can use the `broom.mixed::tidy()` function of the `broom.mixed` package [@R-broom.mixed], as in @fig-individualFantasyPointsRandomLinearFixedQuadraticSlopesInteractionCoefficients.

```{r}
#| label: fig-individualFantasyPointsRandomLinearFixedQuadraticSlopesInteractionCoefficients
#| fig-cap: "Regression Coefficients with 95% Confidence Interval from a Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes in Interaction With Position. Overlaid with the Model-Implied Trajectory by Position."
#| fig-alt: "Regression Coefficients with 95% Confidence Interval from a Mixed Model With Random Intercepts, Random Linear Slopes, and Fixed Quadratic Slopes in Interaction With Position. Overlaid with the Model-Implied Trajectory by Position."
#| code-fold: true

pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience_tidy <- broom.mixed::tidy(
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience,
  effects = "fixed",
  conf.int = TRUE)

ggplot(
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience_tidy,
  aes(
    x = term,
    y = estimate)) +
  geom_point() +
  geom_errorbar(
    aes(
      ymin = conf.low,
      ymax = conf.high),
    width = 0.2) +
  coord_flip() + # flip axes for readability
  labs(
    title = "Regression Coefficients with 95% CI",
    y = "Coefficient Estimate",
    x = "") +
  theme_minimal()
```

##### Plot of Residuals {#sec-fantasyPointsByAgeExperienceModelsResidualPlot}

```{r}
#| label: fig-residualPlotMixedModel
#| fig-cap: "Plot of Residuals."
#| fig-alt: "Plot of Residuals."

plot(
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience,
  type = c("p","smooth"),
  col.line = 1)
```

##### Scale-Location Plot {#sec-fantasyPointsByAgeExperienceModelsScaleLocationPlot}

```{r}
#| label: fig-scaleLocationPlotMixedModel
#| fig-cap: "Scale-Location Plot."
#| fig-alt: "Scale-Location Plot."

plot(
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience,
  sqrt(abs(resid(.))) ~ fitted(.),
  type = c("p","smooth"),
  col.line = 1)
```

##### Quantile–Quantile (QQ) Plot {#sec-fantasyPointsByAgeExperienceModelsQQPlot}

```{r}
#| label: fig-qqPlotMixedModel
#| fig-cap: "Quantile–Quantile (QQ) Plot."
#| fig-alt: "Quantile–Quantile (QQ) Plot."

qqnorm(resid(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience))
qqline(resid(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience))
```

##### Probability–Probability (PP) Plot {#sec-fantasyPointsByAgeExperienceModelsDiagnosticPPPlots}

```{r}
#| label: fig-ppPlotMixedModel
#| fig-cap: "Probability–Probability (PP) Plot."
#| fig-alt: "Probability–Probability (PP) Plot."

petersenlab::ppPlot(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
```

##### Residuals Versus Leverage Plot {#sec-fantasyPointsByAgeExperienceModelsResidualsVsLeveragePlot}

```{r}
#| label: fig-residualsVsLeveragePlotMixedModel
#| fig-cap: "Residuals Versus Leverage Plot."
#| fig-alt: "Residuals Versus Leverage Plot."

plot(
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience,
  rstudent(.) ~ hatvalues(.))
```

##### Other Diagnostic Plots {#sec-fantasyPointsByAgeExperienceModelsOtherDiagnosticPlots}

```{r}
#| label: fig-otherDiagnosticPlotsMixedModel
#| fig-cap: "Other Diagnostic Plots."
#| fig-alt: "Other Diagnostic Plots."
#| fig-height: 10

performance::check_model(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
```

#### Generalized Additive Model {#sec-fantasyPointsByAgeExperienceModelsAgeGAM}

To fit a generalized additive model, we first determine how many cores are available for multicore processing.

```{r}
#| label: parallel-cores-mixed-models
#| include: false

num_cores <- parallelly::availableCores()
num_true_cores <- parallelly::availableCores(logical = FALSE)
```

```{r}
#| eval: false

num_cores <- parallelly::availableCores() - 1
num_true_cores <- parallelly::availableCores(logical = FALSE) - 1
```

We fit a generalized additive model using the `mgcv::bam()` function of the `mgcv` package [@R-mgcv; @Wood2017].

```{r}
#| label: gam-mixed-models

pointsPerSeason_gam <- mgcv::bam( # using bam() instead of gam() for faster estimation due to large size of data
  fantasyPoints ~ positionFactor + s(ageCentered20, by = positionFactor) + years_of_experience + s(player_idFactor, ageCentered20, bs = "re"),
  data = player_stats_seasonal_offense_subset,
  nthreads = num_cores
)
```

```{r}
#| label: gam-mixed-models-summary

pointsPerSeason_gamSummary <- summary(pointsPerSeason_gam)

pointsPerSeason_gamSummary
pointsPerSeason_gamSummary$r.sq
performance::r2(pointsPerSeason_gam)
AIC(pointsPerSeason_gam)
```

##### Compare Models {#sec-fantasyPointsByAgeExperienceModelsAgeGAMcompare}

```{r}
linearMixedModelsVsGAM <- list(
  "fixedLinear" = pointsPerSeason_positionAgeFixedLinearSlopes,
  "randomLinear" = pointsPerSeason_positionAgeRandomLinearSlopes,
  "randomLinearFixedQuadratic" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  "randomLinearFixedQuadraticInteraction" = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction,
  "gam" = pointsPerSeason_gam
)

AIC(
  #pointsPerSeason_nullModel,
  #pointsPerSeason_linearRegression,
  #pointsPerSeason_quadraticRegression,
  #pointsPerSeason_randomIntercepts,
  pointsPerSeason_positionAgeFixedLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction,
  pointsPerSeason_gam
)

#AICcmodavg::aictab(linearMixedModelsVsGAM) # throws error (can't mix bam with lmer)
#bbmle::AICctab(linearMixedModelsVsGAM) # different numbers of observations

MuMIn::AICc(
  pointsPerSeason_nullModel,
  pointsPerSeason_linearRegression,
  pointsPerSeason_quadraticRegression,
  pointsPerSeason_randomIntercepts,
  pointsPerSeason_positionAgeFixedLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction,
  pointsPerSeason_gam
)

BIC(
  #pointsPerSeason_nullModel,
  #pointsPerSeason_linearRegression,
  #pointsPerSeason_quadraticRegression,
  #pointsPerSeason_randomIntercepts,
  pointsPerSeason_positionAgeFixedLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes,
  pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction,
  pointsPerSeason_gam
)

#AICcmodavg::bictab(linearMixedModelsVsGAM) # throws error (can't mix bam with lmer)
#bbmle::AICctab(linearMixedModelsVsGAM) # different numbers of observations

summary(pointsPerSeason_nullModel)$r.squared
summary(pointsPerSeason_linearRegression)$r.squared
summary(pointsPerSeason_quadraticRegression)$r.squared
performance::r2(pointsPerSeason_randomIntercepts)
performance::r2(pointsPerSeason_positionAgeFixedLinearSlopes)
performance::r2(pointsPerSeason_positionAgeRandomLinearSlopes)
performance::r2(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
performance::r2(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
performance::r2(pointsPerSeason_gam)

deviance(pointsPerSeason_nullModel)
deviance(pointsPerSeason_linearRegression)
deviance(pointsPerSeason_quadraticRegression)
deviance(pointsPerSeason_randomIntercepts)
deviance(pointsPerSeason_positionAgeFixedLinearSlopes)
deviance(pointsPerSeason_positionAgeRandomLinearSlopes)
deviance(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
deviance(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
deviance(pointsPerSeason_gam)

logLik(pointsPerSeason_nullModel)
logLik(pointsPerSeason_linearRegression)
logLik(pointsPerSeason_quadraticRegression)
logLik(pointsPerSeason_randomIntercepts)
logLik(pointsPerSeason_positionAgeFixedLinearSlopes)
logLik(pointsPerSeason_positionAgeRandomLinearSlopes)
logLik(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopes)
logLik(pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteraction)
logLik(pointsPerSeason_gam)
```

##### Players Who Were (at Least Once) at the Top of the End-of-Season Depth Chart {#sec-fantasyPointsByAgeExperienceModelsAgeDepth}

```{r}
#| label: players-top-depth-chart-mixed-models

pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience <- lmerTest::lmer(
  fantasyPoints ~ positionFactor + ageCentered20 + ageCentered20Quadratic + positionFactor:ageCentered20 + positionFactor:ageCentered20Quadratic + years_of_experience + (1 + ageCentered20 | player_idFactor),
  data = player_stats_seasonal_offense_subset,
  control = lmerControl(optimizer = "bobyqa")
)

summary(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
petersenlab::lm.beta.lmer(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)

performance::r2(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
emmeans::emmeans(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "positionFactor")
emmeans::emmeans(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "ageCentered20")
emmeans::emmeans(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "ageCentered20Quadratic")
emmeans::emmeans(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience, "years_of_experience")
performance::icc(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
AIC(pointsPerSeasonDepth_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience)
```

```{r}
#| label: gam-players-top-depth-chart-mixed-models

pointsPerSeasonDepth_gam <- bam( # using bam() instead of gam() for faster estimation due to large size of data
  fantasyPoints ~ positionFactor + s(ageCentered20, by = positionFactor) + years_of_experience + s(player_idFactor, ageCentered20, bs = "re"),
  data = player_stats_seasonal_offense_subsetDepth,
  nthreads = num_cores
)
```

```{r}
#| label: gam-players-top-depth-chart-mixed-models-summary

pointsPerSeasonDepth_gamSummary <- summary(pointsPerSeasonDepth_gam)

pointsPerSeasonDepth_gamSummary
pointsPerSeasonDepth_gamSummary$r.sq
performance::r2(pointsPerSeasonDepth_gam)
AIC(pointsPerSeasonDepth_gam)
```

### Bayesian Mixed Models {#sec-fantasyPointsByAgeExperienceModelsBayesian}

#### Determine Response Distribution {#sec-fantasyPointsByAgeExperienceModelsBayesianDistribution}

```{r}
player_stats_seasonal_offense_subset$fantasyPoints_posOnlyNoZeros <- player_stats_seasonal_offense_subset$fantasyPoints_posOnly <- player_stats_seasonal_offense_subset$fantasyPoints
player_stats_seasonal_offense_subset$fantasyPoints_posOnly[player_stats_seasonal_offense_subset$fantasyPoints < 0] <- 0
player_stats_seasonal_offense_subset$fantasyPoints_posOnlyNoZeros[player_stats_seasonal_offense_subset$fantasyPoints <= 0] <- 0.01

fantasyPointsVector <- player_stats_seasonal_offense_subset$fantasyPoints %>% 
    na.omit() %>% 
    as.vector()

fantasyPointsVector_posOnlyNoZeros <- fantasyPointsVector_posOnly <- fantasyPointsVector
fantasyPointsVector_posOnly[fantasyPointsVector < 0] <- 0
fantasyPointsVector_posOnlyNoZeros[fantasyPointsVector <= 0] <- 0.01
```

```{r}
#| label: fig-histogramFantasyPoints
#| fig-cap: "Histogram of Fantasy Points with Overlaid Density Plot and Rug Plot."
#| fig-alt: "Histogram of Fantasy Points with Overlaid Density Plot and Rug Plot."

ggplot2::ggplot(
  data = player_stats_seasonal_offense_subset,
  mapping = aes(
    x = fantasyPoints)
) +
  geom_histogram(
    aes(y = after_stat(density)),
    color = "#000000",
    fill = "#0099F8"
  ) +
  geom_density(
    color = "#000000",
    fill = "#F85700",
    alpha = 0.6 # add transparency
  ) +
  geom_rug() +
  theme_classic()
```

We can use the `fitdistrplus` package [@R-fitdistrplus; @DelignetteMuller2015_packages]:

```{r}
#| label: fig-cullenFreyGraphy
#| fig-cap: "Cullen and Frey Graph."
#| fig-alt: "Cullen and Frey Graph."

fitdistrplus::descdist(
  fantasyPointsVector,
  boot = 1000)
```

```{r}
# all values
fit.norm <- fitdistrplus::fitdist(fantasyPointsVector, "norm")

# positive-only
fit.exp <- fitdist(fantasyPointsVector_posOnly, "exp")

# positive and no zeros
fit.gamma <- fitdistrplus::fitdist(fantasyPointsVector_posOnlyNoZeros, "gamma")
fit.lognormal <- fitdistrplus::fitdist(fantasyPointsVector_posOnlyNoZeros, "lnorm")
fit.weibull <- fitdistrplus::fitdist(fantasyPointsVector_posOnlyNoZeros, "weibull")

# Model fit
AIC(fit.norm)
AIC(fit.exp)
AIC(fit.gamma) # fits best
AIC(fit.lognormal)
AIC(fit.weibull)
```

```{r}
#| label: fig-normalDistribution
#| fig-cap: "Fit of normal distribution to fantasy points."
#| fig-alt: "Fit of normal distribution to fantasy points."

plot(fit.norm)
```

```{r}
#| label: fig-logNormalDistribution
#| fig-cap: "Fit of log normal distribution to fantasy points."
#| fig-alt: "Fit of log normal distribution to fantasy points."

plot(fit.lognormal)
```

```{r}
#| label: fig-gammaDistribution
#| fig-cap: "Fit of gamma distribution to fantasy points."
#| fig-alt: "Fit of gamma distribution to fantasy points."

plot(fit.gamma)
```

```{r}
#| label: fig-exponentialDistribution
#| fig-cap: "Fit of exponential distribution to fantasy points."
#| fig-alt: "Fit of exponential distribution to fantasy points."

plot(fit.exp)
```

```{r}
#| label: fig-weibullDistribution
#| fig-cap: "Fit of Weibull distribution to fantasy points."
#| fig-alt: "Fit of Weibull distribution to fantasy points."

plot(fit.weibull)
```

```{r}
model_normal <- glm(
  formula = fantasyPoints ~ 1,
  family = gaussian(),
  data = player_stats_seasonal_offense_subset,
  maxit = 100000)

model_gamma <- glm(
  formula = fantasyPoints_posOnlyNoZeros ~ 1,
  family = Gamma(),
  data = player_stats_seasonal_offense_subset,
  maxit = 100000)

sum(resid(model_normal)^2)
sum(resid(model_gamma)^2)

AIC(model_normal)
AIC(model_gamma)
```

Because the response distribution (i.e., fantasy points) is positively skewed, we will use a gamma response distribution.

#### Specify Model Formula {#sec-fantasyPointsByAgeExperienceModelsBayesianFormula}

Information about smooth terms in the `mgcv` package [@R-mgcv; @Wood2017] is provided at the following link: <https://stat.ethz.ch/R-manual/R-devel/library/mgcv/html/smooth.terms.html>.
We estimate the model using the `brms` package [@R-brms; @Buerkner2017_packages; @Buerkner2018_packages], which provides an interface to the `rstan` package [@R-rstan].

Specify model formula:

```{r}
#| eval: false

bayesianMixedModelFormula <- brms::bf(
  fantasyPoints_posOnly ~ positionFactor + s(ageCentered20, by = positionFactor) + years_of_experience + s(ageCentered20, player_idFactor, bs = "re") + (1 | player_idFactor)
)
```

#### Run Model {#sec-fantasyPointsByAgeExperienceModelsBayesianRunModel}

Now, we can run the model.

::: {#nte-bayesianMixedModel .callout-note title="Bayesian mixed model"}
Note: the following code that runs the model takes a while.
If you just want to save time and load the model object instead of running the model, you can load the model object (which has already been fit) using this code:

```{r}
#| eval: false

load(url("https://osf.io/download/q6rjf/"))
```
:::

We use the `cmdstanr` backend [@R-cmdstanr] and threading using the `parellely` package [@R-parallelly] for parallel (faster) processing:

```{r}
#| eval: false

bayesianMixedModelFit <- brms::brm(
  formula = bayesianMixedModelFormula,
  data = player_stats_seasonal_offense_subset,
  family = hurdle_gamma(),
  cores = 4,
  save_pars = save_pars(latent = FALSE, all = FALSE),
  threads = threading(parallelly::availableCores()),
  backend = "cmdstanr",
  seed = 52242,
  silent = 0
)
```

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

save(
  bayesianMixedModelFit,
  file = "./data/bayesianMixedModelFit.RData")
```

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

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

```{r}
#| include: false

petersenlab::robust_load("https://osf.io/download/q6rjf/") #load(url("https://osf.io/download/q6rjf/"))
```

#### Model Summary {#sec-fantasyPointsByAgeExperienceModelsBayesianModelSummary}

```{r}
summary(bayesianMixedModelFit)
brms::prior_summary(bayesianMixedModelFit)
```

::: {#nte-bayesianMixedModelR2 .callout-note title="Bayesian mixed model R-squared"}
Note: the following code that estimates the $R^2$ for the Bayesian mixed model takes a while to run.
:::

```{r}
#| eval: false

performance::r2(bayesianMixedModelFit)
```

#### Trace Plots {#sec-fantasyPointsByAgeExperienceModelsBayesianTrace}

```{r}
#| label: fig-bayesianTracePlots
#| fig-cap: "Trace Plots from Bayesian Mixed Model."
#| fig-alt: "Trace Plots from Bayesian Mixed Model."

plot(bayesianMixedModelFit, ask = FALSE)
```

#### Posterior Predictive Check {#sec-fantasyPointsByAgeExperienceModelsBayesianPPcheck}

```{r}
#| label: fig-bayesianPosteriorPredictiveCheck
#| fig-cap: "Posterior Predictive Check from Bayesian Mixed Model."
#| fig-alt: "Posterior Predictive Check from Bayesian Mixed Model."

pp_check(bayesianMixedModelFit) + 
  ggplot2::xlim(0, 600)
```

#### Plot of Regression Coefficients {#sec-fantasyPointsByAgeExperienceModelsBayesianRegressionCoefficients}

To visualize the regression coefficients, we can use the `broom.mixed::tidy()` function of the `broom.mixed` package [@R-broom.mixed], as in @fig-bayesianRegressionCoefficients.

```{r}
#| label: fig-bayesianRegressionCoefficients
#| fig-cap: "Regression Coefficients with 95% Credible Interval from Bayesian Mixed Model."
#| fig-alt: "Regression Coefficients with 95% Credible Interval from Bayesian Mixed Model."

bayesianMixedModelFit_tidy <- broom.mixed::tidy(
  bayesianMixedModelFit,
  effects = "fixed",
  conf.int = TRUE)

ggplot(
  bayesianMixedModelFit_tidy,
  aes(
    x = term,
    y = estimate)) +
  geom_point() +
  geom_errorbar(
    aes(
      ymin = conf.low,
      ymax = conf.high),
    width = 0.2) +
  coord_flip() + # flip axes for readability
  labs(
    title = "Posterior Distributions and 95% Credible Intervals",
    x = "Predictor",
    y = "Posterior Estimate") +
  theme_minimal()
```

To visualize the distribution of regression coefficients, we can use the `tidybayes` package [@R-tidybayes], as in @fig-bayesianRegressionCoefficientsDistributions.

```{r}
#| label: fig-bayesianRegressionCoefficientsDistributions
#| fig-cap: "Distributions of Regression Coefficients with 66% (Thick Line) and 95% (Thin Line) Credible Intervals from Bayesian Mixed Model."
#| fig-alt: "Distributions of Regression Coefficients with 66% (Thick Line) and 95% (Thin Line) Credible Intervals from Bayesian Mixed Model."

bayesianMixedModelFit %>%
  tidybayes::tidy_draws() %>%
  select(starts_with("b_")) %>%
  pivot_longer(
    everything(),
    names_to = "term",
    values_to = "estimate") %>%
  mutate(term = sub("^b_", "", term)) %>%
  ggplot(
    aes(
      x = estimate,
      y = term)) +
  stat_halfeye() + # by default, .width shows the 66% (thick line) and 95% (thin line) credible intervals
  labs(
    title = "Posterior Distributions of Regression Coefficients\nwith 66% and 95% Credible Intervals",
    x = "Predictor",
    y = "Posterior Estimate") +
  theme_minimal()
```

### Plots of Model-Implied Fantasy Points by Position and Age {#sec-fantasyPointsByAgePlots}

```{r}
#| label: generate-model-implied-predictions-fantasy-points-by-position-age

# From Quadratic Model: All Players
pointsPerSeason_positionAge_newData$fantasyPoints_quadratic <- predict(
  object = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

# From Quadratic Model: Players at Top of End-of-Season Depth Chart
pointsPerSeason_positionAge_newData$fantasyPoints_depthQuadratic <- predict(
  object = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience,
  newdata = pointsPerSeason_positionAge_newData,
  re.form = NA
)

# From GAM Model: All Players
pointsPerSeason_positionAge_newData$fantasyPoints_gam <- predict(
  object = pointsPerSeason_gam,
  newdata = pointsPerSeason_positionAge_newData,
  newdata.guaranteed = TRUE,
  exclude = "s(player_idFactor,ageCentered20)"
)

# From GAM Model: Players at Top of End-of-Season Depth Chart
pointsPerSeason_positionAge_newData$fantasyPoints_depthGAM <- predict(
  object = pointsPerSeasonDepth_gam,
  newdata = pointsPerSeason_positionAge_newData,
  newdata.guaranteed = TRUE,
  exclude = "s(player_idFactor,ageCentered20)"
)
```

```{r}
#| label: generate-model-implied-predictions-fantasy-points-by-position-age-bayesian1

# From Bayesian Model: All Players
model_levels <- levels(bayesianMixedModelFit$data$player_idFactor)

pointsPerSeason_positionAge_newData$player_idFactor <- factor(
  rep(model_levels[1], nrow(pointsPerSeason_positionAge_newData)), 
  levels = model_levels
)
```

```{r}
#| label: generate-model-implied-predictions-fantasy-points-by-position-age-bayesian2
#| eval: false

modelPredictions_bayesian <- predict(
  object = bayesianMixedModelFit,
  newdata = pointsPerSeason_positionAge_newData,
  re_formula = NA,
  cores = num_cores
)
```

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

save(
  modelPredictions_bayesian,
  file = "./data/modelPredictions_bayesian.RData")
```

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

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

```{r}
#| include: false

load(file = "./data/modelPredictions_bayesian.RData")
```

```{r}
#| label: generate-model-implied-predictions-fantasy-points-by-position-age-bayesian3

pointsPerSeason_positionAge_newData$fantasyPoints_bayesian <- modelPredictions_bayesian[,"Estimate"]
```

Plots of model-implied fantasy points by position and age are in Figures [-@fig-fantasyPointsByAgeQuadratic]–[-@fig-fantasyPointsByAgeBayesian].

#### Quadratic Model {#sec-fantasyPointsByAgePlotsQuadratic}

```{r}
#| label: fig-fantasyPointsByAgeQuadratic
#| fig-cap: "Plot of Model-Implied Quadratic Trajectories of Fantasy Points by Age."
#| fig-alt: "Plot of Model-Implied Quadratic Trajectories of Fantasy Points by Age."
#| code-fold: true

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_quadratic,
    color = positionFactor
  )
) + 
  geom_smooth() +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Quadratic Model with All Players",
    color = "Position"
  ) +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(fill = NA))) # transparent legend background
```

#### Quadratic Model: Top of Depth Chart {#sec-fantasyPointsByAgePlotsQuadraticDepth}

```{r}
#| label: fig-fantasyPointsByAgeQuadraticDepth
#| fig-cap: "Plot of Model-Implied Quadratic Trajectories of Fantasy Points by Age For Players Who Were Once at the Top of the End-of-Season Depth Chart."
#| fig-alt: "Plot of Model-Implied Quadratic Trajectories of Fantasy Points by Age For Players Who Were Once at the Top of the End-of-Season Depth Chart."
#| code-fold: true

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_depthQuadratic,
    color = positionFactor
  )
) + 
  geom_smooth() +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Quadratic Model with Players Who Were Once at Top of End-of-Season Depth Chart",
    color = "Position"
  ) +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(fill = NA))) # transparent legend background
```

#### Generalized Additive Model {#sec-fantasyPointsByAgePlotsGAM}

```{r}
#| label: fig-fantasyPointsByAgeGAM
#| fig-cap: "Plot of Implied Trajectories of Fantasy Points by Age from a Generalized Additive Model."
#| fig-alt: "Plot of Implied Trajectories of Fantasy Points by Age from a Generalized Additive Model."
#| code-fold: true

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_gam,
    color = positionFactor
  )
) + 
  geom_smooth() +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Generalized Additive Model with All Players",
    color = "Position"
  ) +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(fill = NA))) # transparent legend background
```

#### Generalized Additive Model: Top of Depth Chart {#sec-fantasyPointsByAgePlotsGAMdepth}

```{r}
#| label: fig-fantasyPointsByAgeGAMdepth
#| fig-cap: "Plot of Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, For Players Who Were Once at the Top of the End-of-Season Depth Chart."
#| fig-alt: "Plot of Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, For Players Who Were Once at the Top of the End-of-Season Depth Chart."
#| code-fold: true

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_depthGAM,
    color = positionFactor
  )
) + 
  geom_smooth() +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Generalized Additive Model with Players Who Were Once at Top of End-of-Season Depth Chart",
    color = "Position"
  ) +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(fill = NA))) # transparent legend background
```

#### Bayesian Mixed Model {#sec-fantasyPointsByAgePlotsBayesian}

```{r}
#| label: fig-fantasyPointsByAgeBayesian
#| fig-cap: "Plot of Implied Trajectories of Fantasy Points by Age from a Bayesian Mixed Model."
#| fig-alt: "Plot of Implied Trajectories of Fantasy Points by Age from a Bayesian Mixed Model."
#| code-fold: true

ggplot2::ggplot(
  data = pointsPerSeason_positionAge_newData,
  mapping = aes(
    x = age,
    y = fantasyPoints_bayesian,
    color = positionFactor
  )
) + 
  geom_smooth() +
  labs(
    x = "Player Age (years)",
    y = "Fantasy Points (Season)",
    title = "Fantasy Points (Season) by Player Age and Position",
    subtitle = "Bayesian Mixed Model with All Players",
    color = "Position"
  ) +
  theme_classic() +
  guides(color = guide_legend(override.aes = list(fill = NA))) # transparent legend background
```

### Plots of Individuals' Model-Implied Fantasy Points by Age {#sec-fantasyPointsIndividualByAgePlots}

We provide plots of individuals' model-implied fantasy points by age from the Bayesian mixed model in @sec-timeSeriesBayesianModelImpliedPredictions.

```{r}
player_stats_seasonal_offense_subsetCC$fantasyPoints_quadratic <- predict(
  object = pointsPerSeason_positionAgeRandomLinearFixedQuadraticSlopesInteractionExperience,
  newdata = player_stats_seasonal_offense_subsetCC
)

player_stats_seasonal_offense_subsetCC$fantasyPoints_gam <- predict(
  object = pointsPerSeason_gam,
  newdata = player_stats_seasonal_offense_subsetCC
)
```

```{r}
#| code-fold: true

zeroAge <- pointsPerSeason_positionAge_newData %>% 
  group_by(positionFactor) %>% 
  filter(fantasyPoints_gam < 0) %>% 
  slice(which.min(age))

peakAge <- pointsPerSeason_positionAge_newData %>% 
  group_by(positionFactor) %>% 
  slice(which.max(fantasyPoints_gam))

peakAge2 <- pointsPerSeason_positionAge_newData %>% 
  filter(age > 22) %>% 
  group_by(positionFactor) %>% 
  slice(which.max(fantasyPoints_gam))

qbPeakAge <- round(peakAge$age[which(peakAge$positionFactor == "QB")], 0)
fbPeakAge <- round(peakAge$age[which(peakAge$positionFactor == "FB")], 0)
rbPeakAge <- round(peakAge$age[which(peakAge$positionFactor == "RB")], 0)
wrPeakAge <- round(peakAge$age[which(peakAge$positionFactor == "WR")], 0)
wrPeakAge2 <- round(peakAge2$age[which(peakAge$positionFactor == "WR")], 0)
tePeakAge <- round(peakAge$age[which(peakAge$positionFactor == "TE")], 0)

qbZeroAge <- round(zeroAge$age[which(zeroAge$positionFactor == "QB")], 0)
fbZeroAge <- round(zeroAge$age[which(zeroAge$positionFactor == "FB")], 0)
rbZeroAge <- round(zeroAge$age[which(zeroAge$positionFactor == "RB")], 0)
wrZeroAge <- round(zeroAge$age[which(zeroAge$positionFactor == "WR")], 0)
teZeroAge <- round(zeroAge$age[which(zeroAge$positionFactor == "TE")], 0)
```

#### Quarterbacks {#sec-fantasyPointsIndividualByAgePlotsQB}

A plot of Quarterbacks' model-implied fantasy points by age is in @fig-individualFantasyPointsByAgeQB.
The model-implied peak of Quarterbacks' fantasy points is at age `{r} qbPeakAge`.
The model-predicted value of zero fantasy points for Quarterbacks is at age `{r} qbZeroAge`.

```{r}
#| label: fig-individualFantasyPointsByAgeQB
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, for Quarterbacks. Overlaid with the Model-Implied Trajectory for Quarterbacks in Blue."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, for Quarterbacks. Overlaid with the Model-Implied Trajectory for Quarterbacks in Blue."
#| code-fold: true

plot_individualFantasyPointsByAgeQB <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% filter(position == "QB"),
  mapping = aes(
    x = age,
    y = fantasyPoints_gam,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_smooth(
    mapping = aes(
      x = age,
      y = fantasyPoints_gam
    ),
    data = pointsPerSeason_positionAge_newData %>% filter(positionFactor == "QB"),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_gam","text","label")
)
```

#### Fullbacks {#sec-fantasyPointsIndividualByAgePlotsFB}

A plot of Fullbacks' model-implied fantasy points by age is in @fig-individualFantasyPointsByAgeFB.
The model-implied peak of Fullbacks' fantasy points is at age `{r} fbPeakAge`.
The model-predicted value of zero fantasy points for Fullbacks is at age `{r} fbZeroAge`.

```{r}
#| label: fig-individualFantasyPointsByAgeFB
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, for Fullbacks. Overlaid with the Model-Implied Trajectory for Fullbacks in Blue."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, for Fullbacks. Overlaid with the Model-Implied Trajectory for Fullbacks in Blue."
#| code-fold: true

plot_individualFantasyPointsByAgeFB <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% filter(position == "FB"),
  mapping = aes(
    x = age,
    y = fantasyPoints_gam,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_smooth(
    mapping = aes(
      x = age,
      y = fantasyPoints_gam
    ),
    data = pointsPerSeason_positionAge_newData %>% filter(positionFactor == "FB"),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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: Fullbacks"
  ) +
  theme_classic()

plotly::ggplotly(
  plot_individualFantasyPointsByAgeFB,
  tooltip = c("age","fantasyPoints_gam","text","label")
)
```

#### Running Backs {#sec-fantasyPointsIndividualByAgePlotsRB}

A plot of Running Backs' model-implied fantasy points by age is in @fig-individualFantasyPointsByAgeRB.
The model-implied peak of Running Backs' fantasy points is at age `{r} rbPeakAge`.
The model-predicted value of zero fantasy points for Running Backs is at age `{r} rbZeroAge`.

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

plot_individualFantasyPointsByAgeRB <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% filter(position == "RB"),
  mapping = aes(
    x = age,
    y = fantasyPoints_gam,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_smooth(
    mapping = aes(
      x = age,
      y = fantasyPoints_gam
    ),
    data = pointsPerSeason_positionAge_newData %>% filter(positionFactor == "RB"),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_gam","text","label")
)
```

#### Wide Receivers {#sec-fantasyPointsIndividualByAgePlotsWR}

A plot of Wide Receivers' model-implied fantasy points by age is in @fig-individualFantasyPointsByAgeWR.
The model-implied peaks of Wide Receivers' fantasy points are at ages `{r} wrPeakAge` and `{r} wrPeakAge2`.
The model-predicted value of zero fantasy points for Wide Receivers is at age `{r} wrZeroAge`.

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

plot_individualFantasyPointsByAgeWR <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% filter(position == "WR"),
  mapping = aes(
    x = age,
    y = fantasyPoints_gam,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_smooth(
    mapping = aes(
      x = age,
      y = fantasyPoints_gam
    ),
    data = pointsPerSeason_positionAge_newData %>% filter(positionFactor == "WR"),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_gam","text","label")
)
```

#### Tight Ends {#sec-fantasyPointsIndividualByAgePlotsTE}

A plot of Tight Ends' model-implied fantasy points by age is in @fig-individualFantasyPointsByAgeTE.
The model-implied peak of Tight Ends' fantasy points is at age `{r} tePeakAge`.
The model-predicted value of zero fantasy points for Tight Ends is at age `{r} teZeroAge`.

```{r}
#| label: fig-individualFantasyPointsByAgeTE
#| fig-cap: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, for Tight Ends. Overlaid with the Model-Implied Trajectory for Wide Tight Ends in Blue."
#| fig-alt: "Plot of Individuals' Implied Trajectories of Fantasy Points by Age, from a Generalized Additive Model, for Tight Ends. Overlaid with the Model-Implied Trajectory for Wide Tight Ends in Blue."
#| code-fold: true

plot_individualFantasyPointsByAgeTE <- ggplot(
  data = player_stats_seasonal_offense_subsetCC %>% filter(position == "TE"),
  mapping = aes(
    x = age,
    y = fantasyPoints_gam,
    group = player_id)) +
  geom_smooth(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_smooth(
    mapping = aes(
      x = age,
      y = fantasyPoints_gam
    ),
    data = pointsPerSeason_positionAge_newData %>% filter(positionFactor == "TE"),
    inherit.aes = FALSE,
    se = TRUE,
    linewidth = 2
  ) +
  geom_point(
    aes(
      x = age,
      y = fantasyPoints_gam,
      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_gam","text","label")
)
```

### Summary of Findings {#sec-fantasyPointsByAgeExperienceSummary}

We applied mixed models with random intercepts and random slopes to allow our model estimates to account for the fact that different players have different starting points (intercepts) and different changes over time (slopes) in fantasy points.
A quadratic, inverted-U-shaped form as a function of age fit better than a linear form as a function of age in predicting players' fantasy points.
A generalized additive model that allowed further nonlinearity fit even better than the quadratic model.

Based on the bivariate scatterplots between age and fantasy points, we might conclude that players tend to stay stable or even increase in fantasy points with age.
However, this conclusion would be wrong.
When we account for the longitudinal data (i.e., multiple observations over time for the same player) using mixed models, we observe that fantasy points tend to decrease with age, with the timing and rate of decline differing for each position.
In other words, the association between age and fantasy points differs at the person level versus the group level.
This is an example of [Simpson's paradox](#sec-simpsonsParadox).

The discrepancy between the positive or null association between age and players' fantasy points at the group level, and the negative association at the person level may be due, in part, to the selective attrition of players with age.
The players who play the longest will tend to be the highest performing players, whereas the poorest performing players will retire or get dropped from the team at younger ages.
Thus, the selective attrition of weaker players may make it seem that there is no association between age and performance (or even a positive one!), when in fact, players' performance tends to decrease with age after age 26 or so (with the timing differing from position to position), until the player eventually retires or is dropped from the team.
Selective attrition is common in [longitudinal studies](#sec-longitudinal) (such as this one) and [intervention studies](#sec-intervention).
For instance, attrition may be more likely for individuals from lower socioeconomic status backgrounds because they may face more challenges in continuing in longitudinal studies such as fewer financial resources, greater life stressors, etc.
In addition, people who experience side effects or lack of improvement may be more likely to drop out of [intervention studies](#sec-intervention).
Examining only those who completed treatment (an example of [selection bias](#sec-causalDiagramSelectionBias)) would make the intervention look more effective than it actually was because the people who stay in the study are those who experience the greatest improvement.
Thus, it is important to use approaches such as mixed models or other approaches that account for the multiple observations from the same person, that use all available information, and that do not exclude people who do not complete all portions of the study.

## Conclusion {#sec-mixedModelsConclusion}

Mixed models allow accounting for multiple levels or units of analysis and to include both fixed and random effects.
Inclusion of random effects allows the association between the [predictor variables](#sec-correlationalStudy) (the intercept and age) and the [outcome variable](#sec-correlationalStudy) (fantasy points) to differ for each individual in the grouping level (in this case, each player).
This allows for more accurately predicting phenomena.
Based on the bivariate scatterplots between age and fantasy points, we might conclude that players tend to stay stable or even increase in fantasy points with age.
However, this conclusion would be wrong.
When we account for the longitudinal data using mixed models, we observe that players' fantasy points tend to decrease with age, with the timing and rate of decline differing for each position.
In other words, the association between age and fantasy points differs at the person level versus the group level, which is an example of [Simpson's paradox](#sec-simpsonsParadox).
In sum, mixed models are valuable for examining associations between variables when there are multiple levels of data (i.e., multiple observations within the same unit, known as clustering or nesting).
It is important not to confuse the association at one level (e.g., group level) with the association at another level (e.g., person level), which is an example of the [ecological fallacy](#sec-fallaciesEcological).

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

## Session Info {#sec-mixedModelsSessionInfo}

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

:::