cluster-analysis.qmd

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

# Cluster Analysis {#sec-clusterAnalysis}

This chapter provides an overview of cluster analysis.

## Getting Started {#sec-clusterAnalysisGettingStarted}

### Load Packages {#sec-clusterAnalysisLoadPackages}

```{r}
library("petersenlab")
library("nflreadr")
library("mclust")
library("tidyLPA")
library("plotly")
library("tidyverse")
```

### Load Data {#sec-clusterAnalysisLoadData}

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

load(file = "./data/nfl_players.RData")
load(file = "./data/nfl_combine.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 = "./data/nfl_advancedStatsPFR_seasonal.RData")
load(file = "./data/nfl_actualStats_player_career.RData")
```

```{r}
load(file = "./data/nfl_players.RData")
load(file = "./data/nfl_combine.RData")
load(file = "./data/player_stats_weekly.RData")
load(file = "./data/player_stats_seasonal.RData")
load(file = "./data/nfl_advancedStatsPFR_seasonal.RData")
load(file = "./data/nfl_actualStats_player_career.RData")
```

We created the `player_stats_weekly.RData` and `player_stats_seasonal.RData` objects in @sec-calculatePlayerAge.

### Overview {#sec-clusterAnalysisOverview}

Whereas [factor analysis](#sec-factorAnalysis) evaluates how *variables* do or do not hang together—in terms of their associations and non-associations, cluster analysis evaluates how *people* are or or not similar—in terms of their scores on one or more variables.
The goal of cluster analysis is to identify distinguishable subgroups of people.
The people within a subgroup are expected to be more similar to each other than they are to people in other subgroups.
For instance, we might expect that there are distinguishable subtypes of Wide Receivers: possession, deep threats, and slot-type Wide Receivers.
Possession Wide Receivers tend to be taller and heavier, with good hands who catch the ball at a high rate.
Deep threat Wide Receivers tend to be fast.
Slot-type Wide Receivers tend to be small, quick, and agile.
In order to identify these clusters of Wide Receivers, we might conduct a cluster analysis with variables relating to the players' height, weight, percent of (catchable) targets caught, air yards received, and various metrics from the National Football League (NFL) Combine, including their times in the 40-yard dash, 20-yard shuttle run, and three cone drill.

There are many approaches to cluster analysis, including model-based clustering, density-based clustering, centroid-based clustering (e.g., *k*-means clustering), hierarchical clustering (aka connectivity-based clustering), etc.
In general, cluster analysis intends to maximize intracluster similarities and to minimize intercluster similarities [@Ramasubramanian2016].
Cluster analysis is a form of an [unsupervised approach to machine learning](#sec-machineLearningTypesUnsupervised) in which the groups are not known a priori.
There is a risk of reifying clusters identified in cluster analysis when, in some cases, the clusters may reflect artifacts of the analytic method rather than true distinct groups [@Everitt2011].
In general, the classification should be judged on the extent to which it is *useful* rather than on the extent to which it is true or false [@Everitt2011].

An overview of approaches to cluster analysis in `R` is provided by @Kassambara2017a.
In this chapter, we focus on examples using model-based clustering with the `R` package `mclust` [@R-mclust; @Scrucca2023_packages], which uses Gaussian finite mixture modeling.
Model-based clustering assumes the data are generated by an underlying statistical model, and the goal is to recover the structure of that model.
Model-based clustering represents the data as coming from a mixture of probability distributions.
We also use *k*-means clustering.

There are several model fit criteria we can use to identify the optimal number of clusters in terms of the tradeoff between fit and parsimony.
For the Bayesian Information Criterion (BIC), closer to zero is better.
For the Integrated Complete-data Likelihood (ICL), higher is better.
For the bootstrapped likelihood ratio tests, a significant likelihood ratio test indicates that the model with more clusters fits significantly better than the model with fewer clusters.
Although the model fit criteria can be helpful for identifying the optimal number of clusters, they should only be used as a guide.
Ultimately, cluster analysis is only useful to the extent that the clusters and useful and interpretable, which may mean selecting more (or fewer) clusters than are suggested by the model fit criteria.

For univariate cluster analyses in the `mclust` package [@R-mclust; @Scrucca2023_packages], there are two types of models: ones with equal variance across clusters (`modelName = "E"`) and ones with different variances across clusters (`modelName = "V"`).
Allowing different variances across clusters is a more complex model than forcing the clusters to have equal variances.

For multivariate cluster analyses in the `mclust` package [@R-mclust; @Scrucca2023_packages], models are specified by a combination of volume, shape, and orientation.
The model name is specified such that the first letter corresponds to the volume, the second letter corresponds to shape, and the third letter corresponds to orientation.
In terms of volume, models can be of equal volume (`E`) or variable volume (`V`) across clusters.
In terms of shape, models can be of equal (`E`), variable (`V`), or identity (`I`) shape.
In terms of orientation, models can be of equal (`E`), variable (`V`), or identity (`I`) orientation.
The varying types of volume, shape, and orientation combine to produce three general types of covariance structures: spherical (low complexity), diagonal (moderate complexity), and ellipsoidal (high complexity).
Models that allow greater differences across clusters (i.e., more `V`'s) are more flexible and typically more complex.

The various types of `mclust` models are provided here:
<https://mclust-org.github.io/mclust/reference/mclustModelNames.html>.

### Tiers of Prior Season Fantasy Points {#sec-clusterAnalysisExample}

#### Prepare Data {#sec-clusterAnalysisExamplePrepareData}

```{r}
recentSeason <- max(player_stats_seasonal$season, na.rm = TRUE) # also works: nflreadr::most_recent_season()
recentSeason

player_stats_seasonal_offense_recent <- player_stats_seasonal %>% 
  filter(season == recentSeason) %>% 
  filter(position_group %in% c("QB","RB","WR","TE"))

player_stats_seasonal_offense_recentQB <- player_stats_seasonal_offense_recent %>% 
  filter(position_group == "QB")

player_stats_seasonal_offense_recentRB <- player_stats_seasonal_offense_recent %>% 
  filter(position_group == "RB")

player_stats_seasonal_offense_recentWR <- player_stats_seasonal_offense_recent %>% 
  filter(position_group == "WR")

player_stats_seasonal_offense_recentTE <- player_stats_seasonal_offense_recent %>% 
  filter(position_group == "TE")
```

#### Identify the Optimal Number of Tiers by Position {#sec-clusterAnalysisExampleNumTiers}

We can evaluate how many clusters to keep, based on the tradeoff between fit and parsimony, using the `mclust::mclustBIC()` (BIC), `mclust::mclustICL()` (ICL), and `mclust::mclustBootstrapLRT()` (bootstrapped likelihood ratio tests) functions of the `mclust` package [@R-mclust; @Scrucca2023_packages].
We can also use *k*-means clustering using the `stats::kmeans()` function.

##### Quarterbacks {#sec-clusterAnalysisExampleNumTiersQBs}

###### Model-Based Clustering {#sec-clusterAnalysisExampleNumTiersQBsModel}

```{r}
tiersQB_bic <- mclust::mclustBIC(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
  G = 1:9
)

tiersQB_bic
summary(tiersQB_bic)
plot(tiersQB_bic)

tiersQB_icl <- mclust::mclustICL(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
  G = 1:9
)

tiersQB_icl
summary(tiersQB_icl)
plot(tiersQB_icl)

tiersQB_bootstrap <- mclust::mclustBootstrapLRT(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
  modelName = "V") # variable/unequal variance (for univariate data)

numTiersQB <- as.numeric(summary(tiersQB_bootstrap)[,"Length"][1]) # or could specify the number of teams manually

tiersQB_bootstrap
plot(
  tiersQB_bootstrap,
  G = numTiersQB - 1)
```

###### *k*-Means Clustering {#sec-clusterAnalysisExampleNumTiersQBsKmeans}

```{r}
tiersQB_kMeans <- kmeans(
  x = player_stats_seasonal_offense_recentQB$fantasyPoints,
  centers = numTiersQB)

tiersQB_kMeans
```

##### Running Backs {#sec-clusterAnalysisExampleNumTiersRBs}

###### Model-Based Clustering {#sec-clusterAnalysisExampleNumTiersRBsModel}

```{r}
tiersRB_bic <- mclust::mclustBIC(
  data = player_stats_seasonal_offense_recentRB$fantasyPoints,
  G = 1:9
)

tiersRB_bic
summary(tiersRB_bic)
plot(tiersRB_bic)

tiersRB_icl <- mclust::mclustICL(
  data = player_stats_seasonal_offense_recentRB$fantasyPoints,
  G = 1:9
)

tiersRB_icl
summary(tiersRB_icl)
plot(tiersRB_icl)

numTiersRB <- 3
```

<!--
The model-based bootstrap clustering of Running Backs' fantasy points is unable to run due to an error:
UPDATE: THE ERROR NO LONGER OCCURS
-->

```{r}
tiersRB_bootstrap <- mclust::mclustBootstrapLRT(
  data = player_stats_seasonal_offense_recentRB$fantasyPoints,
  modelName = "V") # variable/unequal variance (for univariate data)
```

<!--
Thus, we cannot use the following code, which would otherwise summarize the model results, specify the number of tiers, and plot model comparisons:
UPDATE: THE ERROR NO LONGER OCCURS
-->

```{r}
numTiersRB <- as.numeric(summary(tiersRB_bootstrap)[,"Length"][1]) # or could specify the number of teams manually

tiersRB_bootstrap
plot(
  tiersRB_bootstrap,
  G = numTiersRB - 1)
```

###### *k*-Means Clustering {#sec-clusterAnalysisExampleNumTiersRBsKmeans}

```{r}
tiersRB_kMeans <- kmeans(
  x = player_stats_seasonal_offense_recentRB$fantasyPoints,
  centers = numTiersRB)

tiersRB_kMeans
```

##### Wide Receivers {#sec-clusterAnalysisExampleNumTiersWRs}

###### Model-Based Clustering {#sec-clusterAnalysisExampleNumTiersWRsModel}

```{r}
tiersWR_bic <- mclust::mclustBIC(
  data = player_stats_seasonal_offense_recentWR$fantasyPoints,
  G = 1:9
)

tiersWR_bic
summary(tiersWR_bic)
plot(tiersWR_bic)

tiersWR_icl <- mclust::mclustICL(
  data = player_stats_seasonal_offense_recentWR$fantasyPoints,
  G = 1:9
)

tiersWR_icl
summary(tiersWR_icl)
plot(tiersWR_icl)

tiersWR_bootstrap <- mclust::mclustBootstrapLRT(
  data = player_stats_seasonal_offense_recentWR$fantasyPoints,
  modelName = "V") # variable/unequal variance (for univariate data)

numTiersWR <- as.numeric(summary(tiersWR_bootstrap)[,"Length"][1]) # or could specify the number of teams manually

tiersWR_bootstrap
plot(
  tiersWR_bootstrap,
  G = numTiersWR - 1)
```

###### *k*-Means Clustering {#sec-clusterAnalysisExampleNumTiersWRsKmeans}

```{r}
tiersWR_kMeans <- kmeans(
  x = player_stats_seasonal_offense_recentWR$fantasyPoints,
  centers = numTiersWR)

tiersWR_kMeans
```

##### Tight Ends {#sec-clusterAnalysisExampleNumTiersTEs}

###### Model-Based Clustering {#sec-clusterAnalysisExampleNumTiersTEsModel}

```{r}
tiersTE_bic <- mclust::mclustBIC(
  data = player_stats_seasonal_offense_recentTE$fantasyPoints,
  G = 1:9
)

tiersTE_bic
summary(tiersTE_bic)
plot(tiersTE_bic)

tiersTE_icl <- mclust::mclustICL(
  data = player_stats_seasonal_offense_recentTE$fantasyPoints,
  G = 1:9
)

tiersTE_icl
summary(tiersTE_icl)
plot(tiersTE_icl)

tiersTE_bootstrap <- mclust::mclustBootstrapLRT(
  data = player_stats_seasonal_offense_recentTE$fantasyPoints,
  modelName = "V") # variable/unequal variance (for univariate data)

numTiersTE <- as.numeric(summary(tiersTE_bootstrap)[,"Length"][1]) # or could specify the number of teams manually

tiersTE_bootstrap
plot(
  tiersTE_bootstrap,
  G = numTiersTE - 1)
```

###### *k*-Means Clustering {#sec-clusterAnalysisExampleNumTiersTEsKmeans}

```{r}
tiersTE_kMeans <- kmeans(
  x = player_stats_seasonal_offense_recentTE$fantasyPoints,
  centers = numTiersTE)

tiersTE_kMeans
```

#### Fit the Cluster Model to the Optimal Number of Tiers {#sec-clusterAnalysisExampleModel}

##### Quarterbacks {#sec-clusterAnalysisExampleModelQBs}

In our data, all of the following models are equivalent—i.e., they result in the same unequal variance model with a 4-cluster solution—but they arrive there in different ways.
We can fit the cluster model to the optimal number of tiers using the `mclust::Mclust()` function.

```{r}
#| eval: false

mclust::Mclust(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
  G = numTiersQB,
)

mclust::Mclust(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
  G = 4,
)

mclust::Mclust(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
)

mclust::Mclust(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
  x = tiersQB_bic
)
```

Let's fit one of these:

```{r}
clusterModelQBs <- mclust::Mclust(
  data = player_stats_seasonal_offense_recentQB$fantasyPoints,
  G = numTiersQB,
)
```

Here are the number of players that are in each of the four clusters (i.e., tiers):

```{r}
table(clusterModelQBs$classification)
```

##### Running Backs {#sec-clusterAnalysisExampleModelRBs}

```{r}
clusterModelRBs <- mclust::Mclust(
  data = player_stats_seasonal_offense_recentRB$fantasyPoints,
  G = numTiersRB,
)
```

Here are the number of players that are in each of the four clusters (i.e., tiers):

```{r}
table(clusterModelRBs$classification)
```

##### Wide Receivers {#sec-clusterAnalysisExampleModelWRs}

```{r}
clusterModelWRs <- mclust::Mclust(
  data = player_stats_seasonal_offense_recentWR$fantasyPoints,
  G = numTiersWR,
)
```

Here are the number of players that are in each of the four clusters (i.e., tiers):

```{r}
table(clusterModelWRs$classification)
```

##### Tight Ends {#sec-clusterAnalysisExampleModelTEs}

```{r}
clusterModelTEs <- mclust::Mclust(
  data = player_stats_seasonal_offense_recentTE$fantasyPoints,
  G = numTiersTE,
)
```

Here are the number of players that are in each of the four clusters (i.e., tiers):

```{r}
table(clusterModelTEs$classification)
```

#### Plot the Tiers {#sec-clusterAnalysisExamplePlotTiers}

We can merge the player's classification into the dataset and plot each player's classification.

##### Quarterbacks {#sec-clusterAnalysisExamplePlotTiersQB}

```{r}
#| label: fig-qbTiers
#| fig-cap: "Quarterback Fantasy Points by Tier."
#| fig-alt: "Quarterback Fantasy Points by Tier."

player_stats_seasonal_offense_recentQB$tier <- clusterModelQBs$classification

player_stats_seasonal_offense_recentQB <- player_stats_seasonal_offense_recentQB %>%
  mutate(
    tier = factor(max(tier, na.rm = TRUE) + 1 - tier)
  )

player_stats_seasonal_offense_recentQB$position_rank <- rank(
  player_stats_seasonal_offense_recentQB$fantasyPoints * -1,
  na.last = "keep",
  ties.method = "min")

plot_qbTiers <- ggplot2::ggplot(
  data = player_stats_seasonal_offense_recentQB,
  mapping = aes(
    x = fantasyPoints,
    y = position_rank,
    color = tier
  )) +
  geom_point(
    aes(
      text = player_display_name # add player name for mouse over tooltip
  )) +
  scale_y_continuous(trans = "reverse") +
  coord_cartesian(clip = "off") +
  labs(
    x = "Projected Points",
    y = "Position Rank",
    title = "Quarterback Fantasy Points by Tier",
    color = "Tier") +
  theme_classic() +
  theme(legend.position = "top")

plotly::ggplotly(plot_qbTiers)
```

##### Running Backs {#sec-clusterAnalysisExamplePlotTiersRBs}

```{r}
#| label: fig-rbTiers
#| fig-cap: "Running Back Fantasy Points by Tier."
#| fig-alt: "Running Back Fantasy Points by Tier."

player_stats_seasonal_offense_recentRB$tier <- clusterModelRBs$classification

player_stats_seasonal_offense_recentRB <- player_stats_seasonal_offense_recentRB %>%
  mutate(
    tier = factor(max(tier, na.rm = TRUE) + 1 - tier)
  )

player_stats_seasonal_offense_recentRB$position_rank <- rank(
  player_stats_seasonal_offense_recentRB$fantasyPoints * -1,
  na.last = "keep",
  ties.method = "min")

plot_rbTiers <- ggplot2::ggplot(
  data = player_stats_seasonal_offense_recentRB,
  mapping = aes(
    x = fantasyPoints,
    y = position_rank,
    color = tier
  )) +
  geom_point(
    aes(
      text = player_display_name # add player name for mouse over tooltip
  )) +
  scale_y_continuous(trans = "reverse") +
  coord_cartesian(clip = "off") +
  labs(
    x = "Projected Points",
    y = "Position Rank",
    title = "Running Back Fantasy Points by Tier",
    color = "Tier") +
  theme_classic() +
  theme(legend.position = "top")

plotly::ggplotly(plot_rbTiers)
```

##### Wide Receivers {#sec-clusterAnalysisExamplePlotTiersWRs}

```{r}
#| label: fig-wrTiers
#| fig-cap: "Quarterback Fantasy Points by Tier."
#| fig-alt: "Quarterback Fantasy Points by Tier."

player_stats_seasonal_offense_recentWR$tier <- clusterModelWRs$classification

player_stats_seasonal_offense_recentWR <- player_stats_seasonal_offense_recentWR %>%
  mutate(
    tier = factor(max(tier, na.rm = TRUE) + 1 - tier)
  )

player_stats_seasonal_offense_recentWR$position_rank <- rank(
  player_stats_seasonal_offense_recentWR$fantasyPoints * -1,
  na.last = "keep",
  ties.method = "min")

plot_wrTiers <- ggplot2::ggplot(
  data = player_stats_seasonal_offense_recentWR,
  mapping = aes(
    x = fantasyPoints,
    y = position_rank,
    color = tier
  )) +
  geom_point(
    aes(
      text = player_display_name # add player name for mouse over tooltip
  )) +
  scale_y_continuous(trans = "reverse") +
  coord_cartesian(clip = "off") +
  labs(
    x = "Projected Points",
    y = "Position Rank",
    title = "Wide Receiver Fantasy Points by Tier",
    color = "Tier") +
  theme_classic() +
  theme(legend.position = "top")

plotly::ggplotly(plot_wrTiers)
```

##### Tight Ends {#sec-clusterAnalysisExamplePlotTiersTEs}

```{r}
#| label: fig-teTiers
#| fig-cap: "Tight End Fantasy Points by Tier."
#| fig-alt: "Tight End Fantasy Points by Tier."

player_stats_seasonal_offense_recentTE$tier <- clusterModelTEs$classification

player_stats_seasonal_offense_recentTE <- player_stats_seasonal_offense_recentTE %>%
  mutate(
    tier = factor(max(tier, na.rm = TRUE) + 1 - tier)
  )

player_stats_seasonal_offense_recentTE$position_rank <- rank(
  player_stats_seasonal_offense_recentTE$fantasyPoints * -1,
  na.last = "keep",
  ties.method = "min")

plot_teTiers <- ggplot2::ggplot(
  data = player_stats_seasonal_offense_recentTE,
  mapping = aes(
    x = fantasyPoints,
    y = position_rank,
    color = tier
  )) +
  geom_point(
    aes(
      text = player_display_name # add player name for mouse over tooltip
  )) +
  scale_y_continuous(trans = "reverse") +
  coord_cartesian(clip = "off") +
  labs(
    x = "Projected Points",
    y = "Position Rank",
    title = "Tight End Fantasy Points by Tier",
    color = "Tier") +
  theme_classic() +
  theme(legend.position = "top")

plotly::ggplotly(plot_teTiers)
```

### Types of Wide Receivers {#sec-clusterAnalysisWRtypes}

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

# names(nfl_players) #gsis_id: height, weight
# names(nfl_combine) #gsis_id: vertical, forty, ht, wt
# names(player_stats_seasonal_offense) #player_id, season: receptions, receiving_air_yards, air_yards_share, target_share
# names(nfl_advancedStatsPFR_seasonal) #gsis_id, season: adot.rec, rec.rec, brk_tkl.rec, drop.rec, drop_percent.rec
# names(nfl_actualStats_career_player_inclPost) #player_id: receptions, targets, receiving_air_yards, air_yards_share, target_share
```

```{r}
# Compute Advanced PFR Stats by Career
pfrVars <- nfl_advancedStatsPFR_seasonal %>% 
  select(pocket_time.pass:cmp_percent.def, g, gs) %>% 
  names()

weightedAverageVars <- c(
  "pocket_time.pass",
  "ybc_att.rush","yac_att.rush",
  "ybc_r.rec","yac_r.rec","adot.rec","rat.rec",
  "yds_cmp.def","yds_tgt.def","dadot.def","m_tkl_percent.def","rat.def"
)

recomputeVars <- c(
  "drop_pct.pass", # drops.pass / pass_attempts.pass
  "bad_throw_pct.pass", # bad_throws.pass / pass_attempts.pass
  "on_tgt_pct.pass", # on_tgt_throws.pass / pass_attempts.pass
  "pressure_pct.pass", # times_pressured.pass / pass_attempts.pass
  "drop_percent.rec", # drop.rec / tgt.rec
  "rec_br.rec", # rec.rec / brk_tkl.rec
  "cmp_percent.def" # cmp.def / tgt.def
)

sumVars <- pfrVars[pfrVars %ni% c(
  weightedAverageVars, recomputeVars,
  "merge_name", "loaded.pass", "loaded.rush", "loaded.rec", "loaded.def")]

nfl_advancedStatsPFR_career <- nfl_advancedStatsPFR_seasonal %>% 
  group_by(pfr_id, merge_name) %>% 
  summarise(
    across(all_of(weightedAverageVars), ~ weighted.mean(.x, w = g, na.rm = TRUE)),
    across(all_of(sumVars), ~ sum(.x, na.rm = TRUE)),
    .groups = "drop") %>% 
  mutate(
    drop_pct.pass = drops.pass / pass_attempts.pass,
    bad_throw_pct.pass = bad_throws.pass / pass_attempts.pass,
    on_tgt_pct.pass = on_tgt_throws.pass / pass_attempts.pass,
    pressure_pct.pass = times_pressured.pass / pass_attempts.pass,
    drop_percent.rec = drop.rec / tgt.rec,
    rec_br.rec = drop.rec / tgt.rec,
    cmp_percent.def = cmp.def / tgt.def
  )

uniqueCases <- nfl_advancedStatsPFR_seasonal %>% select(pfr_id, merge_name, gsis_id) %>% unique()

uniqueCases %>%
  group_by(pfr_id) %>% 
  filter(n() > 1)

nfl_advancedStatsPFR_seasonal <- nfl_advancedStatsPFR_seasonal %>% 
  filter(pfr_id != "WillMa06" | merge_name != "MARCUSWILLIAMS" | !is.na(gsis_id))


nfl_advancedStatsPFR_career <- left_join(
  nfl_advancedStatsPFR_career,
  nfl_advancedStatsPFR_seasonal %>% select(pfr_id, merge_name, gsis_id) %>% unique(),
  by = c("pfr_id", "merge_name")
)

# Compute Player Stats Per Season
player_stats_seasonal_careerWRs <- player_stats_seasonal %>% 
  filter(position == "WR") %>% 
  group_by(player_id) %>% 
  summarise(
    across(all_of(c("targets", "receptions", "receiving_air_yards")), ~ weighted.mean(.x, w = games, na.rm = TRUE)),
    .groups = "drop")

# Drop players with no receiving air yards
player_stats_seasonal_careerWRs <- player_stats_seasonal_careerWRs %>% 
  filter(receiving_air_yards != 0) %>% 
  rename(
    targets_per_season = targets,
    receptions_per_season = receptions,
    receiving_air_yards_per_season = receiving_air_yards
  )

# Merge
playerListToMerge <- list(
  nfl_players %>% select(gsis_id, display_name, position, height, weight),
  nfl_combine %>% select(gsis_id, vertical, forty, ht, wt),
  player_stats_seasonal_careerWRs %>% select(player_id, targets_per_season, receptions_per_season, receiving_air_yards_per_season) %>% 
    rename(gsis_id = player_id),
  nfl_actualStats_career_player_inclPost %>% select(player_id, receptions, targets, receiving_air_yards, air_yards_share, target_share) %>% 
    rename(gsis_id = player_id),
  nfl_advancedStatsPFR_career %>% select(gsis_id, adot.rec, rec.rec, brk_tkl.rec, drop.rec, drop_percent.rec)
)

merged_data <- playerListToMerge %>% 
  reduce(
    full_join,
    by = c("gsis_id"),
    na_matches = "never")
```

Additional processing:

```{r}
merged_data <- merged_data %>% 
  mutate(
    height_coalesced = coalesce(height, ht),
    weight_coalesced = coalesce(weight, wt),
    receptions_coalesced = pmax(receptions, rec.rec, na.rm = TRUE),
    receiving_air_yards_per_rec = receiving_air_yards / receptions
  )

merged_data$receiving_air_yards_per_rec[which(merged_data$receptions == 0)] <- 0

merged_dataWRs <- merged_data %>% 
  filter(position == "WR")

merged_dataWRs_cluster <- merged_dataWRs %>% 
  filter(receptions_coalesced >= 100) %>% # keep WRs with at least 100 receptions
  select(gsis_id, display_name, vertical, forty, height_coalesced, weight_coalesced, adot.rec, drop_percent.rec, receiving_air_yards_per_rec, brk_tkl.rec, receptions_per_season) %>% #targets_per_season, receiving_air_yards_per_season, air_yards_share, target_share
  na.omit()
```

#### Identify the Number of WR Types {#sec-clusterAnalysisNumWRtypes}

##### Model-Based Clustering  {#sec-clusterAnalysisNumWRtypesModel}

```{r}
wrTypes_bic <- mclust::mclustBIC(
  data = merged_dataWRs_cluster %>% select(-gsis_id, -display_name),
  G = 1:9
)

wrTypes_bic
summary(wrTypes_bic)
plot(wrTypes_bic)

wrTypes_icl <- mclust::mclustICL(
  data = merged_dataWRs_cluster %>% select(-gsis_id, -display_name),
  G = 1:9
)

wrTypes_icl
summary(wrTypes_icl)
plot(wrTypes_icl)
```

Based on the cluster analyses, it appears that three clusters are the best fit to the data.

```{r}
numTypesWR <- 3
```

```{r}
#| eval: false

wrTypes_bootstrap <- mclust::mclustBootstrapLRT(
  data = merged_dataWRs_cluster %>% select(-gsis_id, -display_name),
  modelName = "EVE") # ellipsoidal with equal volume, variable shape, and equal orientation (for multivariate data)

wrTypes_bootstrap
plot(
  wrTypes_bootstrap,
  G = numTypesWR - 1)
```

We can also use the `tidyLPA` package [@R-tidyLPA; @Rosenberg2018_packages], which provides an interface to the `mclust` package [@R-mclust; @Scrucca2023_packages].

```{r}
# model 1 (EEI): Equal variances and covariances fixed to 0
# model 2 (VVI): Varying variances and covariances fixed to 0
# model 3 (EEE): Equal variances and equal covariances
# model 4      : Varying variances and equal covariances (not able to be fit w/ mclust but can be fit using package = "mplus", if installed)
# model 5      : Equal variances and varying covariances (not able to be fit w/ mclust but can be fit using package = "mplus", if installed)
# model 6 (VVV): Varying variances and varying covariances

wrTypes_lpa_classes <- tidyLPA::estimate_profiles(
  df = merged_dataWRs_cluster %>% select(-gsis_id, -display_name),
  n_profiles = 1:6,
  models = c(1, 2, 3) #c(1, 2, 3, 6) takes too long to run because of the model with varying variances and varying covariances (model 6)
)

wrTypes_lpa_classes

tidyLPA::compare_solutions(
  wrTypes_lpa_classes,
  statistics = c("AIC", "BIC", "Entropy", "BLRT_p"))

wrTypes_lpa <- tidyLPA::estimate_profiles(
  df = merged_dataWRs_cluster %>% select(-gsis_id, -display_name),
  n_profiles = numTypesWR,
  model = 3 # equal variances and equal covariances
)

wrTypes_lpa
tidyLPA::get_fit(wrTypes_lpa)
```

##### *k*-Means Clustering  {#sec-clusterAnalysisNumWRtypesKmeans}

```{r}
wrTypes_kMeans <- kmeans(
  x = merged_dataWRs_cluster %>% select(-gsis_id, -display_name),
  centers = numTypesWR)

wrTypes_kMeans
```

#### Fit the Cluster Model to the Optimal Number of WR Types {#sec-clusterAnalysisModelWRtypes}

```{r}
clusterModelWRtypes <- mclust::Mclust(
  data = merged_dataWRs_cluster %>% select(-gsis_id, -display_name),
  G = numTypesWR,
)

summary(clusterModelWRtypes)
```

#### Plots of the Cluster Model {#sec-clusterAnalysisPlotsWRtypes}

```{r}
plot(
  clusterModelWRtypes,
  what = "BIC")
```

```{r}
plot(
  clusterModelWRtypes,
  what = "classification")
```

```{r}
plot(
  clusterModelWRtypes,
  what = "uncertainty")
```

```{r}
plot(
  clusterModelWRtypes,
  what = "density")
```

```{r}
tidyLPA::plot_profiles(wrTypes_lpa)
```

#### Interpreting the Clusters {#sec-clusterAnalysisInterpretationWRtypes}

```{r}
table(clusterModelWRtypes$classification)

merged_dataWRs_cluster$type <- clusterModelWRtypes$classification

merged_dataWRs_cluster %>% 
  group_by(type) %>% 
  summarise(across(
    where(is.numeric),
    ~ mean(., na.rm = TRUE)
    )) %>% 
  t() %>% 
  round(., 2)
```

Based on this analysis (and the variables included), there appear to be three types of Wide Receivers.
We examined the following variables: the player's vertical jump in the NFL Combine,40-yard-dash time in the NFL Combine, height, weight, average depth of target, drop percentage, receiving air yards per reception, broken tackles, and receptions per season.

Type 1 Wide Receivers included the Elite WR1s who are strong possession receivers (note: not all players in a given cluster map on perfectly to the typology—i.e., not all Type 1 Wide Receivers are elite WR1s).
They tended to have the lowest drop percentage, the shortest average depth of target, and the fewest receiving air yards per reception.
They tended to have the most receptions per season and break the most tackles.

Type 2 Wide Receivers included the consistent contributor, WR2 types.
They had fewer receptions and fewer broken tackles than Type 1 Wide Receivers.
Their average depth of target was longer than Type 1, and they had more receiving air yards per reception than Type 1.

Type 3 Wide Receivers included the deep threats.
They had the greatest average depth of target and the most receiving yards per reception.
However, they also had the fewest receptions, the highest drop percentage, and the fewest broken tackles.
Thus, they may be considered the boom-or-bust Wide Receivers.

The tiers were not particularly distinguishable based on their height, weight, vertical jump, or forty-yard dash time.

Type 1 ("Elite/WR1") WRs:

```{r}
merged_dataWRs_cluster %>% 
  filter(type == 1) %>% 
  select(display_name)
```

Type 2 ("Consistent Contributor/WR2") WRs:

```{r}
merged_dataWRs_cluster %>% 
  filter(type == 2) %>% 
  select(display_name)
```

Type 3 ("Deep Threat/Boom-or-Bust") WRs:

```{r}
merged_dataWRs_cluster %>% 
  filter(type == 3) %>% 
  select(display_name)
```

And here are the clusters based on the `tidyLPA` [@R-tidyLPA; @Rosenberg2018_packages] solution:

```{r}
merged_dataWRs_cluster <- cbind(
  merged_dataWRs_cluster,
  tidyLPA::get_data(wrTypes_lpa) %>% 
    select(model_number, classes_number, CPROB1:Class)
)

merged_dataWRs_cluster %>% 
  select(display_name, CPROB1:Class)

merged_dataWRs_cluster %>% 
  group_by(Class) %>% 
  summarise(across(
    where(is.numeric),
    ~ mean(., na.rm = TRUE)
    )) %>% 
  t() %>% 
  round(., 2)
```

## Conclusion {#sec-clusterAnalysisConclusion}

The goal of cluster analysis is to identify distinguishable subgroups of people.
There are many approaches to cluster analysis, including model-based clustering, density-based clustering, centroid-based clustering, hierarchical clustering (aka connectivity-based clustering), and others.
The present chapter used model-based clustering to identify tiers of players based on projected points.
Using various performance metrics of Wide Receivers, we identified three subtypes of Wide Receivers: 1) Elite WR1s who are strong possession receivers; 2) Consistent Contributor/WR2s; 3) deep threats/boom-or-bust receivers.
The "Elite WR1s" tended to have the lowest drop percentage, the shortest average depth of target, the fewest receiving air yards per reception, the most receptions per season, and the most broken tackles.
The "Consistent Contributor/WR2s" had fewer receptions and fewer broken tackles than the Elite WR1s; their average depth of target was longer than Elite WR1s, and they had more receiving air yards per reception than Elite WR1s.
The "Deep Threat/Boom-or-Bust" receivers had the greatest average depth of target and the most receiving yards per reception; however, they also had the fewest receptions, the highest drop percentage, and the fewest broken tackles.
In sum, cluster analysis can be a useful way of identifying subgroups of individuals who are more similar to one another on various characteristics.

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

## Session Info {#sec-clusterAnalysisSessionInfo}

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

:::