1 Introduction

1.1 Foreword

In the world of professional powerlifting, the ability to compare athletes of different body weights, ages, and genders is crucial for assessing their strength. To facilitate such comparisons, powerlifting relies on a variety of formulas that standardize performance metrics, producing a comparable score that transcends differences in individual characteristics. However, these formulas can vary significantly in how they adjust for factors like body weight, age, and gender, leading to potentially differing perceptions of an athlete’s relative strength. This discrepancy can influence rankings and outcomes, raising questions about the accuracy and fairness of these formulas.

The objective of this study is to explore and compare four of the most commonly used powerlifting formulas: the Dynamic Objective Team Scoring (DOTS), Wilks Coefficient (Wilks), GoodLift Points (Goodlift), and the Glossbrenner Formula (Glossbrenner). These formulas are all well-established within the sport, with some being in use for decades, and their impact on the powerlifting community cannot be overstated.

1.2 Insights from Prior Studies

Over the years, several studies have examined the strengths and weaknesses of these formulas. Research into the Wilks Coefficient has highlighted its widespread adoption, yet also pointed to biases in its body weight adjustment, particularly for lifters at the extremes of body weight. For example, lifters at either very low or very high body weights may see their strength inaccurately represented due to the limitations of the formula’s polynomial adjustments. As a result, newer formulas like DOTS and GoodLift Points have emerged, offering more sophisticated adjustments based on modern statistical analyses and addressing the perceived flaws in the Wilks system. The Glossbrenner formula, which introduces age and additional factors into the equation, are also important for accounting for the performance of older lifters. These formulas offer insights into how different approaches to formula design can affect the rankings and perceptions of athletes across various age groups and weight classes.

1.3 Objectives of Current Study

This study aims to answer the following questions:

  1. How have the different strength ranking systems improved or degraded over time? Are people getting stronger or weaker?

  2. Are the different strength ranking systems staying similar in their representation of the data or is there divergence?

  3. Is it possible to identify a “best” or “most fair” among the evaluated systems?

By analyzing these factors, the study will shed light on the fairness and efficacy of the formulas, as well as their potential for forecasting future results.

1.4 Predictions

Sharing expectations for each of the above questions.

  1. The average powerlifter has gotten significantly stronger over time, with the elite athletes additionally surpassing that.

  2. All strength systems should provide very similar results from analysis, with the most alike being Wilks and Dots. Goodlift and Glossbrenner should be unique but follow similar trends.

  3. “Best” is very subjective, but the system that has the least variance in the extremes of body weight could be a candidate for the “best” system. There also may be a best for average lifters and one for more competitive lifters that value different things

1.5 Formula Explanation

1.5.1 Wilks Coefficient

\[ \text{Wilks Points} = Result \times \frac{500}{A \times Bw^5 + B \times Bw^4 + C \times Bw^3 + D \times Bw^2 + E \times Bw + F} \]

The Wilks Coefficient was the first widely adopted formula in powerlifting and remains in use today. While it has been largely replaced by more modern adaptations, like DOTS, due to known biases and limitations in its adjustment for body weight, it laid the foundation for future scoring systems. The formula, developed by Dr. Robert Wilks in the early 2000s, helped standardize rankings across different weight classes. It’s worth noting that the DOTS formula is essentially an updated version of the Wilks formula, designed to address some of the biases in the original system.

  • Result refers to the combined total of an athletes best squat, bench press, and deadlift.
  • Bw is the athletes body weight.
  • A, B, C, D, E, F coefficients in the Wilks formula that change with the athletes body weight, and are updated as needed in a similar fashion to the DOTS formula.
  • 500 is a constant factor that acts as a means to standardize the formula.

1.5.1.1 Wilks Scoring Reference Table

Wilks Score Classification Description
<=300 Beginner Indicates that you have been lifting for at least a year and have potential for improvement.
350-400 National Level Puts you in the category of national-level competitors.
400+ Championship Level Considered championship-level status.
500+ Elite Puts you in the top 10 in your category.
550+ World-Class Achieved by only a small percentage of lifters, indicating top-tier performance.

1.5.2 DOTS

\[ \text{DOTS Points} = Result \times \frac{500}{A \times Bw^4 + B \times Bw^3 + C \times Bw^2 + D \times Bw + E} \] DOTS is the most modern and widely used formula. DOTS was developed to overcome some of the limitations of older scoring systems like the Wilks Coefficient by using more precise polynomial adjustment based on body weight.

  • Result refers to the combined total of an athletes best squat, bench press, and deadlift.
  • Bw is the athletes body weight.
  • A, B, C, D, E coefficients in the DOTS formula that change based on the athlete’s body weight. These values are determined through statistical analyses of powerlifting data, which are updated when enough new data is gathered to warrant a recalculation or when performance trends significantly diverge from prior projections. These updates typically happen every few years but are not on a fixed schedule.
  • 500 is a constant factor that acts as a means to standardize the formula.

1.5.2.1 DOTS Scoring Reference Table

DOTS Score Classification Description
<300 Beginner Indicates a need for significant improvement.
300-350 Average Represents a basic level of strength.
350-399 Good Shows solid performance; competitive at local levels.
400-449 Strong Considered strong; competitive at regional levels.
450-499 Very Strong Exceptional performance; competitive at national levels.
500+ Elite Outstanding performance; competitive at international levels.

1.5.3 GoodLift Points

\[ \text{GL Points} = Result \times \frac{100}{A - B \times e^{-C\times Bw}} \]

GLP is a modern formula developed by the Goodlift platform.The formula aims to improve on older systems like Wilks by reducing bias and offering a more balanced scoring system. Goodlift Points is used by the IPF (International Powerlifting Federation) and British Powerlifting, as well as several other federations worldwide. It is unique in that it considers factors such as sex, equipment usage, and type of competition.

  • Result refers to the combined total of an athletes best squat, bench press, and deadlift.
  • Bw is the athletes body weight.
  • A, B, C are constants that are selected based on the category of lift. For example, the constants for Men’s Equipped Bench Press are different than Woman’s Classic Bench Press, etc.
  • e is Euler’s number, the base of the natural logarithm.

1.5.3.1 Goodlift Scoring Reference Table

