NCAA Consulting Project
2026-08-07
Section 0: Load Packages
Load in the packages that will be used throughout this document/RMD
library(dplyr)
library(purrr)
library(stringr)
library(tidyr)
library(lubridate)## Warning: package 'lubridate' was built under R version 4.4.3library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.4.3library(fmsb)## Warning: package 'fmsb' was built under R version 4.4.3library(ggrepel)## Warning: package 'ggrepel' was built under R version 4.4.3Section 1: Data Preparation
Section one involves downloading the NCAA volleyball data, standardizing team and opponent names, verifying complete matches, calculating comparable performance rates and create a single standardized record for each match.
Section 1.1: Loading the Pre-Scraped CSVs
Loading in the pre-scraped CSVs from the ncaavolleyballr github page. The function allows multiple years/seasons worth of data to be pulled.
# The base_url pointing towards the CSVs
base_url <- "https://media.githubusercontent.com/media/JeffreyRStevens/ncaavolleyballr/refs/heads/main/data-csv"
# Function which pulls the pre-scraped CSVs
load_csv_data <- function(sport, data_type, division, years) {
map_dfr(years, function(yr) {
url <- glue::glue("{base_url}/{sport}_{data_type}_div{division}_{yr}.csv")
readr::read_csv(url, show_col_types = FALSE) |>
mutate(
source_year = yr,
Season = as.character(yr)
)
})
}Section 1.2: Data-Preparation Functions
Defining multiple helper functions to be used later.
Section 1.2.1: Team Name and ID Helper Functions
These two helper functions standardize team names throughout the dataset and create consistent text identifiers for the construction of unique match IDs later on.
clean_team_name <- function(x) {
x |>
as.character() |>
str_replace_all("\u00A0", " ") |>
str_replace("^\\s*@\\s*", "") |>
str_replace("\\s+@.*$", "") |>
str_squish()
}
make_slug <- function(x) {
x |>
str_to_lower() |>
str_replace_all("[^a-z0-9]+", "-") |>
str_remove("^-") |>
str_remove("-$")
}Section 1.2.2: Verified Match Pairing Function
This function goes through match results and removes unusable rows, pairs up reciprocal team records and create a dataset which contains one row per match. A match is retained in the dataset only when the two team rows agree on the participant teams, the winner, score and number of sets.
build_verified_matches <- function(df, competition) {
# Parse result strings such as "W 3-1", "W, 3-1", or "L 1-3"
result_parts <- str_match(
coalesce(as.character(df$Result), ""),
regex(
"^\\s*([WL]).*?(\\d+)\\s*-\\s*(\\d+)",
ignore_case = TRUE
)
)
prepared <- df |>
mutate(
source_row_id = row_number(),
competition = .env$competition,
# Do not include an order without a year, because that could silently
# assign the wrong calendar year.
match_date = as.Date(
parse_date_time(
as.character(Date),
orders = c("ymd", "mdy", "dmy"),
quiet = TRUE
)
),
team_clean = clean_team_name(Team),
opponent_clean = clean_team_name(Opponent),
outcome = str_to_upper(result_parts[, 2]),
sets_won = as.integer(result_parts[, 3]),
sets_lost = as.integer(result_parts[, 4]),
win = case_when(
outcome == "W" ~ 1L,
outcome == "L" ~ 0L,
TRUE ~ NA_integer_
),
sets_played = suppressWarnings(
as.integer(as.character(S))
),
row_eligible =
!is.na(match_date) &
!is.na(win) &
!is.na(sets_won) &
!is.na(sets_lost) &
!is.na(team_clean) &
!is.na(opponent_clean) &
team_clean != "" &
opponent_clean != "" &
team_clean != opponent_clean
)
# Cancellations, postponements, malformed scores, and bad dates
excluded_rows <- prepared |>
filter(!row_eligible) |>
mutate(
exclusion_reason = case_when(
!outcome %in% c("W", "L") ~
"Not completed or result could not be parsed",
is.na(match_date) ~
"Date could not be parsed",
is.na(team_clean) | team_clean == "" ~
"Missing team name",
is.na(opponent_clean) | opponent_clean == "" ~
"Missing opponent name",
team_clean == opponent_clean ~
"Team and opponent are identical",
TRUE ~
"Other invalid row"
)
)
eligible <- prepared |>
filter(row_eligible) |>
mutate(
# Alphabetical orientation; identical on both team rows
team_low = pmin(team_clean, opponent_clean),
team_high = pmax(team_clean, opponent_clean),
# Determine the actual winner from either team's perspective
winner_team = if_else(
win == 1L,
team_clean,
opponent_clean
),
loser_team = if_else(
win == 0L,
team_clean,
opponent_clean
),
# Express the score from the winner's perspective
winner_sets = if_else(
win == 1L,
sets_won,
sets_lost
),
loser_sets = if_else(
win == 1L,
sets_lost,
sets_won
),
# Broad key for all rows involving these teams on this date
event_key = paste(
competition,
Season,
match_date,
team_low,
team_high,
sep = " | "
),
# More specific signature for repeated same-day matchups
match_signature = paste(
event_key,
winner_team,
paste0(winner_sets, "-", loser_sets),
sep = " | "
)
)
# Number of rows associated with each broad event
event_sizes <- eligible |>
count(event_key, name = "event_rows")
pair_checks <- eligible |>
group_by(event_key, match_signature) |>
summarise(
n_rows = n(),
n_teams = n_distinct(team_clean),
n_winners = sum(win == 1L),
n_losers = sum(win == 0L),
reciprocal_names = setequal(
team_clean,
opponent_clean
),
score_mirrors = {
winner_position <- which(win == 1L)
loser_position <- which(win == 0L)
length(winner_position) == 1L &&
length(loser_position) == 1L &&
isTRUE(
sets_won[winner_position] ==
sets_lost[loser_position]
) &&
isTRUE(
sets_lost[winner_position] ==
sets_won[loser_position]
)
},
total_sets_agree = all(
is.na(sets_played) |
sets_played == winner_sets + loser_sets
),
.groups = "drop"
) |>
left_join(event_sizes, by = "event_key") |>
mutate(
pair_status = case_when(
n_rows > 2 ~
"Ambiguous duplicate: same teams, date, winner, and score",
n_rows == 1 & event_rows == 1 ~
"Missing reciprocal row or team-name mismatch",
n_rows == 1 & event_rows > 1 ~
"Unmatched row within repeated same-day matchup",
n_rows != 2 ~
"Incorrect number of rows",
n_teams != 2 ~
"Two rows do not represent two different teams",
n_winners != 1 | n_losers != 1 ~
"Outcome labels do not contain one winner and one loser",
!reciprocal_names ~
"Team and opponent names are not reciprocal",
!score_mirrors ~
"Set scores do not mirror each other",
!total_sets_agree ~
"Reported number of sets is inconsistent",
TRUE ~
"Valid pair"
)
)
valid_signatures <- pair_checks |>
filter(pair_status == "Valid pair") |>
select(match_signature)
paired_rows <- eligible |>
semi_join(valid_signatures, by = "match_signature") |>
mutate(
match_id = paste(
competition,
Season,
format(match_date, "%Y%m%d"),
make_slug(team_low),
make_slug(team_high),
make_slug(winner_team),
paste0(winner_sets, "-", loser_sets),
sep = "__"
)
)
# Completed-looking rows that could not be paired reliably
exception_rows <- eligible |>
inner_join(
pair_checks |>
filter(pair_status != "Valid pair") |>
select(match_signature, pair_status),
by = "match_signature"
)
# One row per verified match
match_table <- paired_rows |>
group_by(match_id) |>
summarise(
competition = first(competition),
Season = first(Season),
Date = first(match_date),
# Canonical orientation rather than always placing the winner first
team_a = first(team_low),
team_b = first(team_high),
team_a_win = as.integer(
first(winner_team) == first(team_low)
),
winner = first(winner_team),
loser = first(loser_team),
winner_sets = first(winner_sets),
loser_sets = first(loser_sets),
winner_source_row = source_row_id[win == 1L],
loser_source_row = source_row_id[win == 0L],
.groups = "drop"
)
# Stop immediately if an internal assumption has failed
stopifnot(
!anyDuplicated(match_table$match_id),
nrow(paired_rows) == 2 * nrow(match_table),
n_distinct(paired_rows$source_row_id) == nrow(paired_rows),
all(count(paired_rows, match_id)$n == 2)
)
list(
prepared_rows = prepared,
paired_rows = paired_rows,
matches = match_table,
pair_checks = pair_checks,
exceptions = exception_rows,
excluded = excluded_rows
)
}Section 1.2.3: Canonical Match Feature Function
This function combines the two verified team records into a single standardized match row and then calculates Team A minus Team B performance differentials.
build_match_feature_table <- function(team_rates) {
team_a_rows <- team_rates |>
filter(team_clean == team_low) |>
transmute(
match_id,
Season,
match_date,
team_a = team_low,
team_b = team_high,
team_a_win = win,
team_a_kill_rate = kill_rate_pct,
team_a_error_rate = error_rate_pct,
team_a_block_rate = block_rate_pct,
team_a_aces_per_set = aces_per_set,
team_a_digs_per_set = digs_per_set
)
team_b_rows <- team_rates |>
filter(team_clean == team_high) |>
transmute(
match_id,
team_b_kill_rate = kill_rate_pct,
team_b_error_rate = error_rate_pct,
team_b_block_rate = block_rate_pct,
team_b_aces_per_set = aces_per_set,
team_b_digs_per_set = digs_per_set
)
match_features <- team_a_rows |>
inner_join(team_b_rows, by = "match_id") |>
mutate(
diff_kill_rate =
team_a_kill_rate - team_b_kill_rate,
diff_error_rate =
team_a_error_rate - team_b_error_rate,
diff_block_rate =
team_a_block_rate - team_b_block_rate,
diff_aces =
team_a_aces_per_set - team_b_aces_per_set,
diff_digs =
team_a_digs_per_set - team_b_digs_per_set
)
stopifnot(
!anyDuplicated(match_features$match_id),
nrow(match_features) == n_distinct(team_rates$match_id),
all(match_features$team_a_win %in% c(0L, 1L))
)
match_features
}Section 1.2.4: Safe rate calculation
Helps safely calculate rates while preventing errors caused when dividing by zero or using invalid denominators.
safe_divide <- function(numerator, denominator) {
if_else(
!is.na(denominator) & denominator > 0,
numerator / denominator,
NA_real_
)
}Section 1.3: Season-Level Data
Dealing with the season-level CSVs which provide player and team statistics used in the exploratory data analysis and main analysis sections.
Section 1.3.1: Women’s Volleyball
Loading in the division 1 women’s volleyball player-season and team-season statistics from 2020 to 2025.
# Women's player season stats (2020-2025)
wvb_player_season <- load_csv_data("wvb", "playerseason", 1, 2020:2025)
# Women's team season stats (2020-2025)
wvb_team_season <- load_csv_data("wvb", "teamseason", 1, 2020:2025)Section 1.3.2: Men’s Volleyball
Loading in the division 1 women’s volleyball player-season and team-season statistics from 2020 to 2024. The amount of data provided for the men is a little less than for women, requiring one less year here.
# Men's player season stats (2020-2024)
mvb_player_season <- load_csv_data("mvb", "playerseason", 1, 2020:2024)
# Men's team season stats (2020-2024)
mvb_team_season <- load_csv_data("mvb", "teamseason", 1, 2020:2024)Section 1.3.3: Women’s Volleyball Validation
These checks confirm seasonal data coverage and duplicate player-season and team-season records for women’s volleyball.
Section 1.3.4: Men’s Volleyball Validation
These checks confirm seasonal data coverage and duplicate player-season and team-season records for men’s volleyball.
Section 1.4: Women’s Match-Level Data and Verification
This section focuses on converting the original women’s match data into verified, model-ready match tables.
Section 1.4.1: Loading The Data
Utilizes the csv loading helper function to load in the women’s team-match and player-match records for the indicated years.
wvb_team_match <- load_csv_data(
sport = "wvb",
data_type = "teammatch",
division = 1,
years = 2020:2025
)Section 1.4.2: Pair Verification
Application of the verified matches builder function to create two-row team-match and one-row match tables. It reports unmatched, excluded or inconsistent records.
# Match the two team rows belonging to each completed match
wvb_pairing <- build_verified_matches(
df = wvb_team_match,
competition = "wvb_d1"
)
# Two verified team rows per match.
# Retaining the old object name prevents downstream code from breaking.
wvb_team_match_clean <- wvb_pairing$paired_rows
# One row per verified match
wvb_matches <- wvb_pairing$matches
# Summary of pair-validation results
wvb_pairing$pair_checks |>
count(pair_status, sort = TRUE)## # A tibble: 4 × 2
## pair_status n
## <chr> <int>
## 1 Valid pair 25706
## 2 Missing reciprocal row or team-name mismatch 3051
## 3 Ambiguous duplicate: same teams, date, winner, and score 3
## 4 Reported number of sets is inconsistent 1# Rows that appear to describe completed matches but could not be paired
wvb_pairing$exceptions |>
select(
source_row_id,
Season,
match_date,
Team,
Opponent,
Result,
pair_status
) |>
arrange(Season, match_date, Team)## # A tibble: 3,065 × 7
## source_row_id Season match_date Team Opponent Result pair_status
## <int> <chr> <date> <chr> <chr> <chr> <chr>
## 1 4568 2020 2020-09-19 Lamar University @ Louisi… L 0-3 Ambiguous …
## 2 4569 2020 2020-09-19 Lamar University @ Louisi… L 0-3 Ambiguous …
## 3 4997 2020 2020-09-19 Louisiana Lamar Un… W 3-0 Ambiguous …
## 4 4998 2020 2020-09-19 Louisiana Lamar Un… W 3-0 Ambiguous …
## 5 4948 2020 2020-11-18 Ga. Southern #2 Troy … W 3-2 Missing re…
## 6 4990 2020 2020-11-18 Little Rock #2 Louis… W 3-1 Missing re…
## 7 5015 2020 2020-11-18 Louisiana #5 Littl… L 1-3 Missing re…
## 8 5062 2020 2020-11-18 Texas St. #6 ULM @… W 3-0 Missing re…
## 9 5099 2020 2020-11-18 Troy #5 Ga. S… L 2-3 Missing re…
## 10 5124 2020 2020-11-18 ULM #1 Texas… L 0-3 Missing re…
## # ℹ 3,055 more rows# Rows excluded before matching, such as canceled matches
wvb_pairing$excluded |>
count(exclusion_reason, sort = TRUE)## # A tibble: 2 × 2
## exclusion_reason n
## <chr> <int>
## 1 Not completed or result could not be parsed 998
## 2 Date could not be parsed 394# Final structural checks
stopifnot(
nrow(wvb_team_match_clean) == 2 * nrow(wvb_matches),
!anyDuplicated(wvb_matches$match_id),
all(count(wvb_team_match_clean, match_id)$n == 2)
)Section 1.4.3: Rate Calculations
Calculation of attacking, error, blocking, serving and digging rates for each verified women’s team-match record.
wvb_team_rates <- wvb_team_match_clean |>
mutate(
Block_Solos_clean =
coalesce(as.numeric(`Block Solos`), 0),
Block_Assists_clean =
coalesce(as.numeric(`Block Assists`), 0),
Aces_clean =
coalesce(as.numeric(Aces), 0),
Digs_clean =
coalesce(as.numeric(Digs), 0),
sets_num = as.numeric(S),
kill_rate_pct =
100 * safe_divide(Kills, `Total Attacks`),
error_rate_pct =
100 * safe_divide(Errors, `Total Attacks`),
block_rate_pct =
100 * safe_divide(
Block_Solos_clean + 0.5 * Block_Assists_clean,
`Total Attacks`
),
aces_per_set =
safe_divide(Aces_clean, sets_num),
digs_per_set =
safe_divide(Digs_clean, sets_num)
)Section 1.4.4: Final One-Row Match Table
The creation of the final table featuring one row per match. A check of the distribution of Team A wins and losses is also shown.
wvb_match_features <- build_match_feature_table(
wvb_team_rates
)
wvb_match_features |>
count(team_a_win)## # A tibble: 2 × 2
## team_a_win n
## <int> <int>
## 1 0 13037
## 2 1 12669Section 1.5: Men’s Match-Level Data and Verification
The same process of verification and constructing of features done in section 1.4 is repeated for division 1 men’s volleyball.
Section 1.5.1: Loading The Data
Uses the csv loading helper function to load in the men’s team-match and player-match records for the indicated years.
mvb_team_match <- load_csv_data(
sport = "mvb",
data_type = "teammatch",
division = 1,
years = 2020:2024
)Section 1.5.2: Pair Verification
Applies the verified matches builder function to create two-row team-match and one-row match tables. Includes a summary of unmatched, excluded or inconsistent records.
# Match the two team rows belonging to each completed match
mvb_pairing <- build_verified_matches(
df = mvb_team_match,
competition = "mvb_d1"
)
# Two verified team rows per match
mvb_team_match_clean <- mvb_pairing$paired_rows
# One row per verified match
mvb_matches <- mvb_pairing$matches
mvb_pairing$pair_checks |>
count(pair_status, sort = TRUE)## # A tibble: 3 × 2
## pair_status n
## <chr> <int>
## 1 Valid pair 3012
## 2 Missing reciprocal row or team-name mismatch 974
## 3 Two rows do not represent two different teams 2mvb_pairing$exceptions |>
select(
source_row_id,
Season,
match_date,
Team,
Opponent,
Result,
pair_status
) |>
arrange(Season, match_date, Team)## # A tibble: 978 × 7
## source_row_id Season match_date Team Opponent Result pair_status
## <int> <chr> <date> <chr> <chr> <chr> <chr>
## 1 728 2020 2021-01-13 Quincy Indiana… W 3-0 Missing re…
## 2 411 2020 2021-01-16 Tusculum @ Campb… L 0-3 Missing re…
## 3 284 2020 2021-01-19 Alderson Broaddus Mount U… L 0-3 Missing re…
## 4 614 2020 2021-01-20 Lewis Campbel… W 3-1 Missing re…
## 5 376 2020 2021-01-21 Queens (NC) @ Lourd… L 2-3 Missing re…
## 6 730 2020 2021-01-21 Quincy Georget… W 3-0 Missing re…
## 7 167 2020 2021-01-22 Erskine Maranat… W 3-0 Missing re…
## 8 638 2020 2021-01-22 Lindenwood Georget… W 3-2 Missing re…
## 9 285 2020 2021-01-23 Alderson Broaddus @ St. A… W 3-1 Missing re…
## 10 168 2020 2021-01-23 Erskine Lancast… W 3-0 Missing re…
## # ℹ 968 more rowsmvb_pairing$excluded |>
count(exclusion_reason, sort = TRUE)## # A tibble: 2 × 2
## exclusion_reason n
## <chr> <int>
## 1 Not completed or result could not be parsed 361
## 2 Date could not be parsed 78stopifnot(
nrow(mvb_team_match_clean) == 2 * nrow(mvb_matches),
!anyDuplicated(mvb_matches$match_id),
all(count(mvb_team_match_clean, match_id)$n == 2)
)Section 1.5.3: Rate Calculations
Calculates attacking, error, blocking, serving and digging rates for each verified men’s team-match record.
# 1. Calculate rate and per-set metrics on mvb_team_match_clean
mvb_team_rates <- mvb_team_match_clean |>
mutate(
Block_Solos_clean =
coalesce(as.numeric(`Block Solos`), 0),
Block_Assists_clean =
coalesce(as.numeric(`Block Assists`), 0),
Aces_clean =
coalesce(as.numeric(Aces), 0),
Digs_clean =
coalesce(as.numeric(Digs), 0),
sets_num = as.numeric(S),
kill_rate_pct =
100 * safe_divide(Kills, `Total Attacks`),
error_rate_pct =
100 * safe_divide(Errors, `Total Attacks`),
block_rate_pct =
100 * safe_divide(
Block_Solos_clean + 0.5 * Block_Assists_clean,
`Total Attacks`
),
aces_per_set =
safe_divide(Aces_clean, sets_num),
digs_per_set =
safe_divide(Digs_clean, sets_num)
)Section 1.5.4: Final One-Row Match Table
Creates the final feature table consisting of one row per match and checks Team A’s win/loss distribution.
mvb_match_features <- build_match_feature_table(
mvb_team_rates
)
mvb_match_features |>
count(team_a_win)## # A tibble: 2 × 2
## team_a_win n
## <int> <int>
## 1 0 1603
## 2 1 1409Section 1.6: Validation Summary
Final checks summarize the amount of data which was retained and confirm that the match tables structural requirements.
Section 1.6.1: Overall Summary
Summary of the raw, excluded, exception, verified and final match-record counts for all men’s and women’s volleyball seasons considered in this analysis.
validation_summary <- bind_rows(
tibble(
competition = "Women's Division I",
raw_team_rows = nrow(wvb_team_match),
excluded_rows = nrow(wvb_pairing$excluded),
exception_rows = nrow(wvb_pairing$exceptions),
verified_team_rows = nrow(wvb_team_match_clean),
verified_matches = nrow(wvb_matches),
canonical_matches = nrow(wvb_match_features)
),
tibble(
competition = "Men's Division I",
raw_team_rows = nrow(mvb_team_match),
excluded_rows = nrow(mvb_pairing$excluded),
exception_rows = nrow(mvb_pairing$exceptions),
verified_team_rows = nrow(mvb_team_match_clean),
verified_matches = nrow(mvb_matches),
canonical_matches = nrow(mvb_match_features)
)
) |>
mutate(
verified_row_share =
verified_team_rows / raw_team_rows,
rows_per_verified_match =
verified_team_rows / verified_matches,
percent_raw_rows_verified =
scales::percent(
verified_row_share,
accuracy = 0.1
),
percent_candidate_pairs_valid =
case_when(
competition == "Women's Division I" ~
scales::percent(
mean(
wvb_pairing$pair_checks$pair_status ==
"Valid pair"
),
accuracy = 0.1
),
competition == "Men's Division I" ~
scales::percent(
mean(
mvb_pairing$pair_checks$pair_status ==
"Valid pair"
),
accuracy = 0.1
)
)
)
validation_summary## # A tibble: 2 × 11
## competition raw_team_rows excluded_rows exception_rows verified_team_rows
## <chr> <int> <int> <int> <int>
## 1 Women's Divisio… 55869 1392 3065 51412
## 2 Men's Division I 7441 439 978 6024
## # ℹ 6 more variables: verified_matches <int>, canonical_matches <int>,
## # verified_row_share <dbl>, rows_per_verified_match <dbl>,
## # percent_raw_rows_verified <chr>, percent_candidate_pairs_valid <chr>Section 1.6.2: Pair Status Summary
Shows how candidate match pairs were classified during the verification process.
Section 1.6.3: Exclusion Summary
Summary of why raw match records were excluded before the pairing process.
Section 1.6.4: Row Verification
Checks that every verified match has exactly two team rows, one unique final match row and a valid binary outcome.
stopifnot(
# Exactly two verified team rows per match
nrow(wvb_team_match_clean) == 2 * nrow(wvb_matches),
nrow(mvb_team_match_clean) == 2 * nrow(mvb_matches),
# Every verified match ID appears exactly twice in the team-row data
all(count(wvb_team_match_clean, match_id)$n == 2),
all(count(mvb_team_match_clean, match_id)$n == 2),
# Canonical match IDs are unique
!anyDuplicated(wvb_matches$match_id),
!anyDuplicated(mvb_matches$match_id),
!anyDuplicated(wvb_match_features$match_id),
!anyDuplicated(mvb_match_features$match_id),
# Canonical tables contain one row for every verified match
nrow(wvb_match_features) == nrow(wvb_matches),
nrow(mvb_match_features) == nrow(mvb_matches),
# Outcomes are valid binary indicators
all(wvb_match_features$team_a_win %in% c(0L, 1L)),
all(mvb_match_features$team_a_win %in% c(0L, 1L)),
# Match identities agree across tables
setequal(wvb_matches$match_id, wvb_match_features$match_id),
setequal(mvb_matches$match_id, mvb_match_features$match_id) )Section 1.6.5: Outcome Balance
Confirms that alphabetically assigning Team A doesn’t create a super imbalanced win-loss outcome.
Section 1.6.6: Missingness Check
Counts missing values in the final match-level variables used in later sections.
Section 2: Exploratory Data Analysis
This section explores the patterns in coverage, competitiveness, team performance and player participation in the verified data from section 1. Results from this section also help with decision making in later sections.
Section 2.1: EDA Data Construction
Section 2.1.1: Combined EDA Datasets
Four combined dataset were created in this section: one row per verified match, two team rows per match, one feature row per match and one row per player-season.
# One row per verified match
eda_matches <- bind_rows(
wvb_matches |>
mutate(competition = "Women's Division I"),
mvb_matches |>
mutate(competition = "Men's Division I")
) |>
mutate(
total_sets = winner_sets + loser_sets,
scoreline = paste0(winner_sets, "-", loser_sets)
)
# Two verified team rows per match
eda_team_rates <- bind_rows(
wvb_team_rates |>
mutate(competition = "Women's Division I"),
mvb_team_rates |>
mutate(competition = "Men's Division I")
)
# One canonical feature row per match
eda_match_features <- bind_rows(
wvb_match_features |>
mutate(competition = "Women's Division I"),
mvb_match_features |>
mutate(competition = "Men's Division I")
)
# Player-season data
eda_player_season <- bind_rows(
wvb_player_season |>
mutate(competition = "Women's Division I"),
mvb_player_season |>
mutate(competition = "Men's Division I")
)
eda_matches <- eda_matches |>
mutate(
competition = factor(
competition,
levels = c(
"Women's Division I",
"Men's Division I"
)
)
)Section 2.1.2: Structural Checks
Automated check to confirm that the combined tables contain the expected amount of rows, retain unique match identifiers and calculate total sets correctly.
stopifnot(
nrow(eda_matches) ==
nrow(wvb_matches) + nrow(mvb_matches),
nrow(eda_match_features) ==
nrow(wvb_match_features) + nrow(mvb_match_features),
!anyDuplicated(eda_matches$match_id),
!anyDuplicated(eda_match_features$match_id),
all(eda_matches$total_sets ==
eda_matches$winner_sets +
eda_matches$loser_sets)
)Section 2.2: Data Coverage
Summary of the number of verified matches, participating teams and observed data range for each competition and season.
Section 2.2.1: Verified Matches Plot
Compares the number of verified matches available across seasons and highlights differences in data volume between the women’s and men’s competitions.
match_coverage <- eda_matches |>
group_by(competition, Season) |>
summarise(
verified_matches = n(),
first_match = min(Date, na.rm = TRUE),
last_match = max(Date, na.rm = TRUE),
.groups = "drop"
)
team_coverage <- eda_matches |>
select(
competition,
Season,
team_a,
team_b
) |>
pivot_longer(
cols = c(team_a, team_b),
names_to = "team_position",
values_to = "team"
) |>
distinct(competition, Season, team) |>
count(
competition,
Season,
name = "participating_teams"
)
coverage_summary <- match_coverage |>
left_join(
team_coverage,
by = c("competition", "Season")
)
ggplot(
coverage_summary,
aes(
x = factor(Season),
y = verified_matches,
group = competition,
linetype = competition
)
) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
facet_wrap(
vars(competition),
scales = "free_y"
) +
labs(
title = "Verified Matches by Season",
subtitle = "Only matches passing the Section 1 pairing checks are included",
x = "Season",
y = "Number of verified matches"
) +
theme_minimal() +
theme(
legend.position = "none"
)
coverage_summary## # A tibble: 11 × 6
## competition Season verified_matches first_match last_match
## <fct> <chr> <int> <date> <date>
## 1 Women's Division I 2020 2241 2020-09-04 2021-04-14
## 2 Women's Division I 2021 4723 2021-08-27 2021-12-07
## 3 Women's Division I 2022 4875 2022-08-26 2022-12-05
## 4 Women's Division I 2023 4865 2023-08-25 2023-12-01
## 5 Women's Division I 2024 4867 2024-08-27 2024-12-11
## 6 Women's Division I 2025 4135 2025-08-22 2025-11-29
## 7 Men's Division I 2020 345 2021-01-14 2021-05-04
## 8 Men's Division I 2021 601 2022-01-05 2022-05-03
## 9 Men's Division I 2022 664 2023-01-04 2023-05-06
## 10 Men's Division I 2023 668 2024-01-03 2024-04-19
## 11 Men's Division I 2024 734 2025-01-03 2025-04-26
## # ℹ 1 more variable: participating_teams <int>Section 2.2.2: Player Coverage Summary
Summary of the number of teams and player-seasons represented in each season and measures the amount of missing position information. Position-data coverage is important since the player efficiency index later on is calculated separately for each position group.
player_coverage <- eda_player_season |>
group_by(competition, Season) |>
summarise(
teams = n_distinct(Team),
player_seasons = n_distinct(
paste(Team, Player, Number, sep = " | ")
),
missing_position = sum(
is.na(Pos) | str_squish(Pos) == ""
),
missing_position_pct = mean(
is.na(Pos) | str_squish(Pos) == ""
),
.groups = "drop"
)
player_coverage## # A tibble: 11 × 6
## competition Season teams player_seasons missing_position missing_position_pct
## <chr> <chr> <int> <int> <int> <dbl>
## 1 Men's Divi… 2020 43 644 107 0.166
## 2 Men's Divi… 2021 55 857 63 0.0735
## 3 Men's Divi… 2022 57 941 23 0.0244
## 4 Men's Divi… 2023 59 988 12 0.0121
## 5 Men's Divi… 2024 66 1095 18 0.0164
## 6 Women's Di… 2020 312 4545 33 0.00726
## 7 Women's Di… 2021 340 5400 4 0.000741
## 8 Women's Di… 2022 344 5410 11 0.00203
## 9 Women's Di… 2023 344 5361 5 0.000933
## 10 Women's Di… 2024 346 5472 9 0.00164
## 11 Women's Di… 2025 348 5485 10 0.00182Section 2.3: Match Results and Match Length
Section 2.3.1: Scoreline Summary
Calculates the frequency and proportion of each final match scoreline (3-0, 3-1, 3-2).
scoreline_summary <- eda_matches |>
count(
competition,
scoreline,
total_sets,
name = "matches"
) |>
group_by(competition) |>
mutate(
proportion = matches / sum(matches)
) |>
ungroup() |>
arrange(
competition,
total_sets,
scoreline
)
scoreline_summary## # A tibble: 7 × 5
## competition scoreline total_sets matches proportion
## <fct> <chr> <int> <int> <dbl>
## 1 Women's Division I 2-1 3 1 0.0000389
## 2 Women's Division I 3-0 3 12441 0.484
## 3 Women's Division I 3-1 4 8150 0.317
## 4 Women's Division I 3-2 5 5114 0.199
## 5 Men's Division I 3-0 3 1581 0.525
## 6 Men's Division I 3-1 4 893 0.296
## 7 Men's Division I 3-2 5 538 0.179Section 2.3.2: Match Length By Season
This table summarizes the average match length and seasonal shares of straight-set and five-set matches.
match_length_summary <- eda_matches |>
group_by(competition, Season) |>
summarise(
matches = n(),
mean_sets = mean(total_sets),
median_sets = median(total_sets),
five_set_matches = sum(total_sets == 5),
five_set_share = mean(total_sets == 5),
straight_set_share = mean(total_sets == 3),
.groups = "drop"
)
match_length_summary## # A tibble: 11 × 8
## competition Season matches mean_sets median_sets five_set_matches
## <fct> <chr> <int> <dbl> <dbl> <int>
## 1 Women's Division I 2020 2241 3.74 4 479
## 2 Women's Division I 2021 4723 3.68 4 864
## 3 Women's Division I 2022 4875 3.73 4 1008
## 4 Women's Division I 2023 4865 3.70 4 959
## 5 Women's Division I 2024 4867 3.72 4 995
## 6 Women's Division I 2025 4135 3.72 4 809
## 7 Men's Division I 2020 345 3.66 3 64
## 8 Men's Division I 2021 601 3.66 3 115
## 9 Men's Division I 2022 664 3.69 3 127
## 10 Men's Division I 2023 668 3.65 3 125
## 11 Men's Division I 2024 734 3.62 3 107
## # ℹ 2 more variables: five_set_share <dbl>, straight_set_share <dbl>Section 2.3.3: Scoreline Trends By Season
Shows how final scoreline distributions changed across the seasons considered for women’s and men’s volleyball.
scoreline_by_season <- eda_matches |>
count(
competition,
Season,
scoreline,
name = "matches"
) |>
group_by(competition, Season) |>
mutate(
proportion = matches / sum(matches)
) |>
ungroup()
ggplot(
scoreline_by_season,
aes(
x = factor(Season),
y = proportion,
fill = scoreline
)
) +
geom_col() +
geom_text(
data = scoreline_by_season |>
dplyr::filter(proportion >= 0.03),
aes(
label = scales::percent(proportion, accuracy = 0.1)
),
position = position_stack(vjust = 0.5),
color = "white",
size = 4
) +
scale_fill_manual(
breaks = c("3-0", "3-1", "3-2"),
values = c(
"3-0" = "#4E79A7",
"3-1" = "#59A14F",
"3-2" = "#E15759"
),
drop = TRUE
) +
facet_wrap(vars(competition)) +
scale_y_continuous(
labels = scales::label_percent()
) +
labs(
title = "Match Scorelines by Season",
x = "Season",
y = "Percentage of matches",
fill = "Scoreline"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 20, face = "bold"),
strip.text = element_text(size = 16),
axis.title = element_text(size = 20),
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 12),
axis.text.y = element_text(size = 12),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14)
)
Section 2.4: Team Performance Trends
Section 2.4.1: Seasonal Performance Summary
Calculates season-level attacking, serving, blocking and defensive rates from all verified team-match records.
season_performance <- eda_team_rates |>
group_by(competition, Season) |>
summarise(
team_match_rows = n(),
attack_efficiency = safe_divide(
sum(Kills, na.rm = TRUE) -
sum(Errors, na.rm = TRUE),
sum(`Total Attacks`, na.rm = TRUE)
),
kills_per_set = safe_divide(
sum(Kills, na.rm = TRUE),
sum(sets_num, na.rm = TRUE)
),
attack_errors_per_set = safe_divide(
sum(Errors, na.rm = TRUE),
sum(sets_num, na.rm = TRUE)
),
aces_per_set = safe_divide(
sum(Aces_clean, na.rm = TRUE),
sum(sets_num, na.rm = TRUE)
),
digs_per_set = safe_divide(
sum(Digs_clean, na.rm = TRUE),
sum(sets_num, na.rm = TRUE)
),
blocks_per_set = safe_divide(
sum(
Block_Solos_clean +
0.5 * Block_Assists_clean,
na.rm = TRUE
),
sum(sets_num, na.rm = TRUE)
),
.groups = "drop"
)
season_performance## # A tibble: 11 × 9
## competition Season team_match_rows attack_efficiency kills_per_set
## <chr> <chr> <int> <dbl> <dbl>
## 1 Men's Division I 2020 690 0.254 11.7
## 2 Men's Division I 2021 1202 0.256 11.4
## 3 Men's Division I 2022 1328 0.255 11.4
## 4 Men's Division I 2023 1336 0.259 11.4
## 5 Men's Division I 2024 1468 0.252 11.5
## 6 Women's Division I 2020 4482 0.200 12.2
## 7 Women's Division I 2021 9446 0.205 12.3
## 8 Women's Division I 2022 9750 0.206 12.3
## 9 Women's Division I 2023 9730 0.207 12.2
## 10 Women's Division I 2024 9734 0.208 12.4
## 11 Women's Division I 2025 8270 0.211 12.4
## # ℹ 4 more variables: attack_errors_per_set <dbl>, aces_per_set <dbl>,
## # digs_per_set <dbl>, blocks_per_set <dbl>Section 2.4.2: Reshape The Performance Data
Conversion of the seasonal performance table into a long format. This allows multiple metrics to be displayed together.
season_performance_long <- season_performance |>
select(
-team_match_rows
) |>
pivot_longer(
cols = -c(competition, Season),
names_to = "metric",
values_to = "value"
) |>
mutate(
metric = recode(
metric,
attack_efficiency = "Attack efficiency",
kills_per_set = "Kills per set",
attack_errors_per_set = "Attack errors per set",
aces_per_set = "Aces per set",
digs_per_set = "Digs per set",
blocks_per_set = "Blocks per set"
)
)Section 2.4.3: Team Performance Trend Plot
Compares how attack efficiency, kills, errors, serving, digging, and blocking rates changed across seasons and competitions.
ggplot(
season_performance_long,
aes(
x = Season,
y = value,
group = competition,
linetype = competition
)
) +
geom_line(linewidth = 0.9) +
geom_point(size = 2) +
facet_wrap(
vars(metric),
scales = "free_y"
) +
labs(
title = "Team Performance Trends by Season",
subtitle = "Rates are calculated from pooled opportunities within each season",
x = "Season",
y = NULL,
linetype = "Competition"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 10)
)
Section 2.5: Winner-Loser Performance Differences
Section 2.5.1: Orient Differences From Winner To Loser
Reorients match-level statistics so that each value represents the winner’s performance minus the loser’s performance. Positive values mean the winning team had the higher value whereas negative values indicate the loser had the higher value.
winner_loser_differences <- eda_match_features |>
mutate(
winner_minus_loser_kill_rate = if_else(
team_a_win == 1L,
diff_kill_rate,
-diff_kill_rate
),
winner_minus_loser_error_rate = if_else(
team_a_win == 1L,
diff_error_rate,
-diff_error_rate
),
winner_minus_loser_block_rate = if_else(
team_a_win == 1L,
diff_block_rate,
-diff_block_rate
),
winner_minus_loser_aces = if_else(
team_a_win == 1L,
diff_aces,
-diff_aces
),
winner_minus_loser_digs = if_else(
team_a_win == 1L,
diff_digs,
-diff_digs
)
)Section 2.5.2: Convert Differences To Long Format
Reshape the winner-loser performance differences into long format for metric summaries and visual comparisons.
winner_loser_long <- winner_loser_differences |>
select(
competition,
Season,
match_id,
starts_with("winner_minus_loser")
) |>
pivot_longer(
cols = starts_with("winner_minus_loser"),
names_to = "metric",
values_to = "difference"
) |>
mutate(
metric = recode(
metric,
winner_minus_loser_kill_rate =
"Kill-rate difference",
winner_minus_loser_error_rate =
"Error-rate difference",
winner_minus_loser_block_rate =
"Block-rate difference",
winner_minus_loser_aces =
"Aces-per-set difference",
winner_minus_loser_digs =
"Digs-per-set difference"
)
)Section 2.5.3: Winner-Loser Summary Statistics
Summary of typical winner-loser difference for each performance metric and reports how often winners recorded higher values.
Section 2.5.4: Winner-Loser Boxplots
Comparison of the distributions of winner-loser performance differences across competitions.
ggplot(
winner_loser_long,
aes(
x = competition,
y = difference
)
) +
geom_hline(
yintercept = 0,
linetype = "dashed"
) +
geom_boxplot(
outlier.alpha = 0.15
) +
facet_wrap(
vars(metric),
scales = "free_y",
ncol = 5
) +
labs(
title = "Within-Match Performance Differences",
subtitle = "Positive values indicate that the match winner recorded the higher value",
x = NULL,
y = "Winner minus loser"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 20, face = "bold"),
plot.subtitle = element_text(size = 14),
strip.text = element_text(size = 10),
axis.title = element_text(size = 16),
axis.text.x = element_text(angle = 20, hjust = 1, size = 11),
axis.text.y = element_text(size = 12)
)## Warning: Removed 1027 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
Section 2.6: Player Coverage and Playing Time
Section 2.6.1: Position-Label Summary
Summary of the original position labels in the player-season data and identifies missing or inconsistently coded positions.
position_label_summary <- eda_player_season |>
mutate(
position_label = case_when(
is.na(Pos) | str_squish(Pos) == "" ~
"Missing",
TRUE ~ str_squish(Pos)
)
) |>
count(
competition,
position_label,
sort = TRUE,
name = "player_seasons"
)
position_label_summary## # A tibble: 20 × 3
## competition position_label player_seasons
## <chr> <chr> <int>
## 1 Women's Division I OH 10024
## 2 Women's Division I MB 6862
## 3 Women's Division I S 5130
## 4 Women's Division I L/DS 3158
## 5 Women's Division I DS 2101
## 6 Men's Division I OH 1509
## 7 Women's Division I L 1405
## 8 Women's Division I RS 1390
## 9 Women's Division I OPP 918
## 10 Men's Division I MB 905
## 11 Men's Division I S 661
## 12 Women's Division I MH 613
## 13 Men's Division I L 486
## 14 Men's Division I OPP 405
## 15 Men's Division I Missing 223
## 16 Men's Division I L/DS 106
## 17 Men's Division I MH 86
## 18 Men's Division I RS 84
## 19 Women's Division I Missing 72
## 20 Men's Division I DB 60Section 2.6.2: Clean Player Playing-Time Data
Converts sets played to a numeric variable and removes records with missing or invalid playing-time values.
player_playing_time <- eda_player_season |>
mutate(
sets_played = as.numeric(S)
) |>
filter(
!is.na(sets_played),
sets_played >= 0
)Section 2.6.3: Playing-Time Summary
Summarizes the distribution of player participation and the number of players reaching the proposed minimum playing-time threshold of 20 sets.
playing_time_summary <- player_playing_time |>
group_by(competition, Season) |>
summarise(
player_seasons = n(),
median_sets = median(sets_played),
first_quartile = quantile(
sets_played,
0.25
),
third_quartile = quantile(
sets_played,
0.75
),
players_20_sets = sum(
sets_played >= 20
),
proportion_20_sets = mean(
sets_played >= 20
),
.groups = "drop"
)
playing_time_summary## # A tibble: 11 × 8
## competition Season player_seasons median_sets first_quartile third_quartile
## <chr> <chr> <int> <dbl> <dbl> <dbl>
## 1 Men's Divisi… 2020 644 39 16 63
## 2 Men's Divisi… 2021 857 53 19 81
## 3 Men's Divisi… 2022 941 52 16 86
## 4 Men's Divisi… 2023 988 50 17 81
## 5 Men's Divisi… 2024 1095 48 16 80
## 6 Women's Divi… 2020 4545 40 19 59
## 7 Women's Divi… 2021 5400 72 31 100
## 8 Women's Divi… 2022 5410 78 31 104
## 9 Women's Divi… 2023 5361 78 32 104
## 10 Women's Divi… 2024 5472 79 33 104
## 11 Women's Divi… 2025 5485 79 34 103
## # ℹ 2 more variables: players_20_sets <int>, proportion_20_sets <dbl>Section 2.6.4: Sets-Played Histogram
Shows the distribution of player-season playing time and marks the 20-set eligibility threshold for the PEI later on.
ggplot(
player_playing_time,
aes(x = sets_played)
) +
geom_histogram(
bins = 40
) +
geom_vline(
xintercept = 20,
linetype = "dashed"
) +
facet_wrap(
vars(competition),
scales = "free_y"
) +
labs(
title = "Distribution of Player Sets Played",
subtitle = "The dashed line marks the proposed 20-set PEI threshold",
x = "Sets played",
y = "Player-seasons"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 20, face = "bold"),
plot.subtitle = element_text(size = 14),
strip.text = element_text(size = 10),
axis.title = element_text(size = 16),
axis.text.x = element_text(size = 11),
axis.text.y = element_text(size = 12)
)
Section 2.6.5: Qualification Sensitivity
Shows how many player-seasons would qualify under alternative set thresholds. This shows the trade off between retaining more players and excluding athletes with very limited playing time.
qualification_sensitivity <- player_playing_time |>
group_by(competition, Season) |>
summarise(
total_players = n(),
qualify_20 = sum(sets_played >= 20),
qualify_30 = sum(sets_played >= 30),
qualify_40 = sum(sets_played >= 40),
qualify_50 = sum(sets_played >= 50),
.groups = "drop"
)
qualification_sensitivity## # A tibble: 11 × 7
## competition Season total_players qualify_20 qualify_30 qualify_40 qualify_50
## <chr> <chr> <int> <int> <int> <int> <int>
## 1 Men's Divis… 2020 644 449 382 316 254
## 2 Men's Divis… 2021 857 638 586 522 449
## 3 Men's Divis… 2022 941 680 611 548 487
## 4 Men's Divis… 2023 988 713 629 573 500
## 5 Men's Divis… 2024 1095 790 694 608 540
## 6 Women's Div… 2020 4545 3380 2863 2325 1722
## 7 Women's Div… 2021 5400 4415 4115 3788 3482
## 8 Women's Div… 2022 5410 4390 4113 3824 3547
## 9 Women's Div… 2023 5361 4421 4099 3797 3533
## 10 Women's Div… 2024 5472 4530 4204 3920 3646
## 11 Women's Div… 2025 5485 4570 4252 3955 3655Section 2.7: EDA Findings and Modeling Implications
The exploratory analysis show significant differences in data volume between women’s and men’s competitions while showing recurring patterns in match length, team performance and player participation. Rate-based measures were used to improve comparability across matches of varying lengths. The analysis also supported the use of a 20-set minimum threshold for the PEI calculation. The PEI eventually required participation in at least 40% of the team’s reported sets for the season. The winner-loser performance comparison supported the inclusion of attacking, error, blocking, serving and defensive metrics in match prediction section.
Section 3: Player Efficiency Index
Section 3.1: Purpose and Interpretation
The Player Efficiency Index (PEI) measures how exceptional a player’s statistical production was relative to other qualified players at the same position in the same season. Component weights are estimated separately by competition and position using principal component analysis of qualified player-season performance. The PEI is retrospective and descriptive; it is not a causal estimate of wins produced by an individual player.
A PEI of 50 represents the average qualified player within a season-position group. A score of 60 is approximately one standard deviation above that average, and a score of 40 is approximately one standard deviation below it.
Section 3.2: Position Classification and PEI Components
Section 3.2.1: Position-Classification and Z-Score Functions
Helper functions which convert inconsistent position labels into five standardized position groups and safely standardize player metrics using z-scores. The ‘Other’ category contains missing or unrecognized position labels but excludes those records from position-specific PEI scoring.
classify_position <- function(x) {
position <- x |>
coalesce("") |>
str_to_upper() |>
str_replace_all("[._-]", " ") |>
str_squish()
case_when(
position == "" ~ "Other",
str_detect(position, "LIBERO|(^|[/ ])L($|[/ ])|\\bDS\\b|\\bDB\\b|DEFENSIVE") ~ "L_DS",
str_detect(position, "SETTER|(^|[/ ])S($|[/ ])") ~ "S",
str_detect(position, "MIDDLE|\\bMB\\b|\\bMH\\b") ~ "MB",
str_detect(position, "OPPOSITE|RIGHT SIDE|\\bRS\\b|\\bOPP\\b") ~ "OPP",
str_detect(position, "OUTSIDE|LEFT SIDE|\\bOH\\b|\\bLS\\b") ~ "OH",
TRUE ~ "Other"
)
}
safe_z <- function(x) {
valid <- is.finite(x)
result <- rep(NA_real_, length(x))
if (sum(valid) == 0) {
return(result)
}
if (sum(valid) == 1) {
result[valid] <- 0
return(result)
}
x_sd <- sd(x[valid], na.rm = TRUE)
if (!is.finite(x_sd) || x_sd == 0) {
result[valid] <- 0
} else {
result[valid] <-
(x[valid] - mean(x[valid], na.rm = TRUE)) / x_sd
}
result
}Section 3.2.2: Position Audit
Shows how the original position labels were assigned to the standardized position groups in each competition.
Section 3.2.3: PEI Feature Specification
Defines performance metrics included in the PEI calculation for each position which allows players to be evaluated according to the responsibilities of their role.
pei_feature_spec <- tribble(
~pos_group, ~metric,
"OH", "net_attack_per_set",
"OH", "net_serve_per_set",
"OH", "digs_per_set",
"OH", "blocks_per_set",
"OPP", "net_attack_per_set",
"OPP", "net_serve_per_set",
"OPP", "digs_per_set",
"OPP", "blocks_per_set",
"MB", "net_attack_per_set",
"MB", "net_serve_per_set",
"MB", "digs_per_set",
"MB", "blocks_per_set",
"S", "assists_per_set",
"S", "net_serve_per_set",
"S", "digs_per_set",
"S", "blocks_per_set",
"L_DS", "net_serve_per_set",
"L_DS", "digs_per_set"
) |>
mutate(
feature = paste(
pos_group,
metric,
sep = "__"
)
)Section 3.3: Prepare Player-Season Components
Player-season statistics were converted into a more comparable per-set version, joined to team participation totals and then filtered using the minimum set requirement.
Section 3.3.1: Component-Preparation Function
Cleans the player statistics, calculates position-relevant per-set production, measures each player’s share of team set and identifies players who satisfy the PEI eligibility threshold.
prepare_player_season_components <- function(
player_season,
team_season,
minimum_sets = 20,
minimum_team_share = 0.40
) {
team_sets_lookup <- team_season |>
transmute(
canonical_season = as.character(Season),
team_clean = clean_team_name(Team),
team_sets = suppressWarnings(as.numeric(S))
) |>
filter(
!is.na(team_sets),
team_sets > 0
) |>
group_by(
canonical_season,
team_clean
) |>
summarise(
team_sets = max(team_sets, na.rm = TRUE),
.groups = "drop"
)
player_season |>
mutate(
canonical_season =
as.character(Season),
original_player_season =
as.character(Season),
team_clean =
clean_team_name(Team),
pos_group =
classify_position(Pos),
player_sets = suppressWarnings(
as.numeric(S)
),
kills_num = coalesce(
suppressWarnings(as.numeric(Kills)),
0
),
attack_errors_num = coalesce(
suppressWarnings(as.numeric(Errors)),
0
),
assists_num = coalesce(
suppressWarnings(as.numeric(Assists)),
0
),
aces_num = coalesce(
suppressWarnings(as.numeric(Aces)),
0
),
service_errors_num = coalesce(
suppressWarnings(as.numeric(SErr)),
0
),
digs_num = coalesce(
suppressWarnings(as.numeric(Digs)),
0
),
block_solos_num = coalesce(
suppressWarnings(
as.numeric(`Block Solos`)
),
0
),
block_assists_num = coalesce(
suppressWarnings(
as.numeric(`Block Assists`)
),
0
)
) |>
left_join(
team_sets_lookup,
by = c(
"canonical_season",
"team_clean"
)
) |>
mutate(
team_set_share_raw = safe_divide(
player_sets,
team_sets
),
team_set_share = if_else(
team_set_share_raw >= 0 &
team_set_share_raw <= 1,
team_set_share_raw,
NA_real_
),
qualifies =
!is.na(player_sets) &
player_sets >= minimum_sets &
!is.na(team_set_share) &
team_set_share >= minimum_team_share,
net_attack_per_set = safe_divide(
kills_num - attack_errors_num,
player_sets
),
assists_per_set = safe_divide(
assists_num,
player_sets
),
net_serve_per_set = safe_divide(
aces_num - service_errors_num,
player_sets
),
digs_per_set = safe_divide(
digs_num,
player_sets
),
blocks_per_set = safe_divide(
block_solos_num +
0.5 * block_assists_num,
player_sets
)
)
}Section 3.3.2: Applying the Function to Both Competitions
Application of the function in 3.3.1 to women’s and men’s player-season data separately.
wvb_player_components <-
prepare_player_season_components(
player_season = wvb_player_season,
team_season = wvb_team_season,
minimum_sets = 20,
minimum_team_share = 0.40
)
mvb_player_components <-
prepare_player_season_components(
player_season = mvb_player_season,
team_season = mvb_team_season,
minimum_sets = 20,
minimum_team_share = 0.40
)Section 3.3.3: Eligibility and Participation Validation
Summary of player-season coverage, team-set matching, eligibility and unclassified positions for each competition.
player_component_summary <- bind_rows(
wvb_player_components |>
summarise(
competition = "Women's Division I",
player_seasons = n(),
invalid_raw_shares = sum(
!is.na(team_set_share_raw) &
(
team_set_share_raw < 0 |
team_set_share_raw > 1
)
),
valid_team_shares =
sum(!is.na(team_set_share)),
qualified_players =
sum(qualifies, na.rm = TRUE),
maximum_valid_share =
max(team_set_share, na.rm = TRUE)
),
mvb_player_components |>
summarise(
competition = "Men's Division I",
player_seasons = n(),
invalid_raw_shares = sum(
!is.na(team_set_share_raw) &
(
team_set_share_raw < 0 |
team_set_share_raw > 1
)
),
valid_team_shares =
sum(!is.na(team_set_share)),
qualified_players =
sum(qualifies, na.rm = TRUE),
maximum_valid_share =
max(team_set_share, na.rm = TRUE)
)
)
player_component_summary## # A tibble: 2 × 6
## competition player_seasons invalid_raw_shares valid_team_shares
## <chr> <int> <int> <int>
## 1 Women's Division I 31673 0 31673
## 2 Men's Division I 4525 8 4517
## # ℹ 2 more variables: qualified_players <int>, maximum_valid_share <dbl>Section 3.4: Estimate Position-Specific PCA Weights
Principal Component Analysis (PCA) is estimated separately by competition and position so that the relative importance of each metric is determined from the observed player data.
Section 3.4.1: PCA Weight-Estimation Function
Function which standardizes each metric within season and position, estimating the first principal component and converts the absolute component loading into non-negative weights that sum to one.
estimate_position_pca_weights <- function(
player_components,
feature_spec,
minimum_players = 20
) {
metric_spec <- feature_spec |>
distinct(
pos_group,
metric,
feature
)
eligible_players <- player_components |>
filter(
qualifies,
pos_group != "Other",
!is.na(canonical_season)
) |>
mutate(
pca_row_id = row_number()
)
component_long <- eligible_players |>
pivot_longer(
cols = c(
net_attack_per_set,
assists_per_set,
net_serve_per_set,
digs_per_set,
blocks_per_set
),
names_to = "metric",
values_to = "component_value"
) |>
inner_join(
metric_spec,
by = c(
"pos_group",
"metric"
)
) |>
group_by(
canonical_season,
pos_group,
metric
) |>
mutate(
component_z = safe_z(
component_value
)
) |>
ungroup()
position_groups <- sort(
unique(metric_spec$pos_group)
)
pca_tables <- lapply(
position_groups,
function(position_name) {
requested_spec <- metric_spec |>
filter(
pos_group == position_name
)
requested_metrics <- requested_spec$metric
position_wide <- component_long |>
filter(
pos_group == position_name
) |>
select(
pca_row_id,
metric,
component_z
) |>
pivot_wider(
names_from = metric,
values_from = component_z
)
# Add any expected metric that is entirely absent.
missing_metrics <- setdiff(
requested_metrics,
names(position_wide)
)
if (length(missing_metrics) > 0) {
for (metric_name in missing_metrics) {
position_wide[[metric_name]] <-
NA_real_
}
}
metric_variation <- vapply(
position_wide[
requested_metrics
],
function(x) {
if (sum(is.finite(x)) < 2) {
return(NA_real_)
}
sd(
x,
na.rm = TRUE
)
},
numeric(1)
)
retained_metrics <- names(
metric_variation[
is.finite(metric_variation) &
metric_variation > 0
]
)
# PCA requires at least two variables.
if (length(retained_metrics) >= 2) {
complete_rows <- position_wide |>
filter(
if_all(
all_of(retained_metrics),
is.finite
)
)
} else {
complete_rows <- position_wide[
0,
,
drop = FALSE
]
}
use_pca <-
length(retained_metrics) >= 2 &&
nrow(complete_rows) >= minimum_players
if (!use_pca) {
fallback_metrics <- retained_metrics
if (length(fallback_metrics) == 0) {
fallback_metrics <-
requested_metrics
}
return(
requested_spec |>
mutate(
loading = NA_real_,
oriented_loading = NA_real_,
weight = if_else(
metric %in% fallback_metrics,
1 / length(fallback_metrics),
0
),
pc1_variance_explained =
NA_real_,
players_used =
nrow(complete_rows),
weighting_method =
"Equal-weight fallback"
)
)
}
pca_matrix <- as.matrix(
complete_rows[
retained_metrics
]
)
# Components were already standardized within
# season and position.
pca_fit <- prcomp(
pca_matrix,
center = FALSE,
scale. = FALSE
)
raw_loadings <- pca_fit$rotation[
retained_metrics,
1
]
# The mathematical sign of a PCA component is arbitrary.
orientation <- if (
sum(raw_loadings, na.rm = TRUE) < 0
) {
-1
} else {
1
}
oriented_loadings <-
raw_loadings * orientation
loading_magnitudes <-
abs(oriented_loadings)
pca_weights <-
loading_magnitudes /
sum(loading_magnitudes)
variance_explained <- (
pca_fit$sdev^2 /
sum(pca_fit$sdev^2)
)[1]
loading_table <- tibble(
metric = retained_metrics,
loading =
as.numeric(raw_loadings),
oriented_loading =
as.numeric(oriented_loadings),
weight =
as.numeric(pca_weights)
)
requested_spec |>
left_join(
loading_table,
by = "metric"
) |>
mutate(
weight = coalesce(
weight,
0
),
pc1_variance_explained =
variance_explained,
players_used =
nrow(complete_rows),
weighting_method =
"PCA-derived"
)
}
)
weight_table <- bind_rows(
pca_tables
) |>
group_by(
pos_group
) |>
mutate(
weight_total = sum(
weight,
na.rm = TRUE
),
weight = if_else(
weight_total > 0,
weight / weight_total,
1 / n()
)
) |>
ungroup() |>
select(
pos_group,
metric,
feature,
loading,
oriented_loading,
weight,
pc1_variance_explained,
players_used,
weighting_method
)
list(
weights = weight_table,
component_long = component_long
)
}Section 3.4.2: Estimate Women’s and Men’s PCA Models
Estimation of separate position-specific PCA weights for women’s and men’s volleyball.
wvb_pca_fit <- estimate_position_pca_weights(
player_components =
wvb_player_components,
feature_spec =
pei_feature_spec,
minimum_players = 20
)
mvb_pca_fit <- estimate_position_pca_weights(
player_components =
mvb_player_components,
feature_spec =
pei_feature_spec,
minimum_players = 20
)Section 3.4.3: Combine the Weight Tables
Combines the PCA results into one table for validation and presentation.
estimated_pei_weights <- bind_rows(
wvb_pca_fit$weights |>
mutate(
competition = "Women's Division I"
),
mvb_pca_fit$weights |>
mutate(
competition = "Men's Division I"
)
) |>
select(
competition,
everything()
) |>
arrange(
competition,
pos_group,
desc(weight)
)Section 3.4.4: Display the Estimated Weights
Displays each metrics normalized PEI weight, variance explained by the first principal component and whether PCA or the equal-weight fallback was used.
estimated_pei_weights |>
select(
competition,
pos_group,
metric,
weight,
pc1_variance_explained,
weighting_method
) |>
mutate(
weight = round(weight, 3),
pc1_variance_explained =
round(pc1_variance_explained, 3)
)## # A tibble: 36 × 6
## competition pos_group metric weight pc1_variance_explained weighting_method
## <chr> <chr> <chr> <dbl> <dbl> <chr>
## 1 Men's Divisi… L_DS net_s… 0.5 0.63 PCA-derived
## 2 Men's Divisi… L_DS digs_… 0.5 0.63 PCA-derived
## 3 Men's Divisi… MB block… 0.429 0.332 PCA-derived
## 4 Men's Divisi… MB net_a… 0.297 0.332 PCA-derived
## 5 Men's Divisi… MB digs_… 0.187 0.332 PCA-derived
## 6 Men's Divisi… MB net_s… 0.087 0.332 PCA-derived
## 7 Men's Divisi… OH net_a… 0.307 0.528 PCA-derived
## 8 Men's Divisi… OH block… 0.267 0.528 PCA-derived
## 9 Men's Divisi… OH net_s… 0.221 0.528 PCA-derived
## 10 Men's Divisi… OH digs_… 0.205 0.528 PCA-derived
## # ℹ 26 more rowsSection 3.5: Calculate PEI scores
The estimated PCA weights were applied to the standardized player metrics in order to get the final position-adjusted PEI scores and rankings.
Section 3.5.1: PEI Scoring Function
Function which standardizes each component within season and position, applies the PCA derived weights and converts the resulting score to a PEI scale with mean 50 and standard deviation of 10. It also calculates each player’s position rank and percentile within the same season-position group.
score_player_season_pei <- function(
player_components,
weight_table
) {
eligible_players <- player_components |>
filter(
qualifies,
pos_group != "Other",
!is.na(canonical_season)
) |>
select(
-any_of("Season")
) |>
mutate(
player_row_id = row_number()
)
component_scores <- eligible_players |>
pivot_longer(
cols = c(
net_attack_per_set,
assists_per_set,
net_serve_per_set,
digs_per_set,
blocks_per_set
),
names_to = "metric",
values_to = "component_value"
) |>
inner_join(
weight_table |>
select(
pos_group,
metric,
weight
),
by = c(
"pos_group",
"metric"
)
) |>
group_by(
canonical_season,
pos_group,
metric
) |>
mutate(
component_z = safe_z(
component_value
)
) |>
ungroup() |>
mutate(
weighted_component =
weight * component_z
)
score_summary <- component_scores |>
group_by(
player_row_id
) |>
summarise(
weighted_component_sum = sum(
weighted_component,
na.rm = TRUE
),
available_weight = sum(
weight[
is.finite(component_z)
],
na.rm = TRUE
),
pei_raw = if_else(
available_weight > 0,
weighted_component_sum /
available_weight,
NA_real_
),
.groups = "drop"
)
rankings <- eligible_players |>
left_join(
score_summary,
by = "player_row_id"
) |>
group_by(
canonical_season,
pos_group
) |>
mutate(
PEI = round(
50 + 10 * safe_z(
pei_raw
),
1
),
position_rank = min_rank(
desc(pei_raw)
),
position_percentile = if (
dplyr::n() > 1
) {
round(
100 * percent_rank(
pei_raw
),
1
)
} else {
rep(
50,
dplyr::n()
)
},
players_at_position =
dplyr::n()
) |>
ungroup() |>
rename(
Season = canonical_season
) |>
select(
-player_row_id
)
list(
rankings = rankings,
component_scores =
component_scores
)
}Section 3.5.2: Calculate Women’s and Men’s PEI
Calculates the final women’s and men’s PEI results and saves the player rankings and component-level contributions.
wvb_pei_results <- score_player_season_pei(
player_components = wvb_player_components,
weight_table = wvb_pca_fit$weights
)
mvb_pei_results <- score_player_season_pei(
player_components = mvb_player_components,
weight_table = mvb_pca_fit$weights
)
wvb_player_pei <- wvb_pei_results$rankings
mvb_player_pei <- mvb_pei_results$rankingsSection 3.6: PEI Rankings
Presents the highest-rated players, summarizes PEI coverage and examines which statistical components contributed to selected player scores.
Section 3.6.1: Identify Latest Available Season
Identifies the most recent season with PEI results for each competition.
latest_wvb_season <- wvb_player_pei |>
distinct(Season) |>
mutate(
season_order = readr::parse_number(
as.character(Season)
)
) |>
filter(
!is.na(season_order)
) |>
slice_max(
season_order,
n = 1,
with_ties = FALSE
) |>
pull(Season)
latest_mvb_season <- mvb_player_pei |>
distinct(Season) |>
mutate(
season_order = Season
) |>
filter(
!is.na(season_order)
) |>
slice_max(
season_order,
n = 1,
with_ties = FALSE
) |>
pull(Season)Section 3.6.2: Latest-Season Leaders
Identifies the five highest rated women’s players at each position in the latest season.
Identifies the five highest rated men’s players at each position in the latest season.
Section 3.6.3: PEI Player Coverage
Reports how many qualified players received PEI scores in each season, competition and position group.
Section 3.6.4: Top 5 Players By Season And Position
This table retains the five highest-ranked women’s players within each season-position group.
This table retains the five highest-ranked men’s players within each season-position group.
Section 3.6.5: Overall Leaders Within Each Season
Identifies the women’s players with the highest position-adjusted PEI scores in each season.
Identifies the men’s players with the highest position-adjusted PEI scores in each season.
Since PEI is standardized within each position, the tables identify players who are most exceptional relative to peers in that position rather than showing that one position is more valuable than another.
Section 3.6.6: Latest-Season Leaders Plot
Displays the three highest-rated women’s players at each position in the latest season.
wvb_latest_plot_data <- wvb_player_pei |>
filter(
Season == "2024",
is.finite(PEI)
) |>
group_by(
pos_group
) |>
slice_max(
PEI,
n = 3,
with_ties = FALSE
) |>
ungroup() |>
mutate(
player_label = paste(
Player,
Team,
sep = " — "
),
player_label = reorder(
player_label,
PEI
)
)
ggplot(
wvb_latest_plot_data,
aes(
x = PEI,
y = player_label
)
) +
geom_col() +
geom_text(
aes(
x = PEI / 2,
label = round(PEI, 1)
),
color = "white",
size = 3.5
) +
facet_wrap(
vars(pos_group),
scales = "free_y",
ncol = 1
) +
labs(
title = paste(
"Women's PEI Leaders:",
"2024"
),
subtitle =
"Top three qualified players within each position",
x = "Player Efficiency Index",
y = NULL
) +
theme_minimal() +
theme(
plot.title = element_text(size = 20, face = "bold"),
plot.subtitle = element_text(size = 14),
strip.text = element_text(size = 12),
axis.title = element_text(size = 16),
axis.text.x = element_text(size = 11),
axis.text.y = element_text(size = 12)
)
Displays the three highest-rated men’s players at each position in the latest season.
mvb_latest_plot_data <- mvb_player_pei |>
filter(
Season == "2023",
is.finite(PEI)
) |>
group_by(
pos_group
) |>
slice_max(
PEI,
n = 3,
with_ties = FALSE
) |>
ungroup() |>
mutate(
player_label = paste(
Player,
Team,
sep = " — "
),
player_label = reorder(
player_label,
PEI
)
)
ggplot(
mvb_latest_plot_data,
aes(
x = PEI,
y = player_label
)
) +
geom_col() +
geom_text(
aes(
x = PEI / 2,
label = round(PEI, 1)
),
color = "white",
size = 3.5
) +
facet_wrap(
vars(pos_group),
scales = "free_y",
ncol = 1
) +
labs(
title = paste(
"Men's PEI Leaders:",
"2023"
),
subtitle =
"Top three qualified players within each position",
x = "Player Efficiency Index",
y = NULL
) +
theme_minimal() +
theme(
plot.title = element_text(size = 20, face = "bold"),
plot.subtitle = element_text(size = 14),
strip.text = element_text(size = 12),
axis.title = element_text(size = 16),
axis.text.x = element_text(size = 11),
axis.text.y = element_text(size = 12)
)
Section 3.6.7: Women’s Player Component Profiles
Identifies the latest-season women’s PEI leaders whose component-level contributions will be examined.
Shows how each standardized metric and its PCA weight contributed to the PEI scores of the selected women’s players.
Section 3.6.8: Men’s Player Component Profiles
Identifies the latest-season men’s PEI leaders whose component-level contributions will be examined.
Shows how each standardized metric and its PCA weight contributed to the PEI scores of the selected men’s players.
Section 3.7: Validation and Sensitivity Analysis
These checks evaluate if the PEI is constructed correctly and whether the player rankings are sensitive to the selected weighting model.
Section 3.7.1: Verify That Weights Sum To One
Confirms that the component weights sum to one within every competition-position group and reports whether PCA or the fallback method was used.
pei_weight_validation <- estimated_pei_weights |>
group_by(
competition,
pos_group
) |>
summarise(
weight_sum = sum(
weight,
na.rm = TRUE
),
pc1_variance_explained =
first(pc1_variance_explained),
players_used =
first(players_used),
weighting_method =
first(weighting_method),
.groups = "drop"
) |>
mutate(
valid_weight_sum =
abs(weight_sum - 1) < 1e-8
)
pei_weight_validation## # A tibble: 10 × 7
## competition pos_group weight_sum pc1_variance_explained players_used
## <chr> <chr> <dbl> <dbl> <int>
## 1 Men's Division I L_DS 1 0.630 427
## 2 Men's Division I MB 1 0.332 563
## 3 Men's Division I OH 1 0.528 899
## 4 Men's Division I OPP 1 0.525 254
## 5 Men's Division I S 1 0.528 391
## 6 Women's Division I L_DS 1 0.552 4847
## 7 Women's Division I MB 1 0.342 4928
## 8 Women's Division I OH 1 0.379 7024
## 9 Women's Division I OPP 1 0.419 1495
## 10 Women's Division I S 1 0.518 3615
## # ℹ 2 more variables: weighting_method <chr>, valid_weight_sum <lgl>Section 3.7.2: Variance Explained by PC1
Shows how much of the variation in each position’s standardized player statistics is summarized by the first principal component. Higher values indicate a stronger common statistical performance dimension while lower values indicate that the position’s contributions are more multidimensional.
pca_variance_summary <- estimated_pei_weights |>
distinct(
competition,
pos_group,
pc1_variance_explained,
players_used,
weighting_method
) |>
mutate(
pc1_variance_percent =
scales::percent(
pc1_variance_explained,
accuracy = 0.1
)
) |>
arrange(
competition,
desc(pc1_variance_explained)
)
pca_variance_summary## # A tibble: 10 × 6
## competition pos_group pc1_variance_explained players_used weighting_method
## <chr> <chr> <dbl> <int> <chr>
## 1 Men's Divisio… L_DS 0.630 427 PCA-derived
## 2 Men's Divisio… S 0.528 391 PCA-derived
## 3 Men's Divisio… OH 0.528 899 PCA-derived
## 4 Men's Divisio… OPP 0.525 254 PCA-derived
## 5 Men's Divisio… MB 0.332 563 PCA-derived
## 6 Women's Divis… L_DS 0.552 4847 PCA-derived
## 7 Women's Divis… S 0.518 3615 PCA-derived
## 8 Women's Divis… OPP 0.419 1495 PCA-derived
## 9 Women's Divis… OH 0.379 7024 PCA-derived
## 10 Women's Divis… MB 0.342 4928 PCA-derived
## # ℹ 1 more variable: pc1_variance_percent <chr>Section 3.7.3: Check PEI Distributions
Verifies that the PEI scores are centered near 50 with a standard deviation near 10 within each season-position group.
pei_distribution_check <- bind_rows(
wvb_player_pei |>
group_by(
Season,
pos_group
) |>
summarise(
players = n(),
mean_pei = mean(
PEI,
na.rm = TRUE
),
sd_pei = sd(
PEI,
na.rm = TRUE
),
minimum_pei = min(
PEI,
na.rm = TRUE
),
maximum_pei = max(
PEI,
na.rm = TRUE
),
.groups = "drop"
) |>
mutate(
competition =
"Women's Division I"
),
mvb_player_pei |>
group_by(
Season,
pos_group
) |>
summarise(
players = n(),
mean_pei = mean(
PEI,
na.rm = TRUE
),
sd_pei = sd(
PEI,
na.rm = TRUE
),
minimum_pei = min(
PEI,
na.rm = TRUE
),
maximum_pei = max(
PEI,
na.rm = TRUE
),
.groups = "drop"
) |>
mutate(
competition =
"Men's Division I"
)
) |>
select(
competition,
everything()
)
pei_distribution_check## # A tibble: 55 × 8
## competition Season pos_group players mean_pei sd_pei minimum_pei maximum_pei
## <chr> <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Women's Div… 2020 L_DS 709 50.0 10.0 23.3 82.6
## 2 Women's Div… 2020 MB 726 50.0 10.0 16.3 91.7
## 3 Women's Div… 2020 OH 1020 50.0 10.0 23.6 81.2
## 4 Women's Div… 2020 OPP 199 50.0 10.0 19.3 85
## 5 Women's Div… 2020 S 508 50 10.0 26.8 74.3
## 6 Women's Div… 2021 L_DS 816 50.0 10.0 11.5 79.1
## 7 Women's Div… 2021 MB 846 50 10.0 20 95.6
## 8 Women's Div… 2021 OH 1189 50.0 10.0 25 84.9
## 9 Women's Div… 2021 OPP 242 50.0 10.0 25.9 101.
## 10 Women's Div… 2021 S 604 50.0 10.0 27.7 76
## # ℹ 45 more rowsSection 3.7.4: Equal-Weight Benchmark
Creates an alternative scoring system that assigns every included metric equal importance within each position. This essentially recalculates the rankings using equal weights (instead of PCA).
Section 3.7.5: Compare PCA and Equal-Weight Rankings
Helper function safely calculates Spearman rank correlations while returning missing values for groups that are either too small or have no variation in score.
Table which compares PCA and equal-weight PEI rankings using Spearman correlations and average absolute score differences.
Section 3.7.6: Audit PCA Loading Directions
Identifies metrics whose first-component loading remain negative after orienting the PCA solution which could indicate trade offs among dimensions of player production.
Section 3.7.7: Structural Checks
Automated checks confirm whether PEI values are finite, ranks are valid and player participation shares are non-negative.
stopifnot(
all(
is.na(wvb_player_pei$PEI) |
is.finite(wvb_player_pei$PEI)
),
all(
is.na(mvb_player_pei$PEI) |
is.finite(mvb_player_pei$PEI)
),
all(
wvb_player_pei$position_rank[
!is.na(wvb_player_pei$position_rank)
] >= 1
),
all(
mvb_player_pei$position_rank[
!is.na(mvb_player_pei$position_rank)
] >= 1
),
all(
wvb_player_pei$team_set_share[
!is.na(wvb_player_pei$team_set_share)
] >= 0 &
wvb_player_pei$team_set_share[
!is.na(wvb_player_pei$team_set_share)
] <= 1
),
all(
mvb_player_pei$team_set_share[
!is.na(mvb_player_pei$team_set_share)
] >= 0 &
mvb_player_pei$team_set_share[
!is.na(mvb_player_pei$team_set_share)
] <= 1
)
)Section 3.8: Limitations
The Player Efficiency Index is a purely statistical measure of production instead of wins generated by each player. Principal Component Analysis identified dominant patterns of variation among the selected metrics but doesn’t determine which one(s) cause teams to win. The index and calculations within are limited by the box-score variables, position-label quality and the lack of rotation, lineup, opponent-strength and additional contextual information. PEI should be interpreted as a reproducible, position-adjusted ranking of player production instead of a complete measure of a player’s value.
Section 4: Match Outcome Prediction Model
Development of a logistic regression model that estimates the probability that Team A will win a match. Predictors are calculated from matches completed before the target match to prevent incoming information from unintentionally entering the model.
Section 4.1: Purpose and Model Design
Match outcomes are binary which means that a logistic regression can be an appropriate and interpretable model choice. Separate logistic models were estimated from women’s and men’s volleyball using difference in the team’s prior win percentage, attacking, error, blocking, serving and defensive performance. Team A and Team B are assigned alphabetically. The most recent season is used as a test set while the other years are used as the training set for the model. Due to the limitations of the data, the goal of the model was to create a useful baseline forecast model instead of a comprehensive prediction system which incorporates every possible piece of contextual information.
Section 4.2: Construct Pre-Match Team Features
Calculation of each team’s historical performance entering every match, using only matches played on earlier dates within the same source year.
Section 4.2.1: Prior Cumulative-Mean Helper Function
Helper function which calculates an expanding historical average while excluding current date results from the value used for prediction. This prevents matches from the target date being used to predict one another.
prior_cumulative_mean <- function(value_sum, value_count) {
prior_sum <- dplyr::lag(
cumsum(value_sum),
default = 0
)
prior_count <- dplyr::lag(
cumsum(value_count),
default = 0
)
dplyr::if_else(
prior_count > 0,
prior_sum / prior_count,
NA_real_
)
}Section 4.2.2: Pre-Match Feature Construction Function
Summarizes each team’s prior match count, win percentage, kill rate, error rate, block rate, aces per set, and digs per set before every match date. Historical averages reset at the beginning each season so every team begins each new season without carrying statistics forward from the prior season.
build_prematch_team_features <- function(team_rates) {
daily_team_results <- team_rates |>
group_by(
source_year,
Season,
team_clean,
match_date
) |>
summarise(
matches_on_date = n(),
win_sum = sum(win, na.rm = TRUE),
win_count = sum(is.finite(win)),
kill_sum = sum(
kill_rate_pct,
na.rm = TRUE
),
kill_count = sum(
is.finite(kill_rate_pct)
),
error_sum = sum(
error_rate_pct,
na.rm = TRUE
),
error_count = sum(
is.finite(error_rate_pct)
),
block_sum = sum(
block_rate_pct,
na.rm = TRUE
),
block_count = sum(
is.finite(block_rate_pct)
),
ace_sum = sum(
aces_per_set,
na.rm = TRUE
),
ace_count = sum(
is.finite(aces_per_set)
),
dig_sum = sum(
digs_per_set,
na.rm = TRUE
),
dig_count = sum(
is.finite(digs_per_set)
),
.groups = "drop"
) |>
arrange(
source_year,
team_clean,
match_date
) |>
group_by(
source_year,
team_clean
) |>
mutate(
prior_matches = lag(
cumsum(matches_on_date),
default = 0
),
prior_win_pct = prior_cumulative_mean(
win_sum,
win_count
),
prior_kill_rate = prior_cumulative_mean(
kill_sum,
kill_count
),
prior_error_rate = prior_cumulative_mean(
error_sum,
error_count
),
prior_block_rate = prior_cumulative_mean(
block_sum,
block_count
),
prior_aces_per_set = prior_cumulative_mean(
ace_sum,
ace_count
),
prior_digs_per_set = prior_cumulative_mean(
dig_sum,
dig_count
)
) |>
ungroup() |>
select(
source_year,
Season,
team_clean,
match_date,
starts_with("prior_")
)
team_rates |>
left_join(
daily_team_results,
by = c(
"source_year",
"Season",
"team_clean",
"match_date"
)
)
}Section 4.2.3: Apply Function to Women’s and Men’s
Application of the feature construction function to women’s and men’s data separately.
wvb_prematch_team <- build_prematch_team_features(
wvb_team_rates
)
mvb_prematch_team <- build_prematch_team_features(
mvb_team_rates
)Section 4.3: Construct One Prediction Row Per Match
Team-level historical features are combined into one row per match so that the model can compare the two participating teams directly.
Section 4.3.1: Prediction-Data Function
Function which joins Team A and Team B records, calculates Team A minus Team B differences for each predictor and removes matches in which either team has fewer than three prior matches. A positive predictor difference means Team A had the higher historical value while a negative difference means Team B had the higher value (the exception is error rate where a positive difference means Team A committed errors at a higher rate).
build_prediction_match_data <- function(
prematch_team,
minimum_prior_matches = 3
) {
team_a_rows <- prematch_team |>
filter(
team_clean == team_low
) |>
transmute(
match_id,
source_year = as.integer(source_year),
Season,
match_date,
team_a = team_low,
team_b = team_high,
team_a_win = win,
team_a_prior_matches =
prior_matches,
team_a_prior_win_pct =
prior_win_pct,
team_a_prior_kill_rate =
prior_kill_rate,
team_a_prior_error_rate =
prior_error_rate,
team_a_prior_block_rate =
prior_block_rate,
team_a_prior_aces_per_set =
prior_aces_per_set,
team_a_prior_digs_per_set =
prior_digs_per_set
)
team_b_rows <- prematch_team |>
filter(
team_clean == team_high
) |>
transmute(
match_id,
team_b_prior_matches =
prior_matches,
team_b_prior_win_pct =
prior_win_pct,
team_b_prior_kill_rate =
prior_kill_rate,
team_b_prior_error_rate =
prior_error_rate,
team_b_prior_block_rate =
prior_block_rate,
team_b_prior_aces_per_set =
prior_aces_per_set,
team_b_prior_digs_per_set =
prior_digs_per_set
)
team_a_rows |>
inner_join(
team_b_rows,
by = "match_id"
) |>
mutate(
diff_prior_win_pct =
team_a_prior_win_pct -
team_b_prior_win_pct,
diff_prior_kill_rate =
team_a_prior_kill_rate -
team_b_prior_kill_rate,
diff_prior_error_rate =
team_a_prior_error_rate -
team_b_prior_error_rate,
diff_prior_block_rate =
team_a_prior_block_rate -
team_b_prior_block_rate,
diff_prior_aces_per_set =
team_a_prior_aces_per_set -
team_b_prior_aces_per_set,
diff_prior_digs_per_set =
team_a_prior_digs_per_set -
team_b_prior_digs_per_set
) |>
filter(
team_a_prior_matches >=
minimum_prior_matches,
team_b_prior_matches >=
minimum_prior_matches
) |>
filter(
if_all(
starts_with("diff_prior_"),
~ is.finite(.x)
)
)
}Section 4.3.2: Construct Women’s and Men’s Prediction Datasets
Creates the final women’s and men’s modeling datasets using a minimum requirement of three prior matches for both teams. This minimum is imposed in order to prevent the inclusion of early season statistics which are have too little information to be stable.
wvb_prediction_data <- build_prediction_match_data(
wvb_prematch_team,
minimum_prior_matches = 3
)
mvb_prediction_data <- build_prediction_match_data(
mvb_prematch_team,
minimum_prior_matches = 3
)Section 4.3.3: Prediction-Data Coverage
Compares the number of verified matches with the number that contain sufficient pre-match history for inclusion in the prediction model.
prediction_data_coverage <- bind_rows(
tibble(
competition = "Women's Division I",
verified_matches = nrow(wvb_matches),
usable_matches = nrow(wvb_prediction_data),
usable_share =
usable_matches / verified_matches
),
tibble(
competition = "Men's Division I",
verified_matches = nrow(mvb_matches),
usable_matches = nrow(mvb_prediction_data),
usable_share =
usable_matches / verified_matches
)
)
prediction_data_coverage## # A tibble: 2 × 4
## competition verified_matches usable_matches usable_share
## <chr> <int> <int> <dbl>
## 1 Women's Division I 25706 22263 0.866
## 2 Men's Division I 3012 2454 0.815Section 4.3.4: Usable Matches by Source Year
Reports the number of eligible matches in each source year and confirms that both training and test years are represented.
Section 4.4: Fit and Evaluate Logistic Regression
The logistic regression models estimate how the differences in teams’ historical performance are associated with the probability that Team A wins.
Section 4.4.1: Predictor List
Vector which defines the six pre-match performance differences included as predictors in the separate logistic models.
match_predictors <- c(
"diff_prior_win_pct",
"diff_prior_kill_rate",
"diff_prior_error_rate",
"diff_prior_block_rate",
"diff_prior_aces_per_set",
"diff_prior_digs_per_set"
)Section 4.4.2: Binary Evaluation Metrics Helper Function
Helper function which evaluates predicted probabilities using accuracy, log loss and Brier score. Accuracy measures how often the predicted winner was correct. Log loss and Brier score evaluate the quality of the predicted probabilities. Lower log loss and Brier scores indicate better probability estimates.
calculate_binary_metrics <- function(
actual,
probability
) {
probability <- pmin(
pmax(probability, 1e-10),
1 - 1e-10
)
tibble(
accuracy = mean(
(probability >= 0.50) == actual
),
log_loss = -mean(
actual * log(probability) +
(1 - actual) *
log(1 - probability)
),
brier_score = mean(
(probability - actual)^2
)
)
}Section 4.4.3: Logistic Model Fitting Function
Function which divides the data chronologically, fits a logistic regression model on all years before the most recent one, evaluates predictions on the most recent season and compares the model with a simple historical-prevalence baseline. After the held-out evaluation is completed, the function refits the final model using all eligible years so that the full dataset can support future predictions. The baseline assigns every test match the same probabilities based on Team A’s win rate in the training data. Since Team A is alphabetically assigned, the baseline probability should be close to 50%.
fit_match_logistic_model <- function(
match_data,
predictors
) {
modeling_data <- match_data |>
filter(
!is.na(team_a_win),
!is.na(source_year),
if_all(
all_of(predictors),
~ is.finite(.x)
)
)
available_years <- sort(
unique(modeling_data$source_year)
)
if (length(available_years) < 2) {
stop(
"At least two source years are required."
)
}
test_year <- max(
available_years
)
training_data <- modeling_data |>
filter(
source_year < test_year
)
testing_data <- modeling_data |>
filter(
source_year == test_year
)
if (
nrow(training_data) == 0 ||
nrow(testing_data) == 0
) {
stop(
"Training or testing data contain zero matches."
)
}
model_formula <- reformulate(
predictors,
response = "team_a_win"
)
training_model <- glm(
formula = model_formula,
data = training_data,
family = binomial()
)
predicted_probability <- as.numeric(
predict(
training_model,
newdata = testing_data,
type = "response"
)
)
test_predictions <- testing_data |>
mutate(
predicted_probability =
predicted_probability,
predicted_team_a_win = as.integer(
predicted_probability >= 0.50
),
prediction_correct =
predicted_team_a_win ==
team_a_win
)
baseline_probability <- mean(
training_data$team_a_win
)
logistic_metrics <- calculate_binary_metrics(
actual = test_predictions$team_a_win,
probability =
test_predictions$predicted_probability
) |>
mutate(
method = "Logistic regression"
)
baseline_metrics <- calculate_binary_metrics(
actual = test_predictions$team_a_win,
probability = rep(
baseline_probability,
nrow(test_predictions)
)
) |>
mutate(
method =
"Training prevalence baseline"
)
metrics <- bind_rows(
logistic_metrics,
baseline_metrics
) |>
mutate(
test_year = test_year,
training_matches =
nrow(training_data),
test_matches =
nrow(testing_data)
) |>
select(
method,
test_year,
training_matches,
test_matches,
accuracy,
log_loss,
brier_score
)
coefficient_matrix <-
summary(training_model)$coefficients
coefficient_table <- tibble(
term = rownames(
coefficient_matrix
),
estimate =
coefficient_matrix[, "Estimate"],
standard_error =
coefficient_matrix[, "Std. Error"],
p_value =
coefficient_matrix[, "Pr(>|z|)"]
) |>
mutate(
odds_ratio = exp(estimate)
)
# Refit using all eligible years after evaluation.
final_model <- glm(
formula = model_formula,
data = modeling_data,
family = binomial()
)
list(
training_model = training_model,
final_model = final_model,
predictions = test_predictions,
metrics = metrics,
coefficients = coefficient_table,
test_year = test_year,
converged = training_model$converged
)
}Section 4.4.4: Fit Separate Competition Models
Fits separate logistic models to women’s and men’s volleyball.
wvb_match_model <- fit_match_logistic_model(
match_data = wvb_prediction_data,
predictors = match_predictors
)
mvb_match_model <- fit_match_logistic_model(
match_data = mvb_prediction_data,
predictors = match_predictors
)Section 4.5: Model Performance
Evaluation of model performance based on the most recent source year, which wasn’t used to estimate the training model.
Section 4.5.1: Combined Performance Table
Table which compares held-out accuracy, log loss and Brier score of each logistic model with the historical-prevalence baseline. A useful model should aim to have higher accuracy and lower log loss and brier score than the baseline model.
match_model_performance <- bind_rows(
wvb_match_model$metrics |>
mutate(
competition =
"Women's Division I"
),
mvb_match_model$metrics |>
mutate(
competition =
"Men's Division I"
)
) |>
select(
competition,
everything()
)
match_model_performance |>
mutate(
accuracy = round(accuracy, 3),
log_loss = round(log_loss, 3),
brier_score = round(brier_score, 3)
)## # A tibble: 4 × 8
## competition method test_year training_matches test_matches accuracy log_loss
## <chr> <chr> <int> <int> <int> <dbl> <dbl>
## 1 Women's Divi… Logis… 2025 18682 3581 0.727 0.542
## 2 Women's Divi… Train… 2025 18682 3581 0.513 0.693
## 3 Men's Divisi… Logis… 2024 1851 603 0.755 0.499
## 4 Men's Divisi… Train… 2024 1851 603 0.522 0.692
## # ℹ 1 more variable: brier_score <dbl>Section 4.5.2: The Model in Action
wvb_match_model$predictions |>
mutate(
team_a_probability =
scales::percent(
predicted_probability,
accuracy = 0.1
),
team_b_probability =
scales::percent(
1 - predicted_probability,
accuracy = 0.1
),
predicted_winner = if_else(
predicted_probability >= 0.50,
team_a,
team_b
),
actual_winner = if_else(
team_a_win == 1,
team_a,
team_b
)
) |>
select(
match_id,
Season,
match_date,
team_a,
team_b,
team_a_probability,
team_b_probability,
predicted_winner,
actual_winner,
prediction_correct
) |>
arrange(desc(match_date)) |>
slice_head(n = 20)## # A tibble: 20 × 10
## match_id Season match_date team_a team_b team_a_probability
## <chr> <chr> <date> <chr> <chr> <chr>
## 1 wvb_d1__2025__20251129__b… 2025 2025-11-29 Bosto… Virgi… 61.8%
## 2 wvb_d1__2025__20251129__g… 2025 2025-11-29 Georg… Pitts… 22.9%
## 3 wvb_d1__2025__20251129__l… 2025 2025-11-29 Louis… Stanf… 50.3%
## 4 wvb_d1__2025__20251129__m… 2025 2025-11-29 Miami… North… 50.2%
## 5 wvb_d1__2025__20251129__n… 2025 2025-11-29 NC St… SMU 7.3%
## 6 wvb_d1__2025__20251129__n… 2025 2025-11-29 Notre… Wake … 48.8%
## 7 wvb_d1__2025__20251129__a… 2025 2025-11-29 Arizo… West … 79.7%
## 8 wvb_d1__2025__20251129__b… 2025 2025-11-29 BYU Color… 63.5%
## 9 wvb_d1__2025__20251129__h… 2025 2025-11-29 Houst… TCU 27.8%
## 10 wvb_d1__2025__20251129__i… 2025 2025-11-29 Iowa … Kansa… 64.0%
## 11 wvb_d1__2025__20251129__k… 2025 2025-11-29 Kansas UCF 60.3%
## 12 wvb_d1__2025__20251129__t… 2025 2025-11-29 Texas… Utah 30.2%
## 13 wvb_d1__2025__20251129__i… 2025 2025-11-29 India… Purdue 46.9%
## 14 wvb_d1__2025__20251129__i… 2025 2025-11-29 Iowa Penn … 44.5%
## 15 wvb_d1__2025__20251129__n… 2025 2025-11-29 Nebra… Ohio … 98.8%
## 16 wvb_d1__2025__20251128__c… 2025 2025-11-28 Calif… Duke 57.6%
## 17 wvb_d1__2025__20251128__c… 2025 2025-11-28 Clems… Flori… 44.6%
## 18 wvb_d1__2025__20251128__s… 2025 2025-11-28 Syrac… Virgi… 45.8%
## 19 wvb_d1__2025__20251128__a… 2025 2025-11-28 Arizo… Cinci… 41.0%
## 20 wvb_d1__2025__20251128__i… 2025 2025-11-28 Illin… North… 44.9%
## # ℹ 4 more variables: team_b_probability <chr>, predicted_winner <chr>,
## # actual_winner <chr>, prediction_correct <lgl>Section 4.6: Coefficient Interpretation
Logistic regression coefficients describe how each Team A minus Team B predictor difference is associated with Team A’s estimated odds of winning.
Section 4.6.1: Coefficient Table
Combines women’s and men’s model coefficients and reports their estimated direction, odds ratio and statistical significance. Positive coefficients indicate that a larger Team A advantage in the predictor is associated with a higher probability of winning. Negative coefficients indicate that larger value is associated with a lower probability of winning (expected for error-rate).
Raw coefficient magnitudes should NOT be compared directly because the predictors use different measurement scales. Coefficient signs are directly interpretable but comparisons of coefficient magnitude would require predictor standardization.
Section 4.7: Calibration
Evaluates whether the predicted probabilities correspond to the win rates that were actually observed in the held-out test set.
Section 4.7.1: Combine Held-Out Predictions
Combines the held-out women’s and men’s match predictions into one dataset for calibration analysis.
match_test_predictions <- bind_rows(
wvb_match_model$predictions |>
mutate(
competition =
"Women's Division I"
),
mvb_match_model$predictions |>
mutate(
competition =
"Men's Division I"
)
)Section 4.7.2: Calibration Summary
Divides each competitions predictions into ten similarly size probability groups and compares the average predicted probability with the observed Team A win rate in each group. A well-calibrated group with an average predicted probability of 70% should contain an observed Team A win rate close to 70%.
calibration_summary <- match_test_predictions |>
group_by(
competition
) |>
mutate(
probability_group = ntile(
predicted_probability,
10
)
) |>
group_by(
competition,
probability_group
) |>
summarise(
matches = n(),
average_predicted_probability =
mean(predicted_probability),
observed_team_a_win_rate =
mean(team_a_win),
.groups = "drop"
)
calibration_summary## # A tibble: 20 × 5
## competition probability_group matches average_predicted_probability
## <chr> <int> <int> <dbl>
## 1 Men's Division I 1 61 0.0438
## 2 Men's Division I 2 61 0.140
## 3 Men's Division I 3 61 0.224
## 4 Men's Division I 4 60 0.326
## 5 Men's Division I 5 60 0.417
## 6 Men's Division I 6 60 0.501
## 7 Men's Division I 7 60 0.582
## 8 Men's Division I 8 60 0.677
## 9 Men's Division I 9 60 0.773
## 10 Men's Division I 10 60 0.888
## 11 Women's Division I 1 359 0.0839
## 12 Women's Division I 2 358 0.188
## 13 Women's Division I 3 358 0.279
## 14 Women's Division I 4 358 0.369
## 15 Women's Division I 5 358 0.458
## 16 Women's Division I 6 358 0.543
## 17 Women's Division I 7 358 0.621
## 18 Women's Division I 8 358 0.714
## 19 Women's Division I 9 358 0.809
## 20 Women's Division I 10 358 0.915
## # ℹ 1 more variable: observed_team_a_win_rate <dbl>Section 4.7.3: Calibration Plot
Compares predicted and observed win probabilities with the straight dashed diagonal line representing perfect calibration. Points close to the diagonal indicate that the model’s probabilities are trustworthy. Points above the line indicate an underestimation of Team A’s win probability while points below indicate overestimation.
ggplot(
calibration_summary,
aes(
x = average_predicted_probability,
y = observed_team_a_win_rate
)
) +
geom_point() +
geom_line() +
geom_abline(
slope = 1,
intercept = 0,
linetype = "dashed"
) +
facet_wrap(
vars(competition)
) +
coord_equal(
xlim = c(0, 1),
ylim = c(0, 1)
) +
labs(
title =
"Match Outcome Model Calibration",
subtitle =
"Observed and predicted Team A win rates in probability groups",
x = "Average predicted probability",
y = "Observed win rate"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 20, face = "bold"),
plot.subtitle = element_text(size = 14),
strip.text = element_text(size = 12),
axis.title = element_text(size = 16),
axis.text.x = element_text(size = 11),
axis.text.y = element_text(size = 12)
)
The women’s model closely follows the diagonal across most probability groups , indicating strong calibration. The men’s model is also reasonably calibrated but has more fluctuation. This is consistent with the smaller test sample.
Section 4.8: Limitations
The logistic models used are supposed to strike a balance between complexity and interpretability. While it does improve on the baseline model, it also doesn’t account for many other factors which can influence match outcomes. For example, opponent strength, venue advantages, injuries, lineup changes, travel, etc. Additionally, the predictors are season-to-date averages which can respond very slowly to sudden increases or decreases in performance during a season.
Matches were excluded until both Team A and Team B have completed at least three prior matches. This means the model doesn’t generate predictions for the earliest part of each season. Also, the held-out evaluation uses only one source year so performance may vary in future seasons.
The model should be interpreted as fairly intuitive benchmark for estimating match-win probabilities rather than a comprehensive match forecasting system.
Section 5: Conference Parity and Championship Analysis
Evaluates similarity in offensive and defensive performance within conferences and examines the statistical profiles of recent national champions.
Section 5.1: Construct Conference-Match Records
Section 5.1.1: Conference Record Construction Function
Joins team-season conference affiliations to be verified match records, retains same-conference matches and calculates team’s various metrics. Since the original data contains one row for each team, every conference match contributes one record to each participating team.
build_conference_team_records <- function(
team_rates,
team_season
) {
conference_lookup <- team_season |>
transmute(
source_year = as.integer(source_year),
team_clean = clean_team_name(Team),
Conference = stringr::str_squish(
as.character(Conference)
)
) |>
filter(
!is.na(Conference),
Conference != ""
) |>
distinct(
source_year,
team_clean,
.keep_all = TRUE
)
team_conference_lookup <- conference_lookup |>
rename(
team_conference = Conference
)
opponent_conference_lookup <- conference_lookup |>
transmute(
source_year,
opponent_clean = team_clean,
opponent_conference = Conference
)
team_rates |>
mutate(
source_year = as.integer(source_year),
kills_num = as.numeric(Kills),
errors_num = as.numeric(Errors),
attacks_num = as.numeric(`Total Attacks`)
) |>
left_join(
team_conference_lookup,
by = c(
"source_year",
"team_clean"
)
) |>
left_join(
opponent_conference_lookup,
by = c(
"source_year",
"opponent_clean"
)
) |>
filter(
!is.na(team_conference),
!is.na(opponent_conference),
team_conference == opponent_conference
) |>
group_by(
source_year,
Season,
Conference = team_conference,
team_clean
) |>
summarise(
conference_matches =
n_distinct(match_id),
conference_wins =
sum(win == 1L, na.rm = TRUE),
conference_losses =
sum(win == 0L, na.rm = TRUE),
conference_sets_won =
sum(sets_won, na.rm = TRUE),
conference_sets_lost =
sum(sets_lost, na.rm = TRUE),
total_kills =
sum(kills_num, na.rm = TRUE),
total_errors =
sum(errors_num, na.rm = TRUE),
total_attacks =
sum(attacks_num, na.rm = TRUE),
total_block_points =
sum(
Block_Solos_clean +
0.5 * Block_Assists_clean,
na.rm = TRUE
),
total_aces =
sum(Aces_clean, na.rm = TRUE),
total_digs =
sum(Digs_clean, na.rm = TRUE),
total_sets =
sum(sets_num, na.rm = TRUE),
.groups = "drop"
) |>
mutate(
conference_win_pct =
safe_divide(
conference_wins,
conference_wins +
conference_losses
),
conference_set_pct =
safe_divide(
conference_sets_won,
conference_sets_won +
conference_sets_lost
),
kill_rate_pct =
100 * safe_divide(
total_kills,
total_attacks
),
error_rate_pct =
100 * safe_divide(
total_errors,
total_attacks
),
block_rate_pct =
100 * safe_divide(
total_block_points,
total_attacks
),
aces_per_set =
safe_divide(
total_aces,
total_sets
),
digs_per_set =
safe_divide(
total_digs,
total_sets
)
)
}Section 5.1.2: Apply Function to Both Women’s and Men’s
Constructs women’s and men’s conference performance records and combines them into one analysis table.
wvb_conference_team_records <-
build_conference_team_records(
team_rates = wvb_team_rates,
team_season = wvb_team_season
) |>
mutate(
competition = "Women's Division I"
)
mvb_conference_team_records <-
build_conference_team_records(
team_rates = mvb_team_rates,
team_season = mvb_team_season
) |>
mutate(
competition = "Men's Division I"
)
conference_team_records <- bind_rows(
wvb_conference_team_records,
mvb_conference_team_records
)Section 5.1.3: Conference-Match Coverage
Reports the number of teams, team-match records and estimated unique conference matches available for each conference and season. Unique conference matches are calculated by dividing team-match record count by two since each verified appears once from each team’s perspective.
conference_match_coverage <-
conference_team_records |>
group_by(
competition,
source_year,
Conference
) |>
summarise(
teams = n_distinct(team_clean),
team_match_records =
sum(conference_matches),
unique_conference_matches =
team_match_records / 2,
.groups = "drop"
)
conference_match_coverage## # A tibble: 225 × 6
## competition source_year Conference teams team_match_records
## <chr> <int> <chr> <int> <int>
## 1 Men's Division I 2020 Big West 6 78
## 2 Men's Division I 2020 Conference Carolinas 8 106
## 3 Men's Division I 2020 DI Independent 8 48
## 4 Men's Division I 2020 EIVA 6 96
## 5 Men's Division I 2020 MIVA 8 112
## 6 Men's Division I 2020 MPSF 7 122
## 7 Men's Division I 2021 Big West 6 60
## 8 Men's Division I 2021 Conference Carolinas 8 112
## 9 Men's Division I 2021 DI Independent 10 82
## 10 Men's Division I 2021 EIVA 9 144
## # ℹ 215 more rows
## # ℹ 1 more variable: unique_conference_matches <dbl>Section 5.2: Measure Conference Parity
Section 5.2.1: Calculate Metric Parity Function
Function which finds standard deviation of relevant metrics for parity analysis.
build_conference_metric_parity <- function(
conference_team_records,
minimum_teams = 4,
minimum_conference_matches = 3
) {
metric_variables <- c(
"kill_rate_pct",
"error_rate_pct",
"block_rate_pct",
"aces_per_set",
"digs_per_set"
)
conference_team_records |>
filter(
conference_matches >=
minimum_conference_matches
) |>
pivot_longer(
cols = all_of(metric_variables),
names_to = "metric",
values_to = "team_metric_value"
) |>
filter(
is.finite(team_metric_value)
) |>
group_by(
competition,
source_year,
Season,
Conference,
metric
) |>
summarise(
teams =
n_distinct(team_clean),
metric_mean =
mean(
team_metric_value,
na.rm = TRUE
),
metric_median =
median(
team_metric_value,
na.rm = TRUE
),
metric_sd =
sd(
team_metric_value,
na.rm = TRUE
),
metric_iqr =
IQR(
team_metric_value,
na.rm = TRUE
),
metric_range =
max(
team_metric_value,
na.rm = TRUE
) -
min(
team_metric_value,
na.rm = TRUE
),
.groups = "drop"
) |>
filter(
teams >= minimum_teams,
is.finite(metric_sd)
) |>
group_by(
competition,
source_year,
metric
) |>
mutate(
metric_parity_rank =
min_rank(metric_sd),
metric_parity_score = if (dplyr::n() > 1L) {
100 * (1 - dplyr::percent_rank(metric_sd))
} else {
rep(100, dplyr::n())
}
) |>
ungroup() |>
mutate(
metric_label = recode(
metric,
kill_rate_pct = "Kill rate",
error_rate_pct = "Attack error rate",
block_rate_pct = "Block rate",
aces_per_set = "Aces per set",
digs_per_set = "Digs per set"
)
)
}Section 5.2.2: Average Metric Parity Across Seasons
Applies the conference parity rules and create a summary record for each eligible conference-season.
conference_metric_parity <-
build_conference_metric_parity(
conference_team_records,
minimum_teams = 4,
minimum_conference_matches = 3
)Summary of the average metric standard deviations achieved by various conferences over the course of multiple seasons.
conference_metric_parity_summary <-
conference_metric_parity |>
group_by(
competition,
Conference,
metric,
metric_label
) |>
summarise(
mean_sd = mean(
metric_sd,
na.rm = TRUE
),
seasons = n_distinct(source_year),
.groups = "drop"
)
conference_metric_parity_summary## # A tibble: 210 × 6
## competition Conference metric metric_label mean_sd seasons
## <chr> <chr> <chr> <chr> <dbl> <int>
## 1 Men's Division I Big West aces_per_… Aces per set 0.265 5
## 2 Men's Division I Big West block_rat… Block rate 1.95 5
## 3 Men's Division I Big West digs_per_… Digs per set 0.723 5
## 4 Men's Division I Big West error_rat… Attack erro… 2.41 5
## 5 Men's Division I Big West kill_rate… Kill rate 2.95 5
## 6 Men's Division I Conference Carolinas aces_per_… Aces per set 0.282 5
## 7 Men's Division I Conference Carolinas block_rat… Block rate 1.47 5
## 8 Men's Division I Conference Carolinas digs_per_… Digs per set 0.938 5
## 9 Men's Division I Conference Carolinas error_rat… Attack erro… 2.76 5
## 10 Men's Division I Conference Carolinas kill_rate… Kill rate 3.63 5
## # ℹ 200 more rowsSection 5.2.3: Creating Parity Percentiles
Conference performance parity is measured by within-conference standard deviation of team performance metrics. Lower values indicate that teams produced more similar statistical results. Here, parity is converted into a percentile rank showing how the parity rank of the conference in that metric.
First, we’ll do it for women’s volleyball.
conference_parity_scores_women <- conference_metric_parity_summary |>
filter(
competition == "Women's Division I"
) |>
group_by(
metric_label
) |>
mutate(
parity_score = if (dplyr::n() > 1L) {
100 * (1 - dplyr::percent_rank(mean_sd))
} else {
100
}
) |>
ungroup()
selected_conferences_women <- c("Big Ten", "SEC", "Big 12", "ACC")Then for men’s volleyball.
conference_parity_scores_men <- conference_metric_parity_summary |>
filter(
competition == "Men's Division I"
) |>
group_by(
metric_label
) |>
mutate(
parity_score = if (dplyr::n() > 1L) {
100 * (1 - dplyr::percent_rank(mean_sd))
} else {
100
}
) |>
ungroup()
selected_conferences_men <- c("Big West", "MPSF", "MIVA")Section 5.2.4: Create Radar Data
Pulling together the data to be used in the radar plots.
First for women.
radar_data_women <- conference_parity_scores_women |>
filter(
competition == "Women's Division I",
Conference %in% selected_conferences_women
) |>
select(
Conference,
metric_label,
parity_score
) |>
pivot_wider(
names_from = metric_label,
values_from = parity_score
)Then for men.
radar_data_men <- conference_parity_scores_men |>
filter(
competition == "Men's Division I",
Conference %in% selected_conferences_men
) |>
select(
Conference,
metric_label,
parity_score
) |>
pivot_wider(
names_from = metric_label,
values_from = parity_score
)Section 5.2.5: Radar Plot Function
Function which plots the percentiles for each metric on the radar plot alongside the average percentile for all conferences (which is clearly going to be the 50th percentile).
plot_conference_radar_with_average <- function(
data,
selected_competition,
conference_name
) {
metric_levels <- c(
"Kill rate",
"Attack error rate",
"Block rate",
"Aces per set",
"Digs per set"
)
comp_data <- data |>
filter(
competition == selected_competition
) |>
mutate(
metric_label = factor(
metric_label,
levels = metric_levels
)
) |>
group_by(metric_label) |>
mutate(
parity_score = if (dplyr::n() > 1L) {
100 * (1 - dplyr::percent_rank(mean_sd))
} else {
100
}
) |>
ungroup()
average_profile <- comp_data |>
group_by(metric_label) |>
summarise(
parity_score = mean(parity_score, na.rm = TRUE),
.groups = "drop"
) |>
arrange(metric_label) |>
pull(parity_score)
selected_profile <- comp_data |>
filter(
Conference == conference_name
) |>
arrange(metric_label) |>
pull(parity_score)
radar_df <- rbind(
c(100, 100, 100, 100, 100),
c(0, 0, 0, 0, 0),
average_profile,
selected_profile
) |>
as.data.frame()
colnames(radar_df) <- metric_levels
fmsb::radarchart(
radar_df,
axistype = 1,
pcol = c("grey60", "steelblue"),
pfcol = c(
scales::alpha("grey60", 0.15),
scales::alpha("steelblue", 0.25)
),
plwd = c(2, 2),
plty = c(2, 1),
cglcol = "grey85",
cglty = 1,
axislabcol = "grey40",
caxislabels = c("0", "25", "50", "75", "100"),
vlcex = 0.9,
title = conference_name
)
}Section 5.2.6: Applying the Plotting Function
Applying the plotting function with some visual modifications.
First for women.
par(
mfrow = c(2, 2),
mar = c(1, 1, 3, 1),
oma = c(2, 0, 3, 0)
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Women's Division I",
conference_name = "Big Ten"
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Women's Division I",
conference_name = "SEC"
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Women's Division I",
conference_name = "Big 12"
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Women's Division I",
conference_name = "ACC"
)
mtext(
"Women's Division I Conference Parity Profiles",
side = 3,
outer = TRUE,
line = 1,
cex = 1.8,
font = 2
)
par(xpd = NA)
legend(
x = -1.75,
y = -1.25,
legend = c(
"Conference Average",
"Selected Conference"
),
col = c(
"grey60",
"steelblue"
),
lwd = 2,
horiz = TRUE,
bty = "n",
xjust = 0.5
)
Then for men.
par(
mfrow = c(2, 2),
mar = c(1, 1, 3, 1),
oma = c(2, 0, 3, 0)
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Men's Division I",
conference_name = "Big West"
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Men's Division I",
conference_name = "MPSF"
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Men's Division I",
conference_name = "MIVA"
)
plot_conference_radar_with_average(
conference_metric_parity_summary,
selected_competition = "Men's Division I",
conference_name = "EIVA"
)
mtext(
"Men's Division I Conference Parity Profiles",
side = 3,
outer = TRUE,
line = 1,
cex = 1.8,
font = 2
)
par(xpd = NA)
legend(
x = -1.75,
y = -1.25,
legend = c(
"Conference Average",
"Selected Conference"
),
col = c(
"grey60",
"steelblue"
),
lwd = 2,
horiz = TRUE,
bty = "n",
xjust = 0.5
)
The blue shape represents the selected conference while the gray represents the average for all conferences. Father from the center means greater parity (team within that conference were more similar to one another in that metric). High parity doesn’t equate to high performance. High kill-rate parity means that all of the teams in that conference either had high or low kill rates.
Section 5.3: National Champion Performance Profiles
Creates team-season performance profiles and converts each metric into a within-season percentile so national champions can be compared across seasons and competitions.
Section 5.3.1: Official Champion Look-up Table
Records official nation champion for both women’s and men’s volleyball and for each year considered. Also, team names were standardized.
official_champions <- tribble(
~competition, ~source_year, ~champion_team,
"Women's Division I", 2021, "Wisconsin",
"Women's Division I", 2022, "Texas",
"Women's Division I", 2023, "Texas",
"Women's Division I", 2024, "Penn St.",
"Women's Division I", 2025, "Texas A&M",
"Men's Division I", 2021, "Hawaii",
"Men's Division I", 2022, "Hawaii",
"Men's Division I", 2023, "UCLA",
"Men's Division I", 2024, "UCLA"
) |>
mutate(
source_year = as.integer(source_year),
team_clean = clean_team_name(champion_team)
)Section 5.3.2: Safe-Mean Helper Function
Calculates a mean using only finite observations and returns a missing value when no valid observations are available.
safe_mean <- function(x) {
valid_x <- x[
is.finite(x)
]
if (length(valid_x) == 0) {
return(NA_real_)
}
mean(valid_x)
}Section 5.3.3: Team-Season Metric Function
Aggregates verified team-match data into season-level records containing win percentage, set percentage and attacking, error, blocking, serving and digging performance.
build_team_season_championship_metrics <- function(
team_rates,
competition_label
) {
team_rates |>
mutate(
source_year =
as.integer(source_year)
) |>
group_by(
source_year,
Season,
team_clean
) |>
summarise(
matches =
n_distinct(match_id),
win_pct =
safe_mean(win),
set_pct =
safe_divide(
sum(
sets_won,
na.rm = TRUE
),
sum(
sets_won,
na.rm = TRUE
) +
sum(
sets_lost,
na.rm = TRUE
)
),
kill_rate_pct =
100 * safe_divide(
sum(Kills, na.rm = TRUE),
sum(`Total Attacks`, na.rm = TRUE)
),
error_rate_pct =
100 * safe_divide(
sum(Errors, na.rm = TRUE),
sum(`Total Attacks`, na.rm = TRUE)
),
block_rate_pct =
100 * safe_divide(
sum(Block_Solos_clean + 0.5 * Block_Assists_clean, na.rm = TRUE),
sum(`Total Attacks`, na.rm = TRUE)
),
aces_per_set =
safe_divide(
sum(Aces_clean, na.rm = TRUE),
sum(sets_num, na.rm = TRUE)
),
digs_per_set =
safe_divide(
sum(Digs_clean, na.rm = TRUE),
sum(sets_num, na.rm = TRUE)
),
.groups = "drop"
) |>
mutate(
competition =
competition_label
) |>
select(
competition,
everything()
)
}Section 5.3.4: Apply the Function to Both Women’s and Men’s
Creates and combines the women’s and men’s team-season performance profiles.
wvb_team_season_metrics <-
build_team_season_championship_metrics(
team_rates = wvb_team_rates,
competition_label =
"Women's Division I"
)
mvb_team_season_metrics <-
build_team_season_championship_metrics(
team_rates = mvb_team_rates,
competition_label =
"Men's Division I"
)
team_season_championship_metrics <-
bind_rows(
wvb_team_season_metrics,
mvb_team_season_metrics
)Section 5.3.5: Team-Season Coverage
Reports how many team-season profiles are available.
team_season_championship_metrics |>
count(
competition,
source_year,
name = "teams"
)## # A tibble: 11 × 3
## competition source_year teams
## <chr> <int> <int>
## 1 Men's Division I 2020 43
## 2 Men's Division I 2021 55
## 3 Men's Division I 2022 57
## 4 Men's Division I 2023 59
## 5 Men's Division I 2024 66
## 6 Women's Division I 2020 308
## 7 Women's Division I 2021 340
## 8 Women's Division I 2022 344
## 9 Women's Division I 2023 344
## 10 Women's Division I 2024 346
## 11 Women's Division I 2025 314Section 5.3.6: Convert Team Metrics to Percentiles
Converts each performance metric int a percentile relative to teams in the same gender and source year and create composite offensive and defensive scores. Higher percentiles indicate stronger performance. Error rate is reversed to create an error-control percentile so that higher values represent better results. The offensive composite score averages kill rate, error control and aces per set. The defensive composite score averages block rate and digs per set.
team_season_percentiles <-
team_season_championship_metrics |>
group_by(
competition,
source_year
) |>
mutate(
win_pct_percentile =
100 * percent_rank(
win_pct
),
set_pct_percentile =
100 * percent_rank(
set_pct
),
kill_rate_percentile =
100 * percent_rank(
kill_rate_pct
),
error_control_percentile =
100 * (
1 - percent_rank(
error_rate_pct
)
),
block_rate_percentile =
100 * percent_rank(
block_rate_pct
),
aces_per_set_percentile =
100 * percent_rank(
aces_per_set
),
digs_per_set_percentile =
100 * percent_rank(
digs_per_set
)
) |>
ungroup() |>
mutate(
offense_score = rowMeans(
cbind(
kill_rate_percentile,
error_control_percentile,
aces_per_set_percentile
),
na.rm = TRUE
),
defense_score = rowMeans(
cbind(
block_rate_percentile,
digs_per_set_percentile
),
na.rm = TRUE
)
)Section 5.3.7: Join Champion Profiles
Joins the official champion list to the corresponding team-season percentile profiles.
champion_team_profiles <-
official_champions |>
left_join(
team_season_percentiles,
by = c(
"competition",
"source_year",
"team_clean"
)
)Section 5.4: National Champion Visuals
Section 5.4.1: Reshape Champion Metrics
Reshapes the champion percentile measures into long format and assigns labels to each metric.
champion_metric_profiles <-
champion_team_profiles |>
pivot_longer(
cols = c(
win_pct_percentile,
set_pct_percentile,
kill_rate_percentile,
error_control_percentile,
block_rate_percentile,
aces_per_set_percentile,
digs_per_set_percentile
),
names_to = "metric",
values_to = "percentile"
) |>
mutate(
metric_label = case_when(
metric ==
"win_pct_percentile" ~
"Win percentage",
metric ==
"set_pct_percentile" ~
"Set percentage",
metric ==
"kill_rate_percentile" ~
"Kill rate",
metric ==
"error_control_percentile" ~
"Error control",
metric ==
"block_rate_percentile" ~
"Block rate",
metric ==
"aces_per_set_percentile" ~
"Aces per set",
metric ==
"digs_per_set_percentile" ~
"Digs per set",
TRUE ~ metric
),
metric_label = factor(
metric_label,
levels = c(
"Win percentage",
"Set percentage",
"Kill rate",
"Error control",
"Aces per set",
"Block rate",
"Digs per set"
)
)
)Section 5.4.2: Champion Heatmap
Shows how each national champion ranked relative to other teams in the same gender and source year across seven performance measures. Values near 100 indicate that the champion ranked among the strongest teams in that metric while lower values identify areas where the champion was less dominant.
ggplot(
champion_metric_profiles,
aes(
x = factor(source_year),
y = metric_label,
fill = percentile
)
) +
geom_tile(
linewidth = 0.5
) +
geom_text(
aes(
label = round(
percentile
)
),
size = 3
) +
facet_wrap(
vars(competition),
scales = "free_x"
) +
scale_fill_viridis_c(
limits = c(0, 100),
labels = scales::percent_format(
scale = 1,
accuracy = 1
)
) +
labs(
title =
"Performance Profiles of National Champions",
subtitle =
"Percentile rank relative to teams in the same competition and source year",
x = "Championship year",
y = NULL,
fill = "National\npercentile"
) +
theme_minimal() +
theme(
panel.grid = element_blank(),
axis.text.x = element_text(
angle = 45,
hjust = 1
)
)
Section 5.4.3: Champion Benchmark Summary
Summarizes the distribution of champion percentiles for each metric using the median, IQR and shares reaching the top quartile/decile.
champion_metric_benchmarks <-
champion_metric_profiles |>
filter(
is.finite(percentile)
) |>
group_by(
competition,
metric_label
) |>
summarise(
champions =
n_distinct(source_year),
median_percentile =
median(
percentile,
na.rm = TRUE
),
lower_quartile =
quantile(
percentile,
0.25,
na.rm = TRUE
),
upper_quartile =
quantile(
percentile,
0.75,
na.rm = TRUE
),
top_quartile_share =
mean(
percentile >= 75,
na.rm = TRUE
),
top_decile_share =
mean(
percentile >= 90,
na.rm = TRUE
),
.groups = "drop"
)Section 5.4.4: Champion Benchmark Classification
Classifies each metric as a consistent champion hallmark, common champion strength or a variable characteristic based on its median percentile and top-quartile frequency. The labels are descriptive and aren’t statistically estimated championship requirements.
champion_contender_benchmarks <-
champion_metric_benchmarks |>
mutate(
profile_interpretation = case_when(
median_percentile >= 80 &
top_quartile_share >= 0.75 ~
"Consistent champion hallmark",
median_percentile >= 70 &
top_quartile_share >= 0.50 ~
"Common champion strength",
TRUE ~
"Variable among champions"
)
) |>
arrange(
competition,
desc(median_percentile)
)
champion_contender_benchmarks## # A tibble: 14 × 9
## competition metric_label champions median_percentile lower_quartile
## <chr> <fct> <int> <dbl> <dbl>
## 1 Men's Division I Kill rate 4 97.4 95.0
## 2 Men's Division I Block rate 4 96.1 93.2
## 3 Men's Division I Set percentage 4 95.5 93.1
## 4 Men's Division I Aces per set 4 95.1 93.6
## 5 Men's Division I Win percentage 4 94.6 92.5
## 6 Men's Division I Error control 4 93.7 91.4
## 7 Men's Division I Digs per set 4 30.7 17.2
## 8 Women's Division I Set percentage 5 98.1 97.7
## 9 Women's Division I Win percentage 5 97.9 97.4
## 10 Women's Division I Kill rate 5 97.7 91.8
## 11 Women's Division I Block rate 5 91.7 88.7
## 12 Women's Division I Error control 5 88.5 70.4
## 13 Women's Division I Aces per set 5 63.1 44.3
## 14 Women's Division I Digs per set 5 41.1 37.4
## # ℹ 4 more variables: upper_quartile <dbl>, top_quartile_share <dbl>,
## # top_decile_share <dbl>, profile_interpretation <chr>Section 5.4.5: Champion Identifiers
Creates a compact look-up table identifying the national champions.
champion_keys <- official_champions |>
select(
competition,
source_year,
team_clean,
champion_team
) |>
mutate(
national_champion = TRUE
)Section 5.4.6: Contender Profiles
Adds champion indicators and labels to the complete set of team-season offensive and defensive profiles.
team_contender_profiles <-
team_season_percentiles |>
left_join(
champion_keys,
by = c(
"competition",
"source_year",
"team_clean"
)
) |>
mutate(
national_champion =
coalesce(
national_champion,
FALSE
),
champion_label = if_else(
national_champion,
paste0(
champion_team,
" (",
source_year,
")"
),
NA_character_
)
)
team_contender_profiles## # A tibble: 2,276 × 24
## competition source_year Season team_clean matches win_pct set_pct
## <chr> <int> <chr> <chr> <int> <dbl> <dbl>
## 1 Women's Division I 2020 2020 A&M-Corpus Chr… 16 0.688 0.667
## 2 Women's Division I 2020 2020 Abilene Christ… 14 0.5 0.491
## 3 Women's Division I 2020 2020 Air Force 13 0.538 0.585
## 4 Women's Division I 2020 2020 Akron 22 0.273 0.367
## 5 Women's Division I 2020 2020 Alabama 22 0.318 0.341
## 6 Women's Division I 2020 2020 Alcorn 9 0.333 0.367
## 7 Women's Division I 2020 2020 American 7 0.571 0.556
## 8 Women's Division I 2020 2020 App State 9 0.111 0.194
## 9 Women's Division I 2020 2020 Arizona 21 0.476 0.507
## 10 Women's Division I 2020 2020 Arizona St. 20 0.3 0.392
## # ℹ 2,266 more rows
## # ℹ 17 more variables: kill_rate_pct <dbl>, error_rate_pct <dbl>,
## # block_rate_pct <dbl>, aces_per_set <dbl>, digs_per_set <dbl>,
## # win_pct_percentile <dbl>, set_pct_percentile <dbl>,
## # kill_rate_percentile <dbl>, error_control_percentile <dbl>,
## # block_rate_percentile <dbl>, aces_per_set_percentile <dbl>,
## # digs_per_set_percentile <dbl>, offense_score <dbl>, defense_score <dbl>, …Section 5.4.7: Offense-Defense Plot
Compares every team-season’s offensive and defensive percentile composites and highlights the eventual national champions. The dashed line mark the median offensive and defensive performance. Teams in the upper-right quadrant performed above the national median in both composites.
Most observed champions combined strong offensive performance with at least above-average defensive performance but it varies across teams and seasons. Several non-champions also occupy the upper-right quadrant suggesting that strong regular season performance is sometimes associated with post-season success but doesn’t guarantee a national title.
ggplot(
team_contender_profiles,
aes(
x = offense_score,
y = defense_score
)
) +
geom_point(alpha = 0.15) +
geom_vline(
xintercept = 50,
linetype = "dashed"
) +
geom_hline(
yintercept = 50,
linetype = "dashed"
) +
geom_point(
data = team_contender_profiles |>
filter(national_champion),
size = 3
) +
ggrepel::geom_text_repel(
data = team_contender_profiles |>
filter(national_champion),
aes(label = champion_label),
size = 7,
box.padding = 0.5,
point.padding = 0.3,
min.segment.length = 0,
segment.alpha = 1,
segment.size = 2,
seed = 123
) +
facet_wrap(vars(competition)) +
coord_cartesian(
xlim = c(0, 105),
ylim = c(0, 105)
) +
labs(
title = "Offensive and Defensive Profiles of National Champions",
subtitle = "Dashed lines mark the 50th percentile in each composite",
x = "Offensive Production (Kills, Error Control and Aces)",
y = "Defensive Production (Blocks and Digs)"
) +
theme_minimal() +
theme(
plot.title = element_text(size = 20, face = "bold"),
plot.subtitle = element_text(size = 14),
strip.text = element_text(size = 15),
axis.title = element_text(size = 20),
axis.text.x = element_text(size = 14),
axis.text.y = element_text(size = 14)
)
Section 5.5: Limitations
Conference parity measures similarity measures similarity among teams rather than overall conference strength. Conference membership and scheduling structures can also change across seasons and some conferences are represented in fewer seasons than others.
Additional Code which was used in the application section
mruzik_pei_path <- wvb_player_pei |>
filter(
Player == "Jess Mruzik",
Season %in% c("2022", "2023", "2024")
) |>
arrange(
as.integer(Season)
) |>
mutate(
stage = paste(
Team,
Season,
sep = "\n"
),
stage = factor(
stage,
levels = stage
)
)
mruzik_pei_path## # A tibble: 3 × 64
## Team Conference Number Player Yr Pos Ht GP GS S MS Kills
## <chr> <chr> <dbl> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Mich… Big Ten 5 Jess … Jr OH 6-1 30 30 107 NA 394
## 2 Penn… Big Ten 9 Jess … Sr OH 6-1 32 32 118 NA 519
## 3 Penn… Big Ten 9 Jess … Sr OH 6-1 37 37 130 NA 565
## # ℹ 52 more variables: Errors <dbl>, `Total Attacks` <dbl>, `Hit Pct` <dbl>,
## # Assists <dbl>, Aces <dbl>, SErr <dbl>, Digs <dbl>, RErr <dbl>,
## # `Block Solos` <dbl>, `Block Assists` <dbl>, BErr <dbl>, TB <dbl>,
## # PTS <dbl>, BHE <dbl>, `Trpl Dbl` <dbl>, source_year <int>, SrvAtt <dbl>,
## # RetAtt <dbl>, `Dbl Dbl` <dbl>, Hometown <chr>, `High School` <chr>,
## # TeamID <dbl>, Season <chr>, original_player_season <chr>, team_clean <chr>,
## # pos_group <chr>, player_sets <dbl>, kills_num <dbl>, …ggplot(
mruzik_pei_path,
aes(
x = stage,
y = PEI,
group = 1
)
) +
geom_hline(
yintercept = 50,
linetype = "dashed"
) +
geom_line(
linewidth = 1
) +
geom_point(
size = 5
) +
geom_text(
aes(
label = paste0(
"PEI: ",
PEI,
"\n",
position_percentile,
" percentile"
)
),
vjust = -1,
size = 4
) +
labs(
title = "Jess Mruzik: PEI Across Transfer to Penn State",
subtitle =
"Penn State won the 2024 national championship",
x = NULL,
y = "Player Efficiency Index"
) +
theme_minimal()
wvb_match_model$predictions |>
filter(
team_a == "Texas A&M" |
team_b == "Texas A&M"
) |>
arrange(
desc(match_date)
) |>
select(
match_date,
team_a,
team_b,
predicted_probability,
team_a_win
) |>
slice_head(n = 10)## # A tibble: 10 × 5
## match_date team_a team_b predicted_probability team_a_win
## <date> <chr> <chr> <dbl> <int>
## 1 2025-11-16 Georgia Texas A&M 0.253 0
## 2 2025-11-12 Florida Texas A&M 0.257 0
## 3 2025-11-07 Auburn Texas A&M 0.336 0
## 4 2025-11-02 Tennessee Texas A&M 0.558 0
## 5 2025-10-31 Texas Texas A&M 0.632 0
## 6 2025-10-26 Ole Miss Texas A&M 0.132 0
## 7 2025-10-24 LSU Texas A&M 0.161 0
## 8 2025-10-19 Arkansas Texas A&M 0.109 0
## 9 2025-10-17 Oklahoma Texas A&M 0.293 0
## 10 2025-10-12 Mississippi St. Texas A&M 0.520 0wvb_match_model$predictions |>
filter(
team_a == "Texas A&M" & team_b == "Kentucky" |
team_a == "Kentucky" & team_b == "Texas A&M"
)## # A tibble: 1 × 30
## match_id source_year Season match_date team_a team_b team_a_win
## <chr> <int> <chr> <date> <chr> <chr> <int>
## 1 wvb_d1__2025__20251008… 2025 2025 2025-10-08 Kentu… Texas… 1
## # ℹ 23 more variables: team_a_prior_matches <int>, team_a_prior_win_pct <dbl>,
## # team_a_prior_kill_rate <dbl>, team_a_prior_error_rate <dbl>,
## # team_a_prior_block_rate <dbl>, team_a_prior_aces_per_set <dbl>,
## # team_a_prior_digs_per_set <dbl>, team_b_prior_matches <int>,
## # team_b_prior_win_pct <dbl>, team_b_prior_kill_rate <dbl>,
## # team_b_prior_error_rate <dbl>, team_b_prior_block_rate <dbl>,
## # team_b_prior_aces_per_set <dbl>, team_b_prior_digs_per_set <dbl>, …selected_champion_profile <-
champion_metric_profiles |>
filter(
competition ==
"Women's Division I",
source_year == 2023
)ggplot(
selected_champion_profile,
aes(
x = percentile,
y = reorder(
metric_label,
percentile
)
)
) +
geom_col() +
geom_vline(
xintercept = 75,
linetype = "dashed"
) +
geom_text(
aes(
label = paste0(
round(percentile),
"th"
)
),
hjust = 1.15,
color = "white",
size = 4
) +
scale_x_continuous(
limits = c(0, 100),
breaks = seq(
0,
100,
by = 25
)
) +
labs(
title = paste(
unique(
selected_champion_profile$
champion_team
),
"2023 Championship Profile"
),
subtitle =
"Percentile relative to all Division I teams in the same source year",
x = "National percentile",
y = NULL
) +
theme_minimal()
selected_teams <- c(
"Texas",
"Nebraska"
)championship_radar_data <-
team_season_percentiles |>
filter(
competition == "Women's Division I",
source_year == 2023,
team_clean %in% clean_team_name(selected_teams)
) |>
transmute(
Team = team_clean,
`Win %` = win_pct_percentile,
`Set %` = set_pct_percentile,
`Kill Rate` = kill_rate_percentile,
`Error Control` = error_control_percentile,
`Block Rate` = block_rate_percentile,
`Aces / Set` = aces_per_set_percentile,
`Digs / Set` = digs_per_set_percentile
)plot_team_profile <- function(
radar_data,
team_name
) {
one_team <-
radar_data |>
filter(
Team == team_name
) |>
select(-Team)
radar_df <-
rbind(
rep(100, ncol(one_team)),
rep(0, ncol(one_team)),
one_team
)
colnames(radar_df) <- names(one_team)
radarchart(
radar_df,
axistype = 1,
pcol = "steelblue",
pfcol = scales::alpha(
"steelblue",
0.25
),
plwd = 2,
plty = 1,
cglcol = "grey80",
cglty = 1,
axislabcol = "grey40",
caxislabels =
c(
"0",
"25",
"50",
"75",
"100"
),
vlcex = 0.9,
title = team_name
)
}par(
mfrow = c(1, 2),
mar = c(2, 2, 3, 2),
oma = c(0, 0, 3, 0)
)
plot_team_profile(
championship_radar_data,
"Texas"
)
plot_team_profile(
championship_radar_data,
"Nebraska"
)
mtext(
"2023 National Champion vs Elite Non-Champion",
side = 3,
outer = TRUE,
cex = 1.6,
font = 2
)