Individualized Athlete Monitoring in R, Part 1

Performing SPC in R

It’s always easier to understand this stuff with visuals. If you’d like to follow along, the data can be found here. Let’s see what we’re working with:

These are real flight time values measured on switch mats, although the athletes’ names have been randomized. The data come from 21 vertical jump testing sessions across a fall season and include 1 “baseline” testing session and 20 in-season testing sessions that were carried out ~ 4 hours before most of the season’s matches. Each athlete completed 2 - 3 trials at each testing session, and the two closest trials were stored in the monitoring database.

First things first, we need to obtain session averages for each athlete. I know dplyr tends to be the go-to for data manipulation, so I’ve included both dplyr and data.table (my personal favorite) versions of the code throughout the remainder of the post. One of these days, I’ll get around to discussing the usefulness of using dplyr and data.table together…one day.

# Importing the data
jump.data <- read.csv('jump-data.csv')

# We need to convert the dates from factor to date type
jump.data$jumpDate <- as.Date(jump.data$jumpDate)

# dplyr syntax
dplyr.jump.means <- jump.data %>% group_by(jumpDate, athlete) %>%
  summarise(avg_ft = mean(flightTime))

# data.table syntax
jump.data <- data.table(jump.data)

dt.jump.means <- jump.data[, .(avg_ft = mean(flightTime)), by = .(jumpDate, athlete)]

In either case, we now have mean session values for each athlete:

We need a few more pieces of information to create control charts for our analysis–namely, each athlete’s overall mean flight time and the standard deviation of their flight times:

# dplyr syntax
dplyr.jump.means <- dplyr.jump.means %>% group_by(athlete) %>% 
  mutate(ind_mean = round(mean(avg_ft), 3),
         ind_sd = round(sd(avg_ft), 3))

# data.table syntax
dt.jump.means[, ':=' (ind_mean = round(mean(avg_ft), 3), 
                      ind_sd = round(sd(avg_ft), 3)), by = athlete]

knitr::kable(dt.jump.means[1:10,])

Once we have the athletes’ means and standard deviations, we can calculate their upper and lower control limits. In this example, I’m using limits of mean ± 1.5 SD:

# dplyr syntax
dplyr.jump.means <- dplyr.jump.means %>% group_by(athlete) %>%
  mutate(upper_limit = round(ind_mean + 1.5 * ind_sd, 3),
         lower_limit = round(ind_mean - 1.5 * ind_sd, 3)) %>% ungroup()

# data.table syntax
dt.jump.means[, ':=' (upper_limit = round(ind_mean + 1.5 * ind_sd, 3),
                      lower_limit = round(ind_mean - 1.5 * ind_sd, 3)),
              by = athlete]

knitr::kable(dt.jump.means[1:10,])

Now we have everything we need to plot control charts for our athletes. To start, I’m only going to plot the data from a single athlete because things can get cluttered pretty quickly.

# Plot a single athlete's data
ggplot(dt.jump.means[athlete == 'Steve Galligan'], aes(x = jumpDate, y = avg_ft)) + 
  geom_point() + geom_line() + 
  geom_hline(aes(yintercept = ind_mean)) + 
  geom_hline(aes(yintercept = upper_limit), linetype = 'dotdash') +
  geom_hline(aes(yintercept = lower_limit), linetype = 'dotdash') + theme_bw() +
  labs(title = 'Jumps Over Time', x = 'Testing Date', y = 'Flight Time (s)')

Now we have a very basic control chart worked out. Obviously, we’re examining the data retrospectively for an athlete we managed relatively well so a majority of his flight time values fall within his control limits (the drop you see at the end was after Friday-Sunday double overtime games where he played 110 minutes each game). Examining the plot brings up an important issue with SPC, though: process shifts. If you take a look at the first few testing sessions of the season, there is a noticeable increase in his performance leading into September. From knowing this athlete, I can tell you he came in unfit and having performed zero strength training that summer so his increased jump performance was to be expected. It’s safe to assume that we’ve captured his fitness changes in our monitoring data, and these fitness changes have created a process shift in his vertical jump performance. Not accounting for this shift will artificially inflate the control limits, which could cause us to miss otherwise important decreases in his monitoring data.

So how can we deal with process shifts? Your chosen solution will ultimately boil down to how much free time you have, but I’ll wager most of you reading don’t have a lot. One of the easiest solutions is to use a rolling window approach to calculate your individual means and standard deviations. Like most things in sports science, there are no rules for rolling window size, but I would say whatever your window size is for your chronic rolling window (if you’re using A:C calculations) is a good place to start. For us, that was a 28-day window. So let’s reacalculate our limits and plot the data for the final 28 days of the season:

# dplyr syntax using "chaining"
dplyr.recent.jumps <- dplyr.jump.means %>% filter(jumpDate >= max(jumpDate) - 28) %>%
  group_by(athlete) %>% mutate(ind_mean = mean(avg_ft),
                               ind_sd = sd(avg_ft),
                               upper_limit = ind_mean + 1.5 * ind_sd,
                               lower_limit = ind_mean - 1.5 * ind_sd)

# data.table syntax using chaining; data.table functions a little differently
dt.recent.jumps <- dt.jump.means[jumpDate >= max(jumpDate) - 28]

dt.recent.jumps[, ':=' (ind_mean = mean(avg_ft),
                        ind_sd = sd(avg_ft)),
                by = athlete][, ':=' (upper_limit = ind_mean + 1.5 * ind_sd,
                                      lower_limit = ind_mean - 1.5 * ind_sd),
                              by = athlete]

ggplot(dt.recent.jumps[athlete == 'Steve Galligan'], aes(x = jumpDate, y = avg_ft)) + 
  geom_point() + geom_line() + 
  geom_hline(aes(yintercept = ind_mean)) + 
  geom_hline(aes(yintercept = upper_limit), linetype = 'dotdash') +
  geom_hline(aes(yintercept = lower_limit), linetype = 'dotdash') + theme_bw() +
  labs(title = 'Jumps Over Time', x = 'Testing Date', y = 'Flight Time (s)')

In both cases, the final testing session of the season fell below the athlete’s control limits, so both roads would have led to Rome. While this approach worked in the current case, you may have to adjust your window based on the frequency of your testing. For the year in question, we tested every ~ 4 days, but the next season was closer to every 10 days. Implementing a 28-day rolling window for both seasons would provide drastically different numbers of sessions for our calculations. In the case of a college sport, an alternative may be removal of the pre-season from the calculations. Ultimately, though, this will be an individualized decision based on your specific testing and monitoring setup.

Getting Fancy

There are no hard and fast rules for what threshold constitutes an out of control process. For instance, a threshold of 1 SD (.16 probability of occurring by random chance) may be too reactive to otherwise harmless fluctuations in performance while a threshold of 2 SD (.025 probability of occurring by random chance) may be too strict [5,6]. Personally, I prefer a threshold of 1.5 SD (.065 probability of occurring by random chance), but other methods such as setting an “alert” threshold and a “warning” threshold (e.g. 1 SD = alert and 1.5 or 2 SD = warning) or implementing more sophisticated rules (see Wikipedia for a discussion on different rule sets) are also possible [6].

For example, let’s apply alert (1 SD) and warning (1.5 SD) thresholds to our data from the first example.

# dplyr syntax
dplyr.jump.means <- dplyr.jump.means %>% group_by(athlete) %>%
  rename(upper_limit_1.5SD = upper_limit, lower_limit_1.5SD = lower_limit) %>%
  mutate(upper_limit_1SD = round(ind_mean + ind_sd, 3),
         lower_limit_1SD = round(ind_mean - ind_sd, 3)) %>% ungroup()

# data.table syntax
setnames(dt.jump.means,
         old = c('upper_limit', 'lower_limit'),
         new = c('upper_limit_1.5SD', 'lower_limit_1.5SD'))

dt.jump.means[, ':=' (upper_limit_1SD = round(ind_mean + ind_sd, 3),
                      lower_limit_1SD = round(ind_mean - ind_sd, 3)
                      ), by = athlete]

ggplot(dt.jump.means[athlete == 'Steve Galligan'], aes(x = jumpDate, y = avg_ft)) +
  geom_point() +
  geom_line() +
  geom_hline(aes(yintercept = ind_mean)) +
  geom_hline(aes(yintercept = lower_limit_1.5SD)) +
  geom_hline(aes(yintercept = upper_limit_1.5SD)) +
  geom_hline(aes(yintercept = lower_limit_1SD), linetype = 'dashed') +
  geom_hline(aes(yintercept = upper_limit_1SD), linetype = 'dashed') +
  theme_bw() +
  labs(title = 'Jumps Over Time', x = 'Testing Date', y = 'Flight Time (s)')

It’s also possible to color code the plots based on whether they violate your chosen thresholds or more advanced rule sets. The code gets very messy (I originally included example code but removed it because it was difficult to follow), so if you want to use color coding or more advanced rule sets, I would seriously recommend exploring one of the specially-built SPC packages available on CRAN. qcc, for instance, can create both individual control charts (“I” charts) and EWMA-based charts. Selection on chart type will depend on the amount of data you have and the window size you’ve chosen.