Goodlift Score Classification Description
<55 Poor Beginner lifter with minimal training experience.
55-75 Average Recreational lifter with some consistent training.
76-95 Good Dedicated lifter competitive at local level.
96-115 Strong Experienced lifter competitive at regional level.
116-140 Very Strong Advanced lifter competitive at national level.
141+ Elite Top-tier lifter competitive at international level.

1.5.4 Glossbrenner Formula

\[ \text{Product Number} = Result \times A \times B \]

The Glossbrenner formula is another scoring system similar to Wilks or DOTS. The formula takes into account body weight and age (unlike Wilks and DOTS) and applies coefficients accordingly.

  • Product Number is the number that is used to compare lifters.
  • Result refers to the combined total of an athletes best squat, bench press, and deadlift.
  • A is the Glossbrenner Bodyweight Coefficient. This constant is selected from a table based on body weight, ranging from 1.32435 at 40kg and 0.47853 at 231.90kg.
  • B is the Master’s Age Multiple, another constant ranging from 1.000 at <=40 years old to 2.050 at >=80 years old.

2 Data and Methods

2.1 Data Source and Description

The data set that was used is the Open Powerlifting data base. This data set is sizable, with 3,550,611 rows, and contains entries dating back to the 1960’s. This analysis will utilize several variables throughout, mostly focusing on the scores from the different scoring systems and date. In the latter end of the analysis, we will look at how different systems relate to different lifts, bodyweight, age, and total lifted.

2.1.1 Initialization Code

# Boilerplate, uncomment install.packages() as necessary.

#install.packages("tidyverse")
#install.packages("lubridate")
#install.packages("patchwork")
#install.packages("corrplot")
#install.packages("ggpubr")
#install.packages("car")
library(tidyverse)
library(lubridate)
library(patchwork)
library(corrplot)
library(ggpubr)
library(car)

if (!exists("raw_data")) {
  raw_data <- read_csv("openpowerlifting-2025-03-01-e67afafd.csv")
}
# print(prettyNum(nrow(raw_data), big.mark = ","))

2.2 Methods

For each of the scoring systems, a similar format of analysis will be run through. Firstly, data will be wrangled from an unmodified version of the data set into a summary of all rows associated with that scoring system; secondly, a summary will be wrangled that excludes anyone not within a point range that is considered “professional”, which in this case means someone at or above a level that is considered competitive at a local level; finally, a year-over-year rate of change summary will be wrangled. After the data is wrangled into these three categories for their respective scoring system, a plethora of visualizations will be generated from them:

Visualization 1: Distribution histogram

Visualization 2: Trend analysis for all rows within that scoring category

Visualization 3: Trend analysis for “professional” rows within that scoring category

Visualization 4: Trend analysis for the maximum recorded score for each year

Visualization 5: A score variability by year chart

Visualization 6: A bar chart showing year-over-year absolute change

Visualizations 2 and 3 utilize the mean score per year, and visualization 4 uses the maximum. Visualizations 2, 3, and 4 have standard deviation charted as well for additional context.

After running through each scoring system, a pair of visualizations will compare the yearly averages and maximums of each system for a way to compare them visually.

A series of correlation matrices will be generated among the scoring systems, all recorded lifts, age, and body weight to build regression models, VIF analyses, and associated visualizations.

Through this plethora of analyses, we can expect to see patterns within each scoring system and gather greater insight into the relationships between scores and different key variables. With the first block of visualizations, we will be able to clearly see an answer to our first and second research questions, which aim to identify whether lifters are getting stronger or weaker, by how much, and if all of the systems agree. The comparative analyses, correlation matrices, and regression models will directly address our final research question: Is it possible to identify a “best” or “most fair” among the evaluated systems? By examining how each system correlates with body weight and total weight lifted, we can identify potential biases where certain weight classes or strength levels might be systematically advantaged or disadvantaged.

3 Results

3.1 Wilks

3.1.1 Wilks Score Wrangling

# Wilks Generic Scoring Data
wilks_year_summary_all <- raw_data %>% 
  filter(!is.na(Wilks)) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_wilks = mean(Wilks, na.rm = TRUE),
    max_wilks = max(Wilks, na.rm = TRUE),
    sd_wilks = sd(Wilks, na.rm = TRUE),
    row_count = n()
  )

# Wilks Professional Scoring Data
wilks_year_summary_pro <- raw_data %>%
  filter(Wilks >= 400) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_wilks = mean(Wilks, na.rm = TRUE),
    max_wilks = max(Wilks, na.rm = TRUE),
    sd_wilks = sd(Wilks, na.rm = TRUE),
    row_count = n()
  )

# Wilks YoY Rate of Change
wilks_rate_of_change <- wilks_year_summary_all %>%
  arrange(Year) %>%
  mutate(
    previous_year_avg = lag(avg_wilks),
    absolute_change = avg_wilks - previous_year_avg,
  ) %>%
  filter(!is.na(absolute_change))

3.1.2 Wilks Score Visualizations

# Wilks Visualization 1: Distribution of Wilks scores
ggplot(raw_data %>% filter(!is.na(Wilks)), aes(x = Wilks)) +
  geom_histogram(binwidth = 20, fill = "steelblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Wilks Scores",
       x = "Wilks Score",
       y = "Frequency") +
  theme_minimal()

