I want your feedback to make the book better for you and other readers. If you find typos, errors, or places where the text may be improved, please let me know. The best ways to provide feedback are by GitHub or hypothes.is annotations.
You can leave a comment at the bottom of the page/chapter, or open an issue or submit a pull request on GitHub: https://github.com/isaactpetersen/Fantasy-Football-Analytics-Textbook
Alternatively, you can leave an annotation using hypothes.is.
To add an annotation, select some text and then click the
symbol on the pop-up menu.
To see the annotations of others, click the
symbol in the upper right-hand corner of the page.
21 Cluster Analysis
This chapter provides an overview of cluster analysis.
21.1 Getting Started
21.1.1 Load Packages
21.1.2 Load Data
We created the player_stats_weekly.RData and player_stats_seasonal.RData objects in Section 4.4.3.
21.1.3 Overview
Whereas factor analysis 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 (Ramasubramanian & Singh, 2016). Cluster analysis is a form of an unsupervised approach to machine learning 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 (Everitt et al., 2011). 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 (Everitt et al., 2011).
An overview of approaches to cluster analysis in R is provided by Kassambara (2017). In this chapter, we focus on examples using model-based clustering with the R package mclust (Fraley et al., 2024; Scrucca et al., 2023), 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 (Fraley et al., 2024; Scrucca et al., 2023), 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 (Fraley et al., 2024; Scrucca et al., 2023), 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.
21.1.4 Tiers of Prior Season Fantasy Points
21.1.4.1 Prepare Data
Code
[1] 2025Code
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")21.1.4.2 Identify the Optimal Number of Tiers by Position
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 (Fraley et al., 2024; Scrucca et al., 2023). We can also use k-means clustering using the stats::kmeans() function.
21.1.4.2.1 Quarterbacks
21.1.4.2.1.1 Model-Based Clustering
Code
Bayesian Information Criterion (BIC):
E V
1 -1011.1047 -1011.1047
2 -986.5400 -947.7960
3 -979.7968 -924.9390
4 -988.5831 -933.2050
5 -979.4763 -943.8205
6 -983.3936 -943.8137
7 -992.1762 -955.4531
8 -999.6852 -951.8021
9 -1008.4518 -964.2137
Top 3 models based on the BIC criterion:
V,3 V,4 V,6
-924.9390 -933.2050 -943.8137 Best BIC values:
V,3 V,4 V,6
BIC -924.939 -933.20505 -943.81375
BIC diff 0.000 -8.26608 -18.87478
Code
Integrated Complete-data Likelihood (ICL) criterion:
E V
1 -1011.1047 -1011.1047
2 -992.3740 -953.2231
3 -989.9729 -935.6124
4 -1056.4678 -948.2056
5 -1042.1246 -970.0466
6 -1043.7012 -953.0216
7 -1080.9527 -971.2156
8 -1089.2292 -960.5491
9 -1107.1317 -976.4461
Top 3 models based on the ICL criterion:
V,3 V,4 V,6
-935.6124 -948.2056 -953.0216 Best ICL values:
V,3 V,4 V,6
ICL -935.6124 -948.2056 -953.02162
ICL diff 0.0000 -12.5932 -17.40924
Code
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-------------------------------------------------------------
Bootstrap sequential LRT for the number of mixture components
-------------------------------------------------------------
Model = V
Replications = 999
LRTS bootstrap p-value
1 vs 2 76.492030 0.001
2 vs 3 36.040421 0.001
3 vs 4 4.917267 0.237
21.1.4.2.1.2 k-Means Clustering
Code
K-means clustering with 3 clusters of sizes 20, 44, 17
Cluster means:
[,1]
1 291.87700
2 21.15818
3 147.02824
Clustering vector:
[1] 1 2 2 2 2 1 1 2 2 2 3 1 3 1 3 2 2 2 2 1 1 2 2 1 2 2 3 3 1 2 1 2 2 1 2 1 3 2
[39] 3 3 2 1 1 2 2 3 1 2 2 3 2 2 1 3 2 3 2 1 2 3 2 2 1 2 2 2 2 1 2 3 3 2 2 1 2 3
[77] 2 3 2 2 2
Within cluster sum of squares by cluster:
[1] 47123.00 25503.57 19020.71
(between_SS / total_SS = 91.8 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" 21.1.4.2.2 Running Backs
21.1.4.2.2.1 Model-Based Clustering
Code
Bayesian Information Criterion (BIC):
E V
1 -1824.214 -1824.214
2 -1768.968 -1674.774
3 -1779.083 -1666.523
4 -1789.222 -1664.721
5 -1761.238 -1652.250
6 -1771.348 -1663.521
7 -1744.171 -1678.616
8 -1754.287 -1687.786
9 -1751.696 -1687.019
Top 3 models based on the BIC criterion:
V,5 V,6 V,4
-1652.250 -1663.521 -1664.721 Best BIC values:
V,5 V,6 V,4
BIC -1652.25 -1663.52140 -1664.72118
BIC diff 0.00 -11.27187 -12.47165
Code
Integrated Complete-data Likelihood (ICL) criterion:
E V
1 -1824.214 -1824.214
2 -1778.797 -1692.619
3 -1948.646 -1708.755
4 -2032.017 -1715.950
5 -2013.956 -1697.932
6 -2064.069 -1718.785
7 -2014.630 -1761.661
8 -2050.462 -1748.965
9 -2037.163 -1739.977
Top 3 models based on the ICL criterion:
V,2 V,5 V,3
-1692.619 -1697.932 -1708.755 Best ICL values:
V,2 V,5 V,3
ICL -1692.619 -1697.932128 -1708.75475
ICL diff 0.000 -5.312959 -16.13558
Code
-------------------------------------------------------------
Bootstrap sequential LRT for the number of mixture components
-------------------------------------------------------------
Model = V
Replications = 999
LRTS bootstrap p-value
1 vs 2 164.609086 0.001
2 vs 3 23.419296 0.001
3 vs 4 16.971055 0.004
4 vs 5 27.640391 0.001
5 vs 6 3.896865 0.325
21.1.4.2.2.2 k-Means Clustering
Code
K-means clustering with 5 clusters of sizes 18, 20, 86, 7, 26
Cluster means:
[,1]
1 186.95889
2 115.68600
3 11.71372
4 292.65714
5 60.02692
Clustering vector:
[1] 3 2 3 3 5 3 3 3 3 1 5 3 3 5 5 4 2 3 3 5 1 5 3 3 2 3 3 2 1 3 3 2 4 2 3 3 3
[38] 3 1 3 3 3 3 2 4 3 4 5 2 3 5 3 5 5 3 3 3 3 3 3 3 5 5 2 3 3 2 4 3 3 3 4 3 1
[75] 3 3 1 5 5 3 4 2 3 1 3 3 5 3 3 2 5 3 3 3 1 1 3 2 3 2 1 3 3 5 5 3 3 3 3 5 3
[112] 5 2 3 3 3 2 1 2 3 3 3 5 3 1 1 3 3 3 5 1 5 3 3 5 3 1 1 3 3 1 3 3 3 3 5 5 2
[149] 3 3 2 3 3 2 1 3 3
Within cluster sum of squares by cluster:
[1] 10035.805 5124.514 8032.044 3455.057 5195.611
(between_SS / total_SS = 96.7 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" 21.1.4.2.3 Wide Receivers
21.1.4.2.3.1 Model-Based Clustering
Code
Bayesian Information Criterion (BIC):
E V
1 -2588.178 -2588.178
2 -2511.626 -2448.005
3 -2522.606 -2391.282
4 -2533.593 -2393.442
5 -2506.194 -2401.111
6 -2517.163 -2412.029
7 -2496.389 -2428.269
8 -2507.377 NA
9 -2507.671 NA
Top 3 models based on the BIC criterion:
V,3 V,4 V,5
-2391.282 -2393.442 -2401.111 Best BIC values:
V,3 V,4 V,5
BIC -2391.282 -2393.441648 -2401.110804
BIC diff 0.000 -2.159824 -9.828981
Code
Integrated Complete-data Likelihood (ICL) criterion:
E V
1 -2588.178 -2588.178
2 -2526.398 -2492.180
3 -2773.174 -2459.017
4 -2898.343 -2491.014
5 -2902.060 -2522.908
6 -2985.838 -2514.893
7 -2927.266 -2608.115
8 -2987.186 NA
9 -2945.822 NA
Top 3 models based on the ICL criterion:
V,3 V,4 V,2
-2459.017 -2491.014 -2492.180 Best ICL values:
V,3 V,4 V,2
ICL -2459.017 -2491.01428 -2492.18041
ICL diff 0.000 -31.99739 -33.16352
Code
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-------------------------------------------------------------
Bootstrap sequential LRT for the number of mixture components
-------------------------------------------------------------
Model = V
Replications = 999
LRTS bootstrap p-value
1 vs 2 156.627525 0.001
2 vs 3 73.177288 0.001
3 vs 4 14.294566 0.004
4 vs 5 8.785234 0.050
5 vs 6 5.536546 0.177
21.1.4.2.3.2 k-Means Clustering
Code
K-means clustering with 5 clusters of sizes 114, 5, 58, 35, 29
Cluster means:
[,1]
1 8.60614
2 220.78000
3 40.20379
4 81.86286
5 135.87586
Clustering vector:
[1] 5 1 3 5 1 1 2 1 3 1 1 1 1 1 3 1 3 1 4 1 1 3 3 3 1 3 1 5 1 4 1 3 3 5 3 4 1
[38] 1 4 5 1 5 1 5 1 1 1 1 1 4 3 5 1 3 1 3 4 5 1 5 3 1 3 1 1 3 1 1 3 5 3 1 1 4
[75] 1 5 1 1 4 2 3 3 1 4 3 1 3 2 3 1 3 1 1 4 1 1 3 1 1 4 1 3 1 1 1 1 5 1 1 5 2
[112] 4 3 5 3 3 4 1 1 3 1 4 1 1 4 3 1 3 5 1 4 1 3 1 4 1 1 4 3 4 1 1 4 1 3 1 5 1
[149] 1 3 4 4 1 1 1 3 1 4 1 4 4 3 1 1 3 1 1 5 5 3 1 1 1 5 1 1 3 5 1 3 2 5 1 4 4
[186] 3 1 3 1 1 3 4 5 1 4 1 1 1 1 1 1 1 5 3 1 1 1 1 5 4 5 4 1 1 3 1 3 3 4 1 1 5
[223] 3 1 3 3 4 3 3 3 1 5 1 4 3 1 3 1 4 5 1
Within cluster sum of squares by cluster:
[1] 6115.286 2411.788 5731.184 7061.902 7416.050
(between_SS / total_SS = 95.4 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" 21.1.4.2.4 Tight Ends
21.1.4.2.4.1 Model-Based Clustering
Code
Bayesian Information Criterion (BIC):
E V
1 -1393.191 -1393.191
2 -1351.232 -1302.337
3 -1361.086 -1264.700
4 -1358.460 -1263.340
5 -1368.291 -1270.694
6 -1368.182 -1288.864
7 -1378.019 -1289.230
8 -1387.858 -1295.137
9 -1397.710 -1307.786
Top 3 models based on the BIC criterion:
V,4 V,3 V,5
-1263.340 -1264.700 -1270.694 Best BIC values:
V,4 V,3 V,5
BIC -1263.34 -1264.700320 -1270.694175
BIC diff 0.00 -1.359878 -7.353732
Code
Integrated Complete-data Likelihood (ICL) criterion:
E V
1 -1393.191 -1393.191
2 -1360.937 -1325.426
3 -1505.560 -1293.634
4 -1494.452 -1300.262
5 -1577.023 -1308.663
6 -1558.473 -1345.081
7 -1594.663 -1314.113
8 -1620.635 -1315.866
9 -1673.929 -1333.826
Top 3 models based on the ICL criterion:
V,3 V,4 V,5
-1293.634 -1300.262 -1308.663 Best ICL values:
V,3 V,4 V,5
ICL -1293.634 -1300.262248 -1308.66264
ICL diff 0.000 -6.628154 -15.02854
Code
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-------------------------------------------------------------
Bootstrap sequential LRT for the number of mixture components
-------------------------------------------------------------
Model = V
Replications = 999
LRTS bootstrap p-value
1 vs 2 105.613956 0.001
2 vs 3 52.396927 0.001
3 vs 4 16.119821 0.005
4 vs 5 7.406211 0.111
21.1.4.2.4.2 k-Means Clustering
Code
K-means clustering with 4 clusters of sizes 20, 60, 26, 31
Cluster means:
[,1]
1 108.715000
2 4.606667
3 22.400000
4 52.147742
Clustering vector:
[1] 1 3 2 3 2 4 2 3 2 2 2 2 4 3 1 3 4 2 2 2 2 3 2 4 2 1 4 2 1 1 1 1 4 4 4 3 4
[38] 4 4 2 2 3 2 2 3 4 3 4 2 2 1 2 3 4 4 1 2 2 1 3 2 4 3 2 2 3 1 4 2 2 2 3 3 2
[75] 2 4 4 2 2 3 2 1 2 2 1 3 2 3 3 3 2 1 4 4 4 3 3 2 2 2 2 4 3 1 4 2 2 2 3 2 2
[112] 2 4 2 2 4 2 3 4 2 4 1 2 2 4 1 2 1 2 1 2 4 1 2 2 4 2
Within cluster sum of squares by cluster:
[1] 10970.5655 947.2973 909.4400 4747.3317
(between_SS / total_SS = 91.0 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" 21.1.4.3 Fit the Cluster Model to the Optimal Number of Tiers
21.1.4.3.1 Quarterbacks
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.
Code
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:
Here are the number of players that are in each of the four clusters (i.e., tiers):
21.1.4.3.2 Running Backs
Here are the number of players that are in each of the four clusters (i.e., tiers):
21.1.4.3.3 Wide Receivers
Here are the number of players that are in each of the four clusters (i.e., tiers):
21.1.4.3.4 Tight Ends
Here are the number of players that are in each of the four clusters (i.e., tiers):
21.1.4.4 Plot the Tiers
We can merge the player’s classification into the dataset and plot each player’s classification.
21.1.4.4.1 Quarterbacks
Code
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)21.1.4.4.2 Running Backs
Code
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)21.1.4.4.3 Wide Receivers
Code
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)21.1.4.4.4 Tight Ends
Code
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)21.1.5 Types of Wide Receivers
Code
# 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)Code
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:
Code
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()21.1.5.1 Identify the Number of WR Types
21.1.5.1.1 Model-Based Clustering
Code
Bayesian Information Criterion (BIC):
EII VII EEI VEI EVI VVI EEE
1 -9400.530 -9400.530 -5730.938 -5730.938 -5730.938 -5730.938 -5464.898
2 -8943.268 -8922.692 -5669.928 -5667.602 -5534.637 -5651.493 -5500.031
3 -8790.536 -8717.746 -5651.011 -5651.736 -5493.176 -5595.425 -5527.990
4 -8664.901 -8562.951 -5573.086 -5529.508 -5492.457 -5495.322 -5505.751
5 -8566.286 -8501.066 -5550.684 -5515.871 -5500.458 -5498.740 -5486.218
6 -8580.175 -8491.948 -5577.373 -5509.697 -5535.023 -5508.465 -5510.000
7 -8486.659 -8477.386 -5600.141 -5534.533 -5575.827 -5562.060 -5531.315
8 -8504.350 -8438.208 -5603.559 -5551.912 -5569.165 -5603.975 -5536.610
9 -8487.594 NA -5562.715 NA NA NA -5453.268
VEE EVE VVE EEV VEV EVV VVV
1 -5464.898 -5464.898 -5464.898 -5464.898 -5464.898 -5464.898 -5464.898
2 -5369.084 -5241.748 -5245.311 -5425.118 -5398.129 -5408.345 -5413.176
3 -5396.509 -5156.126 -5157.589 -5422.941 -5420.354 -5557.936 -5538.562
4 -5376.467 -5214.781 -5216.954 -5597.554 -5606.722 -5593.255 -5670.587
5 NA NA NA -5679.078 -5632.101 NA NA
6 NA NA NA -5807.691 -5849.058 NA NA
7 NA NA NA -6004.158 -6017.775 NA NA
8 NA NA NA -6107.885 -6187.817 NA NA
9 NA NA NA -6289.598 NA NA NA
Top 3 models based on the BIC criterion:
EVE,3 VVE,3 EVE,4
-5156.126 -5157.589 -5214.781 Best BIC values:
EVE,3 VVE,3 EVE,4
BIC -5156.126 -5157.588921 -5214.78133
BIC diff 0.000 -1.462805 -58.65521
Code
Integrated Complete-data Likelihood (ICL) criterion:
EII VII EEI VEI EVI VVI EEE
1 -9400.530 -9400.530 -5730.938 -5730.938 -5730.938 -5730.938 -5464.898
2 -8949.602 -8931.122 -5696.540 -5687.787 -5552.652 -5668.797 -5519.325
3 -8803.583 -8729.418 -5680.128 -5673.845 -5515.425 -5614.686 -5562.933
4 -8680.871 -8572.452 -5597.446 -5544.316 -5518.498 -5514.924 -5532.120
5 -8576.692 -8511.696 -5564.227 -5532.269 -5523.810 -5520.105 -5513.130
6 -8597.852 -8503.416 -5606.650 -5529.532 -5560.443 -5527.201 -5543.614
7 -8495.352 -8491.292 -5631.853 -5561.195 -5599.048 -5578.247 -5556.782
8 -8520.829 -8444.475 -5630.370 -5572.200 -5584.268 -5618.975 -5562.428
9 -8502.101 NA -5587.757 NA NA NA -5483.625
VEE EVE VVE EEV VEV EVV VVV
1 -5464.898 -5464.898 -5464.898 -5464.898 -5464.898 -5464.898 -5464.898
2 -5370.942 -5247.483 -5250.248 -5426.495 -5401.170 -5409.550 -5414.403
3 -5411.768 -5175.795 -5175.798 -5425.826 -5421.241 -5565.027 -5549.886
4 -5391.203 -5238.893 -5237.832 -5610.236 -5617.903 -5605.702 -5681.133
5 NA NA NA -5689.620 -5644.062 NA NA
6 NA NA NA -5811.401 -5853.960 NA NA
7 NA NA NA -6007.680 -6020.933 NA NA
8 NA NA NA -6109.909 -6189.465 NA NA
9 NA NA NA -6290.656 NA NA NA
Top 3 models based on the ICL criterion:
EVE,3 VVE,3 VVE,4
-5175.795 -5175.798 -5237.832 Best ICL values:
EVE,3 VVE,3 VVE,4
ICL -5175.795 -5.175798e+03 -5237.83232
ICL diff 0.000 -3.167186e-03 -62.03714
Based on the cluster analyses, it appears that three clusters are the best fit to the data.
We can also use the tidyLPA package (Rosenberg et al., 2018; Rosenberg & van Lissa, 2021), which provides an interface to the mclust package (Fraley et al., 2024; Scrucca et al., 2023).
Code
# 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_classestidyLPA analysis using mclust:
Model Classes AIC BIC Entropy prob_min prob_max n_min n_max BLRT_p
1 1 1 5677.99 5730.94 1.00 1.00 1.00 1.00 1.00
2 1 2 5587.56 5669.93 0.72 0.92 0.92 0.49 0.51 0.01
3 1 3 5539.23 5651.01 0.79 0.88 0.92 0.16 0.45 0.01
4 1 4 5431.89 5573.09 0.86 0.88 0.95 0.14 0.37 0.01
5 1 5 5380.07 5550.68 0.92 0.95 1.00 0.02 0.42 0.01
6 1 6 5377.34 5577.37 0.87 0.88 1.00 0.04 0.29 0.17
7 2 1 5677.99 5730.94 1.00 1.00 1.00 1.00 1.00
8 2 2 5542.65 5651.49 0.81 0.86 0.98 0.26 0.74 0.01
9 2 3 5430.69 5595.42 0.85 0.89 0.97 0.24 0.41 0.01
10 2 4 5274.70 5495.32 0.88 0.93 0.95 0.18 0.35 0.01
11 2 5 5222.23 5498.74 0.89 0.90 0.95 0.12 0.30 0.01
12 2 6 5176.06 5508.47 0.92 0.90 0.98 0.13 0.22 0.01
13 3 1 5378.05 5642.80 1.00 1.00 1.00 1.00 1.00
14 3 2 5383.77 5677.93 0.77 0.94 0.95 0.49 0.51 0.81
15 3 3 5382.31 5705.89 0.76 0.88 0.91 0.29 0.38 0.28
16 3 4 5330.65 5683.65 0.85 0.90 0.95 0.09 0.41 0.01
17 3 5 5281.70 5664.12 0.86 0.88 1.00 0.04 0.39 0.01
18 3 6 5276.07 5687.90 0.85 0.86 1.00 0.04 0.35 0.20Code
Compare tidyLPA solutions:
Model Classes AIC BIC Entropy BLRT_p
1 1 5677.988 5730.938 1.000
1 2 5587.562 5669.928 0.719 0.010
1 3 5539.228 5651.011 0.788 0.010
1 4 5431.887 5573.086 0.856 0.010
1 5 5380.069 5550.684 0.924 0.010
1 6 5377.341 5577.373 0.872 0.168
2 1 5677.988 5730.938 1.000
2 2 5542.652 5651.493 0.807 0.010
2 3 5430.693 5595.425 0.852 0.010
2 4 5274.699 5495.322 0.884 0.010
2 5 5222.226 5498.740 0.892 0.010
2 6 5176.060 5508.465 0.915 0.010
3 1 5378.049 5642.797 1.000
3 2 5383.766 5677.930 0.768 0.812
3 3 5382.309 5705.890 0.760 0.277
3 4 5330.653 5683.650 0.848 0.010
3 5 5281.703 5664.117 0.859 0.010
3 6 5276.069 5687.899 0.849 0.198
Best model according to AIC is Model 2 with 6 classes.
Best model according to BIC is Model 2 with 4 classes.
Best model according to Entropy is Model NA with NA classes.
Best model according to BLRT_p is Model NA with NA classes.
An analytic hierarchy process, based on the fit indices AIC, AWE, BIC, CLC, and KIC (Akogul & Erisoglu, 2017), suggests the best solution is Model 2 with 6 classes.Code
tidyLPA analysis using mclust:
Model Classes AIC BIC Entropy prob_min prob_max n_min n_max BLRT_p
1 3 3 5382.31 5705.89 0.76 0.88 0.91 0.29 0.38 0.2521.1.5.1.2 k-Means Clustering
Code
K-means clustering with 3 clusters of sizes 48, 41, 51
Cluster means:
vertical forty height_coalesced weight_coalesced adot.rec drop_percent.rec
1 36.07292 4.429167 70.77083 186.3958 10.62468 0.05175256
2 36.20732 4.476829 73.09756 205.4146 10.17858 0.04402159
3 36.40196 4.491176 73.92157 210.0196 11.24639 0.05172073
receiving_air_yards_per_rec brk_tkl.rec receptions_per_season
1 16.62610 8.479167 44.83089
2 15.96242 25.780488 76.58586
3 20.79698 4.666667 35.79269
Clustering vector:
[1] 2 2 1 3 1 3 2 2 2 3 3 2 1 2 3 2 3 2 3 1 3 2 3 3 2 3 2 1 3 2 2 1 3 2 2 1 1
[38] 1 3 2 2 1 3 2 3 3 3 2 3 3 3 1 1 3 2 2 2 1 1 2 1 2 3 1 1 2 3 1 1 3 1 1 1 3
[75] 1 1 2 2 2 3 3 3 3 3 1 2 3 1 3 3 1 1 1 3 1 2 2 3 2 1 3 2 2 3 1 3 3 1 1 3 2
[112] 1 1 1 3 1 3 1 3 3 2 1 3 1 3 2 3 1 1 1 3 2 1 2 1 1 1 3 2 3
Within cluster sum of squares by cluster:
[1] 15709.88 24181.95 15738.22
(between_SS / total_SS = 54.9 %)
Available components:
[1] "cluster" "centers" "totss" "withinss" "tot.withinss"
[6] "betweenss" "size" "iter" "ifault" 21.1.5.2 Fit the Cluster Model to the Optimal Number of WR Types
Code
----------------------------------------------------
Gaussian finite mixture model fitted by EM algorithm
----------------------------------------------------
Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components:
log-likelihood n df BIC ICL
-2355.689 140 90 -5156.126 -5175.795
Clustering table:
1 2 3
15 77 48 21.1.5.3 Plots of the Cluster Model
21.1.5.4 Interpreting the Clusters
1 2 3
15 77 48 Code
[,1] [,2] [,3]
type 1.00 2.00 3.00
vertical 36.47 36.51 35.72
forty 4.46 4.46 4.48
height_coalesced 73.13 72.36 72.81
weight_coalesced 206.13 198.22 202.60
adot.rec 12.01 10.73 10.30
drop_percent.rec 0.07 0.05 0.04
receiving_air_yards_per_rec 22.74 18.06 16.28
brk_tkl.rec 0.47 6.19 25.38
receptions_per_season 39.89 40.58 70.72Based 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:
Type 2 (“Consistent Contributor/WR2”) WRs:
Type 3 (“Deep Threat/Boom-or-Bust”) WRs:
And here are the clusters based on the tidyLPA (Rosenberg et al., 2018; Rosenberg & van Lissa, 2021) solution:
Code
Code
[,1] [,2] [,3]
Class 1.00 2.00 3.00
vertical 36.96 35.99 35.90
forty 4.40 4.42 4.55
height_coalesced 71.67 73.19 72.77
weight_coalesced 202.22 195.02 204.25
adot.rec 10.06 12.42 9.71
drop_percent.rec 0.05 0.06 0.04
receiving_air_yards_per_rec 15.29 21.24 17.04
brk_tkl.rec 16.52 10.19 10.60
receptions_per_season 49.82 50.39 52.01
type 2.33 2.19 2.21
model_number 3.00 3.00 3.00
classes_number 3.00 3.00 3.00
CPROB1 0.87 0.03 0.07
CPROB2 0.07 0.93 0.04
CPROB3 0.06 0.05 0.9021.2 Conclusion
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.
21.3 Session Info
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
time zone: UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0 dplyr_1.2.1
[5] purrr_1.2.2 readr_2.2.0 tidyr_1.3.2 tibble_3.3.1
[9] tidyverse_2.0.0 plotly_4.12.1 ggplot2_4.0.3 tidyLPA_2.0.2
[13] tidySEM_0.2.12 mclust_6.1.3 nflreadr_1.5.1 petersenlab_1.2.3
loaded via a namespace (and not attached):
[1] Rdpack_2.6.6 DBI_1.3.0 mnormt_2.1.2
[4] gridExtra_2.3.1 sandwich_3.1-3 rlang_1.3.0
[7] magrittr_2.0.5 otel_0.2.0 compiler_4.6.1
[10] vctrs_0.7.3 reshape2_1.4.5 quadprog_1.5-8
[13] crayon_1.5.3 pkgconfig_2.0.3 fastmap_1.2.0
[16] backports_1.5.1 labeling_0.4.3 pbivnorm_0.6.0
[19] pander_0.6.6 rmarkdown_2.32 tzdb_0.5.0
[22] nloptr_2.2.1 xfun_0.61 cachem_1.1.0
[25] jsonlite_2.0.0 progress_1.2.3 psych_2.6.5
[28] prettyunits_1.2.0 parallel_4.6.1 lavaan_0.7-2
[31] cluster_2.1.8.2 R6_2.6.1 stringi_1.8.9
[34] RColorBrewer_1.1-3 parallelly_1.48.0 car_3.1-5
[37] boot_1.3-32 rpart_4.1.27 Rcpp_1.1.2
[40] knitr_1.52 future.apply_1.20.2 zoo_1.9-0
[43] base64enc_0.1-6 timechange_0.4.0 Matrix_1.7-5
[46] splines_4.6.1 nnet_7.3-20 tidyselect_1.2.1
[49] rstudioapi_0.19.0 abind_1.4-8 yaml_2.3.12
[52] codetools_0.2-20 listenv_1.0.0 lattice_0.22-9
[55] nonnest2_0.5-9 plyr_1.8.9 withr_3.0.3
[58] S7_0.2.2 coda_0.19-4.1 evaluate_1.0.5
[61] foreign_0.8-91 future_1.75.0 fastDummies_1.7.6
[64] CompQuadForm_1.4.4 texreg_1.40 pillar_1.11.1
[67] carData_3.0-6 checkmate_2.3.4 stats4_4.6.1
[70] reformulas_0.4.4 generics_0.1.4 dbscan_1.2.6
[73] hms_1.1.4 mix_1.0-13 scales_1.4.0
[76] minqa_1.2.8 globals_0.19.1 xtable_1.8-8
[79] glue_1.8.1 Hmisc_5.3-0 tools_4.6.1
[82] data.table_1.18.6.1 lme4_2.0-6 gsubfn_0.7
[85] RANN_2.6.3 mvtnorm_1.4-2 grid_4.6.1
[88] mitools_2.7 crosstalk_1.2.2 MplusAutomation_1.3
[91] rbibutils_2.4.1 colorspace_2.1-3 nlme_3.1-169
[94] htmlTable_2.5.0 proto_1.0.0 Formula_1.2-6
[97] cli_3.6.6 viridisLite_0.4.3 gtable_0.3.6
[100] digest_0.6.39 progressr_1.0.0 htmlwidgets_1.6.4
[103] farver_2.1.2 memoise_2.0.1 htmltools_0.5.9
[106] lifecycle_1.0.5 httr_1.4.9 MASS_7.3-65