# Wilks Visualization 2: Trend analysis - All Levels
wilks_all <- ggplot(wilks_year_summary_all, aes(x = Year)) +
  geom_line(aes(y = avg_wilks, color = "Average Wilks"), size = 1) +
  geom_ribbon(aes(ymin = avg_wilks - sd_wilks, ymax = avg_wilks + sd_wilks), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_wilks), color = "steelblue", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_wilks), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Wilks Score",
    sec.axis = sec_axis(
      ~./max(wilks_year_summary_all$avg_wilks)*max(wilks_year_summary_all$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Wilks" = "steelblue")) +
  scale_fill_manual(values = c("Entry Count" = "lightblue")) +
  labs(title = "Wilks Scores and Entries - All Levels",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Wilks Visualization 3: Trend analysis - Professional Levels
wilks_pro <- ggplot(wilks_year_summary_pro, aes(x = Year)) +
  geom_line(aes(y = avg_wilks, color = "Average Wilks"), size = 1) +
  geom_ribbon(aes(ymin = avg_wilks - sd_wilks, ymax = avg_wilks + sd_wilks), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_wilks), color = "darkgreen", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_wilks), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Wilks Score",
    sec.axis = sec_axis(
      ~./max(wilks_year_summary_pro$avg_wilks)*max(wilks_year_summary_pro$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Wilks" = "darkgreen")) +
  scale_fill_manual(values = c("Entry Count" = "lightgreen")) +
  labs(title = "Wilks Scores and Entries - Professional Level (400+)",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Wilks Visualization 4: Max Wilks Score per Year (excluding 2025)
wilks_max <- wilks_year_summary_all %>%
  filter(Year < 2025) %>%
  ggplot(aes(x = Year, y = max_wilks)) +
  geom_line(color = "red", size = 1) +
  geom_ribbon(aes(ymin = avg_wilks - sd_wilks, ymax = avg_wilks + sd_wilks), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(color = "darkred", size = 2) +
  labs(title = "Maximum Wilks Score Per Year (2025 Removed)",
       x = "Year",
       y = "Maximum Wilks Score") +
  theme_minimal()

# Wilks Visualization 5: Wilks Score Variability
wilks_variability <- raw_data %>%
  filter(!is.na(Wilks)) %>%
  mutate(Year = as.factor(year(Date))) %>%
  ggplot(aes(x = Year, y = Wilks, group = Year)) +
  geom_boxplot(fill = "lightgreen", alpha = 0.7) +
  labs(
    title = "Wilks Score Variability by Year",
    x = "Year",
    y = "Wilks Score"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5)
  ) +
  scale_x_discrete(breaks = levels(as.factor(year(raw_data$Date)))[seq(1, length(levels(as.factor(year(raw_data$Date)))), by = 2)])

# Absolute Change YoY for Wilks
wilks_absolute_change <- ggplot(wilks_rate_of_change, aes(x = Year, y = absolute_change)) +
  geom_col(aes(fill = ifelse(absolute_change >= 0, "Positive", "Negative")), width = 0.7) +
  scale_fill_manual(values = c("Positive" = "darkgreen", "Negative" = "darkred")) +
  geom_text(aes(label = sprintf("%.1f", absolute_change), 
                vjust = ifelse(absolute_change >= 0, -0.5, 1.5)),
            size = 3) +
  labs(
    title = "Absolute Change in Average Wilks Scores Year-over-Year",
    x = "Year",
    y = "Change in Score",
    fill = "Direction"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(wilks_all)

print(wilks_pro)

print(wilks_max)

print(wilks_variability)

print(wilks_absolute_change)

3.1.3 Wilks Takeaways

3.1.3.1 Distribution of Wilks Scores

  • Similarly peaked bimodal distribution
  • Potential Causes:
    • Male vs Female
    • Raw vs Equipped
    • Experience Levels
  • First peak of Wilks scores are below 200 which would be considered an extreme beginner
  • Wilks <= 200 @ 200lbs = ~600 Total (generally distributed into Deadlift (~50%) > Squat (~30%) > Bench (20%))
  • Second peak of Wilks scores is between 300-400, which would be considered competitive at regional levels or national levels

3.1.3.2 Average Score per Year, All Levels

  • Extremely wide standard deviation
  • Notable down trend (appears to be stabilizing or increasing)
  • Huge explosion of count in 2010’s

3.1.3.3 Average Score per Year, Professional

  • Consistent within ~50 points for ~60 years
  • Much tighter standard deviation, but widening

3.1.3.4 Max Score per Year

  • Extremely strong up trend for entirety of data set
  • Almost doubled performance
  • Well outside of standard deviation

3.1.3.5 Score Variability

  • Volatile early years
  • Score variability is growing downwards, top side staying fairly stable excluding outliers.

3.1.3.6 Absolute Change

  • Volatile early years
  • 1980 - 2000, tightening range but still fairly volatile
  • 2000 - 2010, downward trend
  • 2010 - 2020+, recovering

3.2 DOTS

3.2.1 DOTS Score Wrangling

# Dots Generic Scoring Data
dots_year_summary_all <- raw_data %>% 
  filter(!is.na(Dots)) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_dots = mean(Dots, na.rm = TRUE),
    max_dots = max(Dots, na.rm = TRUE),
    sd_Dots = sd(Dots, na.rm = TRUE),
    row_count = n()
  )

# Dots Professional Scoring Data
dots_year_summary_pro <- raw_data %>%
  filter(Dots >= 400) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_dots = mean(Dots, na.rm = TRUE),
    max_dots = max(Dots, na.rm = TRUE),
    sd_Dots = sd(Dots, na.rm = TRUE),
    row_count = n()
  )

# Dots YoY Rate of Change
dots_rate_of_change <- dots_year_summary_all %>%
  arrange(Year) %>%
  mutate(
    previous_year_avg = lag(avg_dots),
    absolute_change = avg_dots - previous_year_avg,
  ) %>%
  filter(!is.na(absolute_change))

3.2.2 DOTS Score Visualizations

# Dots Visualization 1: Distribution of Dots scores
ggplot(raw_data %>% filter(!is.na(Dots)), aes(x = Dots)) +
  geom_histogram(binwidth = 20, fill = "steelblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Dots Scores",
       x = "Dots Score",
       y = "Frequency") +
  theme_minimal()

# Dots Visualization 2: Trend analysis - All Levels
dots_all <- ggplot(dots_year_summary_all, aes(x = Year)) +
  geom_line(aes(y = avg_dots, color = "Average Dots"), size = 1) +
  geom_ribbon(aes(ymin = avg_dots - sd_Dots, ymax = avg_dots + sd_Dots), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_dots), color = "steelblue", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_dots), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Dots Score",
    sec.axis = sec_axis(
      ~./max(dots_year_summary_all$avg_dots)*max(dots_year_summary_all$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Dots" = "steelblue")) +
  scale_fill_manual(values = c("Entry Count" = "lightblue")) +
  labs(title = "Dots Scores and Entries - All Levels",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Dots Visualization 3: Trend analysis - Professional Levels
dots_pro <- ggplot(dots_year_summary_pro, aes(x = Year)) +
  geom_line(aes(y = avg_dots, color = "Average Dots"), size = 1) +
  geom_ribbon(aes(ymin = avg_dots - sd_Dots, ymax = avg_dots + sd_Dots), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_dots), color = "darkgreen", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_dots), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Dots Score",
    sec.axis = sec_axis(
      ~./max(dots_year_summary_pro$avg_dots)*max(dots_year_summary_pro$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Dots" = "darkgreen")) +
  scale_fill_manual(values = c("Entry Count" = "lightgreen")) +
  labs(title = "Dots Scores and Entries - Professional Level (400+)",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Dots Visualization 4: Max Dots Score per Year (excluding 2025)
dots_max <- dots_year_summary_all %>%
  filter(Year < 2025) %>%
  ggplot(aes(x = Year, y = max_dots)) +
  geom_line(color = "red", size = 1) +
  geom_ribbon(aes(ymin = avg_dots - sd_Dots, ymax = avg_dots + sd_Dots), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(color = "darkred", size = 2) +
  labs(title = "Maximum Dots Score Per Year (Excluding 2025)",
       x = "Year",
       y = "Maximum Dots Score") +
  theme_minimal()

# Dots Visualization 5: Dots Score Variability
dots_variability <- raw_data %>%
  filter(!is.na(Dots)) %>%
  mutate(Year = as.factor(year(Date))) %>%
  ggplot(aes(x = Year, y = Dots, group = Year)) +
  geom_boxplot(fill = "skyblue", alpha = 0.7) +
  labs(
    title = "Dots Score Variability by Year",
    x = "Year",
    y = "Dots Score"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5)
  ) +
  scale_x_discrete(breaks = levels(as.factor(year(raw_data$Date)))[seq(1, length(levels(as.factor(year(raw_data$Date)))), by = 2)])

# Absolute Change YoY for Dots
dots_absolute_change <- ggplot(dots_rate_of_change, aes(x = Year, y = absolute_change)) +
  geom_col(aes(fill = ifelse(absolute_change >= 0, "Positive", "Negative")), width = 0.7) +
  scale_fill_manual(values = c("Positive" = "darkgreen", "Negative" = "darkred")) +
  geom_text(aes(label = sprintf("%.1f", absolute_change), 
                vjust = ifelse(absolute_change >= 0, -0.5, 1.5)),
            size = 3) +
  labs(
    title = "Absolute Change in Average Dots Scores Year-over-Year",
    x = "Year",
    y = "Change in Score",
    fill = "Direction"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(dots_all)

print(dots_pro)

print(dots_max)

print(dots_variability)

print(dots_absolute_change)

3.2.3 Dots Takeaways

Dots is almost identical to Wilks, which is expected as Dots is an attempt to mimic Wilks but remove or dampen the variance in the extremes of body weight.

3.3 Goodlift

3.3.1 Goodlift Score Wrangling

# Goodlift Generic Scoring Data
goodlift_year_summary_all <- raw_data %>% 
  filter(!is.na(Goodlift)) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_goodlift = mean(Goodlift, na.rm = TRUE),
    max_goodlift = max(Goodlift, na.rm = TRUE),
    sd_goodlift = sd(Goodlift, na.rm = TRUE),
    row_count = n()
  )

# Goodlift Professional Scoring Data
goodlift_year_summary_pro <- raw_data %>%
  filter(Goodlift >= 75) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_goodlift = mean(Goodlift, na.rm = TRUE),
    max_goodlift = max(Goodlift, na.rm = TRUE),
    sd_goodlift = sd(Goodlift, na.rm = TRUE),
    row_count = n()
  )

# Goodlift YoY Rate of Change
goodlift_rate_of_change <- goodlift_year_summary_all %>%
  arrange(Year) %>%
  mutate(
    previous_year_avg = lag(avg_goodlift),
    absolute_change = avg_goodlift - previous_year_avg,
  ) %>%
  filter(!is.na(absolute_change))

3.3.2 Goodlift Score Visualizations

# Goodlift Visualization 1: Distribution of Goodlift scores
ggplot(raw_data %>% filter(!is.na(Goodlift)), aes(x = Goodlift)) +
  geom_histogram(binwidth = 5, fill = "steelblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Goodlift Scores",
       x = "Goodlift Score",
       y = "Frequency") +
  theme_minimal()

# Goodlift Visualization 2: Trend analysis - All Levels
goodlift_all <- ggplot(goodlift_year_summary_all, aes(x = Year)) +
  geom_line(aes(y = avg_goodlift, color = "Average Goodlift"), size = 1) +
  geom_ribbon(aes(ymin = avg_goodlift - sd_goodlift, ymax = avg_goodlift + sd_goodlift), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_goodlift), color = "steelblue", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_goodlift), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Goodlift Score",
    sec.axis = sec_axis(
      ~./max(goodlift_year_summary_all$avg_goodlift)*max(goodlift_year_summary_all$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Goodlift" = "steelblue")) +
  scale_fill_manual(values = c("Entry Count" = "lightblue")) +
  labs(title = "Goodlift Scores and Entries - All Levels",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Goodlift Visualization 3: Trend analysis - Professional Levels
goodlift_pro <- ggplot(goodlift_year_summary_pro, aes(x = Year)) +
  geom_line(aes(y = avg_goodlift, color = "Average Goodlift"), size = 1) +
  geom_ribbon(aes(ymin = avg_goodlift - sd_goodlift, ymax = avg_goodlift + sd_goodlift), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_goodlift), color = "darkgreen", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_goodlift), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Goodlift Score",
    sec.axis = sec_axis(
      ~./max(goodlift_year_summary_pro$avg_goodlift)*max(goodlift_year_summary_pro$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Goodlift" = "darkgreen")) +
  scale_fill_manual(values = c("Entry Count" = "lightgreen")) +
  labs(title = "Goodlift Scores and Entries - Professional Level (75+)",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Goodlift Visualization 4: Max Wilks Score per Year (excluding 2025)
goodlift_max <- goodlift_year_summary_all %>%
  filter(Year < 2025) %>%
  ggplot(aes(x = Year, y = max_goodlift)) +
  geom_ribbon(aes(ymin = avg_goodlift - sd_goodlift, ymax = avg_goodlift + sd_goodlift), 
              alpha = 0.2, fill = "steelblue") +
  geom_line(color = "red", size = 1) +
  geom_point(color = "darkred", size = 2) +
  labs(title = "Maximum Goodlift Score Per Year",
       x = "Year",
       y = "Maximum Goodlift Score") +
  theme_minimal()

# Goodlift Visualization 5: Wilks Score Variability
goodlift_variability <- raw_data %>%
  filter(!is.na(Goodlift)) %>%
  mutate(Year = as.factor(year(Date))) %>%
  ggplot(aes(x = Year, y = Goodlift, group = Year)) +
  geom_boxplot(fill = "skyblue", alpha = 0.7) +
  labs(
    title = "Goodlift Score Variability by Year",
    x = "Year",
    y = "Goodlift Score"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5)
  ) +
  scale_x_discrete(breaks = levels(as.factor(year(raw_data$Date)))[seq(1, length(levels(as.factor(year(raw_data$Date)))), by = 2)])

# Absolute Change YoY for Goodlift
goodlift_absolute_change <- ggplot(goodlift_rate_of_change, aes(x = Year, y = absolute_change)) +
  geom_col(aes(fill = ifelse(absolute_change >= 0, "Positive", "Negative")), width = 0.7) +
  scale_fill_manual(values = c("Positive" = "darkgreen", "Negative" = "darkred")) +
  geom_text(aes(label = sprintf("%.1f", absolute_change), 
                vjust = ifelse(absolute_change >= 0, -0.5, 1.5)),
            size = 3) +
  labs(
    title = "Absolute Change in Average Goodlift Scores Year-over-Year",
    x = "Year",
    y = "Change in Score",
    fill = "Direction"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(goodlift_all)

print(goodlift_pro)

print(goodlift_max)

print(goodlift_variability)

print(goodlift_absolute_change)

3.3.3 Goodlift Takeaways

3.3.3.1 Distribution of Goodlift Scores

  • Much more reasonable bell curve distribution
  • This may be expected as Goodlift is specifically designed for greater weight class normalization

3.3.3.2 Average Score per Year, All Levels

  • Tighter sd than Wilks / Dots
  • More consistent scoring

3.3.3.3 Average Score per Year, Professional

  • Tighter sd than Wilks / Dots
  • More entries (Top of scale is 50k for Wilks/Dots, 80k for Goodlift)

3.3.3.4 Max Score per Year

  • More “spikey” both up and down than Wilks / Dots

3.3.3.5 Score Variability

  • Much more consistent
  • Outliers are closer
  • Range is smaller

3.3.3.6 Absolute Change

  • Very spikey in early yeras, then fairly steady

3.4 Glossbrenner

3.4.1 Glossbrenner Score Wrangling

# Glossbrenner Generic Scoring Data
glossbrenner_year_summary_all <- raw_data %>% 
  filter(!is.na(Glossbrenner)) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_glossbrenner = mean(Glossbrenner, na.rm = TRUE),
    max_glossbrenner = max(Glossbrenner, na.rm = TRUE),
    sd_glossbrenner = sd(Glossbrenner, na.rm = TRUE),
    row_count = n()
  )

# Glossbrenner Professional Scoring Data
glossbrenner_year_summary_pro <- raw_data %>%
  filter(Glossbrenner >= 400) %>%
  mutate(Year = year(Date)) %>%
  group_by(Year) %>%
  summarise(
    avg_glossbrenner = mean(Glossbrenner, na.rm = TRUE),
    max_glossbrenner = max(Glossbrenner, na.rm = TRUE),
    sd_glossbrenner = sd(Glossbrenner, na.rm = TRUE),
    row_count = n()
  )

# Glossbrenner YoY Rate of Change
glossbrenner_rate_of_change <- glossbrenner_year_summary_all %>%
  arrange(Year) %>%
  mutate(
    previous_year_avg = lag(avg_glossbrenner),
    absolute_change = avg_glossbrenner - previous_year_avg,
  ) %>%
  filter(!is.na(absolute_change))

3.4.2 Glossbrenner Score Visualizations

# Glossbrenner Visualization 1: Distribution of Glossbrenner scores
ggplot(raw_data %>% filter(!is.na(Glossbrenner)), aes(x = Glossbrenner)) +
  geom_histogram(binwidth = 20, fill = "steelblue", color = "black", alpha = 0.7) +
  labs(title = "Distribution of Glossbrenner Scores",
       x = "Glossbrenner Score",
       y = "Frequency") +
  theme_minimal()

# Glossbrenner Visualization 2: Trend analysis - All Levels
glossbrenner_all <- ggplot(glossbrenner_year_summary_all, aes(x = Year)) +
  geom_line(aes(y = avg_glossbrenner, color = "Average Glossbrenner"), size = 1) +
  geom_ribbon(aes(ymin = avg_glossbrenner - sd_glossbrenner, ymax = avg_glossbrenner + sd_glossbrenner), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_glossbrenner), color = "steelblue", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_glossbrenner), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Glossbrenner Score",
    sec.axis = sec_axis(
      ~./max(glossbrenner_year_summary_all$avg_glossbrenner)*max(glossbrenner_year_summary_all$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Glossbrenner" = "steelblue")) +
  scale_fill_manual(values = c("Entry Count" = "lightblue")) +
  labs(title = "Glossbrenner Scores and Entries - All Levels",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Glossbrenner Visualization 3: Trend analysis - Professional Levels
glossbrenner_pro <- ggplot(glossbrenner_year_summary_pro, aes(x = Year)) +
  geom_line(aes(y = avg_glossbrenner, color = "Average Glossbrenner"), size = 1) +
  geom_ribbon(aes(ymin = avg_glossbrenner - sd_glossbrenner, ymax = avg_glossbrenner + sd_glossbrenner), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(aes(y = avg_glossbrenner), color = "darkgreen", size = 2) +
  geom_col(aes(y = row_count / max(row_count) * max(avg_glossbrenner), fill = "Entry Count"), alpha = 0.3) +
  scale_y_continuous(
    name = "Average Glossbrenner Score",
    sec.axis = sec_axis(
      ~./max(glossbrenner_year_summary_pro$avg_glossbrenner)*max(glossbrenner_year_summary_pro$row_count), 
      name = "Number of Entries"
    )
  ) +
  scale_color_manual(values = c("Average Glossbrenner" = "darkgreen")) +
  scale_fill_manual(values = c("Entry Count" = "lightgreen")) +
  labs(title = "Glossbrenner Scores and Entries - Professional Level (400+)",
       x = "Year") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Glossbrenner Visualization 4: Max Glossbrenner Score per Year (excluding 2025)
glossbrenner_max <- glossbrenner_year_summary_all %>%
  filter(Year < 2025) %>%
  ggplot(aes(x = Year, y = max_glossbrenner)) +
  geom_line(color = "red", size = 1) +
  geom_ribbon(aes(ymin = avg_glossbrenner - sd_glossbrenner, ymax = avg_glossbrenner + sd_glossbrenner), 
              alpha = 0.2, fill = "steelblue") +
  geom_point(color = "darkred", size = 2) +
  labs(title = "Maximum Glossbrenner Score Per Year",
       x = "Year",
       y = "Maximum Glossbrenner Score") +
  theme_minimal()

# Glossbrenner Visualization 5: Glossbrenner Score Variability
glossbrenner_variability <- raw_data %>%
  filter(!is.na(Glossbrenner)) %>%
  mutate(Year = as.factor(year(Date))) %>%
  ggplot(aes(x = Year, y = Glossbrenner, group = Year)) +
  geom_boxplot(fill = "skyblue", alpha = 0.7) +
  labs(
    title = "Glossbrenner Score Variability by Year",
    x = "Year",
    y = "Glossbrenner Score"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    plot.title = element_text(hjust = 0.5)
  ) +
  scale_x_discrete(breaks = levels(as.factor(year(raw_data$Date)))[seq(1, length(levels(as.factor(year(raw_data$Date)))), by = 2)])

# Absolute Change YoY for Glossbrenner
glossbrenner_absolute_change <- ggplot(glossbrenner_rate_of_change, aes(x = Year, y = absolute_change)) +
  geom_col(aes(fill = ifelse(absolute_change >= 0, "Positive", "Negative")), width = 0.7) +
  scale_fill_manual(values = c("Positive" = "darkgreen", "Negative" = "darkred")) +
  geom_text(aes(label = sprintf("%.1f", absolute_change), 
                vjust = ifelse(absolute_change >= 0, -0.5, 1.5)),
            size = 3) +
  labs(
    title = "Absolute Change in Average Glossbrenner Scores Year-over-Year",
    x = "Year",
    y = "Change in Score",
    fill = "Direction"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

print(glossbrenner_all)

print(glossbrenner_pro)

print(glossbrenner_max)

print(glossbrenner_variability)

print(glossbrenner_absolute_change)

3.4.3 Glossbrenner Takeaways

3.4.3.1 Distribution of Goodlift Scores

  • Back to bimodal distribution

3.4.3.2 Average Score per Year, All Levels

  • Very similar to Dots / Wilks, both in trend and standard deviation

3.4.3.3 Average Score per Year, Professional

  • Also similar to Dots and Wilks

3.4.3.4 Max Score per Year

  • Again, similar to Dots and Wilks, although some data points are placed differently

3.4.3.5 Score Variability

  • Score variability is expectedly very similar to Dots and Wilks

3.4.3.6 Absolute Change

  • Absolute change is similar to Dots / Wilks

3.5 Average and Maximum Comparison Across Systems

wilks_summary <- wilks_year_summary_all %>%
  select(Year, avg_score = avg_wilks, max_score = max_wilks, row_count) %>%
  mutate(system = "Wilks")

dots_summary <- dots_year_summary_all %>%
  select(Year, avg_score = avg_dots, max_score = max_dots, row_count) %>%
  mutate(system = "Dots")

goodlift_summary <- goodlift_year_summary_all %>%
  select(Year, avg_score = avg_goodlift, max_score = max_goodlift, row_count) %>%
  mutate(system = "Goodlift")

glossbrenner_summary <- glossbrenner_year_summary_all %>%
  select(Year, avg_score = avg_glossbrenner, max_score = max_glossbrenner, row_count) %>%
  mutate(system = "Glossbrenner")

all_systems <- bind_rows(wilks_summary, dots_summary, goodlift_summary, glossbrenner_summary)

# Normalize scores within each system (min-max scaling to 0-100 range)
normalized_systems <- all_systems %>%
  group_by(system) %>%
  mutate(
    normalized_avg = (avg_score - min(avg_score, na.rm = TRUE)) / 
      (max(avg_score, na.rm = TRUE) - min(avg_score, na.rm = TRUE)) * 100,
    normalized_max = (max_score - min(max_score, na.rm = TRUE)) / 
      (max(max_score, na.rm = TRUE) - min(max_score, na.rm = TRUE)) * 100
  ) %>%
  ungroup()

combined_all <- ggplot(normalized_systems %>% filter(Year < 2025), aes(x = Year, y = normalized_avg, color = system)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  labs(
    title = "Comparison of All Scoring Systems (Normalized)",
    subtitle = "Average scores normalized to 0-100 scale for comparison",
    x = "Year",
    y = "Normalized Average Score",
    color = "Scoring System"
  ) +
  theme_minimal() +
  scale_color_brewer(palette = "Set1") +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

combined_max <- ggplot(normalized_systems %>% filter(Year < 2025), aes(x = Year, y = normalized_max, color = system)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  labs(
    title = "Comparison of Maximum Scores Across All Systems (Normalized)",
    subtitle = "Maximum scores normalized to 0-100 scale for comparison",
    x = "Year",
    y = "Normalized Maximum Score",
    color = "Scoring System"
  ) +
  theme_minimal() +
  scale_color_brewer(palette = "Set1") +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

print(combined_all)

print(combined_max)

3.5.1 Average & Maximum Comparison Take-away

  • Expectedly Dots/Wilks/Glossbrenner are very similar, but Goodlift is fairly different but maintains trend
  • Goodlift better for Average lifter, worse for Max lifter

3.6 Correlation Matrices and Regression Analysis

scoring_systems <- c("Wilks", "Dots", "Goodlift", "Glossbrenner")
squat_columns <- c("Squat1Kg", "Squat2Kg", "Squat3Kg", "Best3SquatKg")
bench_columns <- c("Bench1Kg", "Bench2Kg", "Bench3Kg", "Best3BenchKg")
deadlift_columns <- c("Deadlift1Kg", "Deadlift2Kg", "Deadlift3Kg", "Best3DeadliftKg")

analysis_data <- raw_data %>%
  na.omit() %>%
  select(Age, BodyweightKg, all_of(squat_columns), all_of(bench_columns), 
         all_of(deadlift_columns), all_of(scoring_systems))

model_results <- list()
for (metric in scoring_systems) {
  cat("Correlation Matrix for", metric, ":\n")
  correlation_matrix <- cor(analysis_data[, c(metric, "BodyweightKg", 
                                              "Best3SquatKg", "Best3BenchKg", 
                                              "Best3DeadliftKg", "Age")])
  print(correlation_matrix)

  analysis_data_scaled <- analysis_data %>%
    mutate(across(c(Best3SquatKg, Best3DeadliftKg, Best3BenchKg, BodyweightKg, Age), scale))
  
  formula <- as.formula(paste(metric, "~ Best3SquatKg + Best3DeadliftKg + Best3BenchKg + BodyweightKg + Age"))
  model <- lm(formula, data = analysis_data_scaled)
  
  cat("\n----- Model Summary for", metric, "-----\n")
  model_summary <- summary(model)
  model_results[[metric]] <- model_summary
  print(model_summary)
  
  # Check for multicollinearity
  cat("\n----- Variance Inflation Factor (VIF) for", metric, "-----\n")
  vif_results <- vif(model)
  print(vif_results)
  
  par(mfrow = c(2, 2))
  plot(model, main = paste("Model Diagnostics for", metric))
}
## Correlation Matrix for Wilks :
##                      Wilks BodyweightKg Best3SquatKg Best3BenchKg Best3DeadliftKg         Age
## Wilks            1.0000000    0.1089929    0.6425968   0.46499225       0.6632655 -0.10179190
## BodyweightKg     0.1089929    1.0000000    0.6983792   0.69943957       0.6479281 -0.10817485
## Best3SquatKg     0.6425968    0.6983792    1.0000000   0.86532717       0.9524503 -0.20817479
## Best3BenchKg     0.4649923    0.6994396    0.8653272   1.00000000       0.8525015 -0.06257753
## Best3DeadliftKg  0.6632655    0.6479281    0.9524503   0.85250149       1.0000000 -0.17165566
## Age             -0.1017919   -0.1081749   -0.2081748  -0.06257753      -0.1716557  1.00000000
## 
## ----- Model Summary for Wilks -----
## 
## Call:
## lm(formula = formula, data = analysis_data_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -61.602 -18.449  -5.685  20.260 101.639 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      315.871      4.111  76.827  < 2e-16 ***
## Best3SquatKg      44.830     15.283   2.933  0.00477 ** 
## Best3DeadliftKg   26.757     13.982   1.914  0.06052 .  
## Best3BenchKg     -12.072      9.000  -1.341  0.18495    
## BodyweightKg     -33.762      6.060  -5.571 6.57e-07 ***
## Age                3.898      4.368   0.892  0.37588    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.15 on 59 degrees of freedom
## Multiple R-squared:  0.6677, Adjusted R-squared:  0.6395 
## F-statistic: 23.71 on 5 and 59 DF,  p-value: 5.502e-13
## 
## 
## ----- Variance Inflation Factor (VIF) for Wilks -----
##    Best3SquatKg Best3DeadliftKg    Best3BenchKg    BodyweightKg             Age 
##       13.605158       11.387201        4.717811        2.138948        1.111566

## Correlation Matrix for Dots :
##                        Dots BodyweightKg Best3SquatKg Best3BenchKg Best3DeadliftKg         Age
## Dots             1.00000000   0.09996338    0.6330862   0.44033734       0.6498749 -0.07407893
## BodyweightKg     0.09996338   1.00000000    0.6983792   0.69943957       0.6479281 -0.10817485
## Best3SquatKg     0.63308615   0.69837916    1.0000000   0.86532717       0.9524503 -0.20817479
## Best3BenchKg     0.44033734   0.69943957    0.8653272   1.00000000       0.8525015 -0.06257753
## Best3DeadliftKg  0.64987486   0.64792811    0.9524503   0.85250149       1.0000000 -0.17165566
## Age             -0.07407893  -0.10817485   -0.2081748  -0.06257753      -0.1716557  1.00000000
## 
## ----- Model Summary for Dots -----
## 
## Call:
## lm(formula = formula, data = analysis_data_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -62.979 -20.202  -8.875  17.752 104.996 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      317.215      4.129  76.821  < 2e-16 ***
## Best3SquatKg      50.237     15.350   3.273  0.00178 ** 
## Best3DeadliftKg   25.139     14.043   1.790  0.07855 .  
## Best3BenchKg     -16.669      9.039  -1.844  0.07019 .  
## BodyweightKg     -33.524      6.086  -5.508 8.33e-07 ***
## Age                5.997      4.387   1.367  0.17684    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.29 on 59 degrees of freedom
## Multiple R-squared:  0.6675, Adjusted R-squared:  0.6394 
## F-statistic: 23.69 on 5 and 59 DF,  p-value: 5.579e-13
## 
## 
## ----- Variance Inflation Factor (VIF) for Dots -----
##    Best3SquatKg Best3DeadliftKg    Best3BenchKg    BodyweightKg             Age 
##       13.605158       11.387201        4.717811        2.138948        1.111566

## Correlation Matrix for Goodlift :
##                    Goodlift BodyweightKg Best3SquatKg Best3BenchKg Best3DeadliftKg         Age
## Goodlift         1.00000000    0.1344140    0.6163768   0.41908477       0.6174662 -0.02585276
## BodyweightKg     0.13441397    1.0000000    0.6983792   0.69943957       0.6479281 -0.10817485
## Best3SquatKg     0.61637676    0.6983792    1.0000000   0.86532717       0.9524503 -0.20817479
## Best3BenchKg     0.41908477    0.6994396    0.8653272   1.00000000       0.8525015 -0.06257753
## Best3DeadliftKg  0.61746618    0.6479281    0.9524503   0.85250149       1.0000000 -0.17165566
## Age             -0.02585276   -0.1081749   -0.2081748  -0.06257753      -0.1716557  1.00000000
## 
## ----- Model Summary for Goodlift -----
## 
## Call:
## lm(formula = formula, data = analysis_data_scaled)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.8854  -5.5180  -0.6677   5.4532  21.9667 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      63.8891     0.9965  64.115  < 2e-16 ***
## Best3SquatKg     12.7333     3.7041   3.438  0.00108 ** 
## Best3DeadliftKg   3.6446     3.3888   1.075  0.28654    
## Best3BenchKg     -4.5372     2.1812  -2.080  0.04187 *  
## BodyweightKg     -6.2253     1.4687  -4.239 8.02e-05 ***
## Age               2.0039     1.0588   1.893  0.06332 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.034 on 59 degrees of freedom
## Multiple R-squared:  0.5996, Adjusted R-squared:  0.5656 
## F-statistic: 17.67 on 5 and 59 DF,  p-value: 1.157e-10
## 
## 
## ----- Variance Inflation Factor (VIF) for Goodlift -----
##    Best3SquatKg Best3DeadliftKg    Best3BenchKg    BodyweightKg             Age 
##       13.605158       11.387201        4.717811        2.138948        1.111566

## Correlation Matrix for Glossbrenner :
##                 Glossbrenner BodyweightKg Best3SquatKg Best3BenchKg Best3DeadliftKg         Age
## Glossbrenner       1.0000000    0.1284344    0.7211190   0.58852795       0.7513779 -0.16177699
## BodyweightKg       0.1284344    1.0000000    0.6983792   0.69943957       0.6479281 -0.10817485
## Best3SquatKg       0.7211190    0.6983792    1.0000000   0.86532717       0.9524503 -0.20817479
## Best3BenchKg       0.5885279    0.6994396    0.8653272   1.00000000       0.8525015 -0.06257753
## Best3DeadliftKg    0.7513779    0.6479281    0.9524503   0.85250149       1.0000000 -0.17165566
## Age               -0.1617770   -0.1081749   -0.2081748  -0.06257753      -0.1716557  1.00000000
## 
## ----- Model Summary for Glossbrenner -----
## 
## Call:
## lm(formula = formula, data = analysis_data_scaled)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -39.055 -14.015  -2.434  11.189  68.835 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     291.8182     2.8257 103.272  < 2e-16 ***
## Best3SquatKg     33.8844    10.5038   3.226  0.00205 ** 
## Best3DeadliftKg  28.1946     9.6096   2.934  0.00476 ** 
## Best3BenchKg      4.0625     6.1854   0.657  0.51387    
## BodyweightKg    -38.1001     4.1648  -9.148 6.48e-13 ***
## Age              -0.4396     3.0024  -0.146  0.88410    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 22.78 on 59 degrees of freedom
## Multiple R-squared:  0.8253, Adjusted R-squared:  0.8105 
## F-statistic: 55.74 on 5 and 59 DF,  p-value: < 2.2e-16
## 
## 
## ----- Variance Inflation Factor (VIF) for Glossbrenner -----
##    Best3SquatKg Best3DeadliftKg    Best3BenchKg    BodyweightKg             Age 
##       13.605158       11.387201        4.717811        2.138948        1.111566

3.6.1 Correlation Matrices and Regression Analysis Important Components

  • Glossbrenner has the highest explanatory power (R² = 0.825) of all scoring systems, suggesting it may provide the most comprehensive adjustment for different lifter characteristics. However, its diagnostic plots show slightly more deviation from ideality compared to the other models, particularly in residual patterns and Q-Q normality

  • Wilks and Dots perform very similarly (R² ≈ 0.67) and show nearly identical patterns in their relationships with lift types

  • All scoring systems strongly adjust for bodyweight (significant negative coefficients), fulfilling their primary purpose of leveling the playing field across weight classes

  • Squat performance is consistently the most significant lift predictor across all systems, with deadlift following closely behind. Bench press is either non-significant or negatively associated with scores when controlling for other lifts

  • The Goodlift system uniquely shows a positive relationship with age (marginally significant), potentially giving some advantage to older lifters. It’s diagnostic plots show good normality in residuals but possible mild heteroscedasticity

  • The correlation analysis shows that deadlift performance has the strongest correlation with Glossbrenner scores (r = 0.75), while having moderate to strong correlations with all other scoring systems.

4 Discussion

I think we promptly answered our first two questions, and made a good start on the third. We identified a consistent and recurring trend of year-over-year similarity amongst scores in both all-levels and professional while also highlighting that there are still amazing feats of strength being achieved every year as well. In relation to these questions, I would love to go back through and experiment with things like removing one-off rows (removing people that are not recurrent participants) and seeing what changes. There is likely a lot of fiddling to be done with this because my assumption is now that the data set is being heavily disrupted by a constant flood of new athletes that are competitive enough to think they should participate in a powerlifting meet, but not enough to progress quickly or consistently. Removing these one-off or low impact anomalies could be very interesting.

For the last question, I think we identified some interesting crumbs to follow with the Glossbrenner explanatory power discovery, and the Goodlift system’s positive relationship with age. This could be expanded more or more effort applied to find similarly strange or unexpected relationships and follow the threads they create.

These analyses showed many interesting trends within competitive powerlifting. We went into the analysis with the firm assumption that we would see fairly boring upward trends across all systems, but what we saw was a fairly confusing mix of upward and downward trends without obvious cause. With what we have here, there are several interesting implications for stakeholders that are focused on competive advantage in the way of statistically verified predictors and their relevance. Additionally, Glossbrenner is revealed to be the system with the highest explanatory power, which may be used as a way to further refine other systems or develop new ones with even more granularity.

Outside of the practical uses and competitive lens, the results of this data seems to imply that powerlifting is the most open and equitable it has ever been, and perhaps some of the people that are putting up amazing performances today were the 200 Wilks level lifters we see dragging the average down in years past. Expansion of interest and participation is always good, and I think there is strong proof here that people are coming in droves to participate, despite their capabilities to be truly competitive.

There are several limitations within this analysis, first and foremost I think the multicollinearity can and should be handled in a more nuanced way to provide for more accurate and stable models. Additionally, a more detailed analysis could be done on different elements of the data set, for example, we could analyze federations specifically and their athletes and see what we can find. Finally, I think doing more specific and advanced wrangling could be really beneficial, trying to either identify potentially junk rows or correcting strange errors with federation inconsistency (some federation will subtract a failed lift from the total, for example).

Overall, I think this analysis answered the majority of the research questions while also providing an avenue for future exploration.