Playing with R, ggplot2 and knitr

Over the last month I have been learning R programming language that is completely FREE and used extensively in statistics and data analysis. I have been also playing with ggplot2 which is a great graphics package for R, along with knitr which was used to produce this blog post/report.

In this short demo, besides bragging with my newly acquired skillz, I wanted to give an example for power calculation in uniform (same speed) circular motion for my friend Michael Clarkson from UK. This was a part of our discussion on difference between running at the same velocity on a straight line, and running on that same velocity in the circle (air resistanc disregarded) in terms of Power production. In the simple dot mass model, in which we are going to represent human body with a dot with a mass of m we are going to disregard air resistance, body oscilations and limb movements. According to this simple model, because the velocity is same all the time, the Power output will be zero. Of course, this is not the case in real life.

The goal of this exercise is to play with vector multiplication and to check if the power is different than zero in uniform circular motion.

Let's generate some data in a data.frame called sample

# create time vector that contains our samplings time.
sample <- data.frame(time = seq(from = 0, to = 1, by = 1/1e+05))

# Calculate time difference or dt used for velocity and acceleration
# differentiation
sample$dt <- with(sample, c(NA, diff(time)))

# Create angular velocity of the object in rad/s
sample$angular.velocity <- 10

# Create radius of the circle
sample$radius <- 2

# Create mass of the object
sample$mass <- 10

This how the sample data.frame looks like

str(sample)

## 'data.frame':    100001 obs. of  5 variables:
##  $ time            : num  0e+00 1e-05 2e-05 3e-05 4e-05 5e-05 6e-05 7e-05 8e-05 9e-05 ...
##  $ dt              : num  NA 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 1e-05 ...
##  $ angular.velocity: num  10 10 10 10 10 10 10 10 10 10 ...
##  $ radius          : num  2 2 2 2 2 2 2 2 2 2 ...
##  $ mass            : num  10 10 10 10 10 10 10 10 10 10 ...

From these parameters that are already in sample data.frame we will generate x and y cordinates of the object.

# Calculate X position
sample$x <- with(sample, sin(time * angular.velocity) * radius)

# Calculate Y position
sample$y <- with(sample, cos(time * angular.velocity) * radius)

On the following graph we can se the path of the object

# Load the ggplot2 package
library(ggplot2)

ggplot(sample, aes(x = x, y = y)) + geom_point(color = "steelblue", size = 2, 
    alpha = 1/2) + labs(title = "Motion path") + geom_hline() + geom_vline()

Now we have positional data and dt or time difference between sampling. We will use this to get velocity for x component and y component (velocity is a vector and we will represent it with x and y components). I beleive this is called reverse dynamics (where you estimate velocity, acceleration and force from positional and time data )

velocity_x = dx / dt

Where dx is chnage in x position between two samples (x2 - x1) and dt is change in time (we have already calculated this when designing the data.frame)

# Calculate velocity for x component
sample$vx <- with(sample, c(NA, diff(x))/dt)

# Calculate velocity for y component
sample$vy <- with(sample, c(NA, diff(y))/dt)

We can now plot velocity (x component)

ggplot(sample, aes(x = time, y = vx)) + geom_path(color = "steelblue", size = 2, 
    alpha = 1/2, na.rm = TRUE) + labs(y = "Velocity x") + geom_hline()

Having velocity data we can calculate acceleration and force, where

acceleration_x = dvelocity_x / dt

Force_x = acceleration_x X mass

# Calculate acceleration for x component
sample$ax <- with(sample, c(NA, diff(vx))/dt)

# Calculate acceleration for y component
sample$ay <- with(sample, c(NA, diff(vy))/dt)

# Calculate force for x and y directions
sample$fx <- with(sample, ax * mass)
sample$fy <- with(sample, ay * mass)

One the following picture is the graph of acceleration (x component)

ggplot(sample, aes(x = time, y = ax)) + geom_path(color = "steelblue", size = 2, 
    alpha = 1/2, na.rm = TRUE) + labs(y = "Acceleration x") + geom_hline()

Now it is time to calculate power on x and on y (scalar product)

# Calculate power for x component and y component
sample$px <- with(sample, vx * fx)
sample$py <- with(sample, vy * fy)

Here is the graph of Power (x component)

ggplot(sample, aes(x = time, y = px)) + geom_path(color = "steelblue", size = 2, 
    alpha = 1/2, na.rm = TRUE) + labs(y = "Power x") + geom_hline()

It seems that the power is fluctuating, but what happens when we calculate total power (sum of x and y component power; power is scalar so we can just add them) and visualize it?

sample$total.power <- with(sample, px + py)
ggplot(sample, aes(x = time, y = total.power)) + geom_path(color = "steelblue", 
    size = 2, alpha = 1/2, na.rm = TRUE) + labs(y = "Total Power") + geom_hline()

Because centripetal force is acting on the object, and that force is perpendicular to the velocity vector, there is no work done by this force and hence power is zero. On the graph above you can see that Total Power has some value - this is constant error due sampling frequency. Higher the sample frequency, lower the error (will see the relationship later)

Another way to calculate total power would be to use scalar velocity

# Calculate scalar velocity from angular velocity and radius
sample$v = with(sample, angular.velocity * radius)

# Calculate acceleration from change in scalar velocity divided by dt
sample$a = with(sample, c(NA, diff(v))/dt)

# Calculate force acting on the object
sample$f = with(sample, a * mass)

# Calculate power
sample$p = with(sample, f * v)

# Plot power calculated
ggplot(sample, aes(x = time, y = p)) + geom_path(color = "steelblue", size = 2, 
    alpha = 1/2, na.rm = TRUE) + labs(y = "Total Power") + geom_hline()

Now we can see that the total power is zero (in this simple model, NOT the real life). Let's code a function that calculates the mean total power for sampling frequency. This will involve copying the whole code above, but changing sampling frequency

mean.total.power <- function(frequency = 100) {

    sample <- data.frame(time = seq(from = 0, to = 1, by = 1/frequency))
    sample$dt <- with(sample, c(NA, diff(time)))
    sample$angular.velocity <- 10
    sample$radius <- 2
    sample$mass <- 10
    sample$x <- with(sample, sin(time * angular.velocity) * radius)
    sample$y <- with(sample, cos(time * angular.velocity) * radius)
    sample$vx <- with(sample, c(NA, diff(x))/dt)
    sample$vy <- with(sample, c(NA, diff(y))/dt)
    sample$ax <- with(sample, c(NA, diff(vx))/dt)
    sample$ay <- with(sample, c(NA, diff(vy))/dt)
    sample$fx <- with(sample, ax * mass)
    sample$fy <- with(sample, ay * mass)
    sample$px <- with(sample, vx * fx)
    sample$py <- with(sample, vy * fy)
    sample$total.power <- with(sample, px + py)
    return(mean(sample$total.power, na.rm = TRUE))
}

Now we can see how does the error in Total Power behaves with different sampling frequencies. We are going to create a vector (series of data) that contains sampling frequencies - 10, 100, 1000, etc and pass that to the function. We are going to graph the error against sampling frequency. If you run this on your computer it is going to take some time.

probing.data <- data.frame(frequency = seq(from = 500, to = 1e+05, by = 1000))  # Different sampling frequencies
probing.data$power = with(probing.data, numeric(length(frequency)))  # This contains the results

for (i in seq_along(probing.data$frequency)) {
    probing.data$power[i] <- mean.total.power(probing.data$frequency[i])
}

ggplot(probing.data, aes(x = frequency, y = power)) + geom_path(color = "steelblue", 
    size = 2, alpha = 1/2, na.rm = TRUE) + labs(y = "Total Power") + geom_hline()

If we do this on other indices, like velocity and acceleration we won't see such a change based on sampling frequency. Let's try - we just need to modify the written function to return the maximum of acceleration on x.

max.x.acceleration <- function(frequency = 100) {

    sample <- data.frame(time = seq(from = 0, to = 1, by = 1/frequency))
    sample$dt <- with(sample, c(NA, diff(time)))
    sample$angular.velocity <- 10
    sample$radius <- 2
    sample$mass <- 10
    sample$x <- with(sample, sin(time * angular.velocity) * radius)
    sample$y <- with(sample, cos(time * angular.velocity) * radius)
    sample$vx <- with(sample, c(NA, diff(x))/dt)
    sample$vy <- with(sample, c(NA, diff(y))/dt)
    sample$ax <- with(sample, c(NA, diff(vx))/dt)
    sample$ay <- with(sample, c(NA, diff(vy))/dt)
    return(max(sample$ax, na.rm = TRUE))
}

And code the same probing procedure as we did for error in total power calculus

probing.data <- data.frame(frequency = seq(from = 500, to = 1e+05, by = 1000))  # Different sampling frequencies
probing.data$acceleration = with(probing.data, numeric(length(frequency)))  # This contains the results

for (i in seq_along(probing.data$frequency)) {
    probing.data$acceleration[i] <- max.x.acceleration(probing.data$frequency[i])
}

ggplot(probing.data, aes(x = frequency, y = acceleration)) + geom_path(color = "steelblue", 
    size = 2, alpha = 1/2, na.rm = TRUE) + labs(y = "Max calculated acceleration") + 
    geom_hline()

As can be seen from the plot, maximal acceleretaion (x comp) is always the same regardles of the sampling frequency.

I hope that this demonstration of R, ggplot and knitr and my simulation knowledge didn't cause major headache.

I will be using R more and more and my next project is something that will be of interest for sport scientists interested in metabolic power measurement using GPS with the goal of real-time tracking of fatigue.

Stay tuned….

Read More

Interview with Chris Carling

Interview with Chris Carling

I have following work by Chris Carling for last couple of years now and one of his recent papers was a staple in my RSA is overrate article. Chris also wrote two books: Performance Assessment for Field Sports and Handbook of Soccer Match Analysis: A Systematic Approach to Improving Performance.

Chris latest article entitled Interpreting Physical Performance in Professional Soccer Match-Play: Should We be More Pragmatic in Our Approach? and is very important and needed.

Chris was kind enough to take some of his free time to answer my questions that might interest a lot of readers, especially those working in team sports, especially soccer.

Mladen: Chris, thank you very much for taking your time to do this interview.  Before I start picking your brain on some topics can you please share to the readers who you are, what do you currently do and where, along with our professional interests.

Chris: I have a BSC degree in Sports Science from Liverpool John Moors University & a PhD in Sports Science from the University of Central Lancashire. I currently work as a Sports Scientist for Lille FC (French Ligue 1) and am Senior Research Fellow in Sports Science at the University of Central Lancashire. I previously worked on the AMISCO Pro game analysis system and conducted research for the Clairefontaine National Football Centre in France.

Mladen: Let’s nail this daemon first since I believe it is very important and often misunderstood by coaches. When coaches read scientific papers (or only abstracts which is even worse) and see significant improvements or differences they usually think of large differences or improvements between treatments and/or groups. What they don’t get is that this statistical significance researchers refer to is probability of making Type I errorand has nothing to do with magnitudes of effects. Thus coaches tend to jump (or not jump) to conclusions based on statistical significance, while they actually think of practical significance.

One example might be using the difference in time motion analysis between positions to make position-specific conditioning based on statistical significance in distance run, while that distance might be only couple of meters and bear no practical significance. What is your take on this and how we can bridge this gap between researchers and coaches? Do we need new statistical approach (i.e. magnitude based statistics)?

Chris: I honestly think that even practical significance type statistics can be misleading. Why? Because differences or changes in data hence performance need to be placed in the real-world context of professional football. Coaches interpret differences in their own way, according to what they might or might not expect, in the context of current form and the quality of the players they have or don’t have at their disposal. Even a difference that is considered low in practical significance and non significant can be considered positive as it might mean that while a team has not improved, its performance has stabilized especially as recent games were against higher standard opposition for example.

Mladen: Continuing with previous question, one way to make decisions based on data might be to know smallest worthwhile change (SWC) and typical error (TE) of estimate. Since the game related performance tend to vary a lot between games for an individual (up to 30% CV), TE usually gets a lot higher than SWC which doesn’t make game related performance good test per se, right? Couple of recent paper stated that teams finishing higher in rank tend to run less than teams finishing lowest in rank. I wonder would that data have any practical significance if viewed with SWC/TE lens?

Chris: I think we need to relate the stats to the type of data we work with – perhaps SWC might be more suited to interpreting changes for example in Repeated sprint test ability (e.g., mean time) after a training intervention rather than match and time motion analysis type data that vary greatly and naturally depend upon many factors that simply cannot be controlled for. I personally prefer simple descriptive statistics (means, totals, percentage changes and differences in these) and in my experience these speak more to practitioners who can attempt to interpret the drop or improvement or even lack of change again according to the context the data were collected in.

Mladen: We touched a bit on the reliability of data with previous question, lets deal with validity for a moment.  Clubs and researchers tend to use GPS data more and more (which is great), but I wonder how much that data is really representing what is happening on the field. Are we missing a lot by only taking velocity into consideration? For example if a player make quick burst for 2m towards the opponent from standing still he won’t reach higher velocity zones for that action to be classified as high-intensity although his power output might be tremendous. Roberto Colli, Osgnach and di Prampero wrote about using power zones instead of velocity zones for this sake (see the translationof one of the articles). What is your take on this and do you think this approach might yield some practical significance between positions, players and levels of play?

Chris: This is an area of research currently being explored in various clubs across the world. Yes the data could be useful to determine the position specific loads experienced in match play but in my opinion we need to take a hard look at the practical usefulness of the data in training and preparation for competition. If differences are observed, this means these already exist and that the player is capable of doing them anyway! Once could say that performing more of this specific training might be useful in developing a players ability to accelerate quicker or perform more of these actions. However, will genetic limitations limit a player’s capacity to improve anyway and the tactical requirements of his/her position might mean that there is no need to perform more of these actions anyway! Running more doesn’t mean a better ability to score or prevent goals which are the two main aims of soccer.

Mladen: Recently coaches started evaluating and training Repeat Sprint Ability (RSA) more and more. Couple of research papers including yours showed that Repeat Sprint Sequence (RSS) doesn’t happen that often in a game, thus decreasing its importance. Do you think these results might change when power-based time motion analysis might be used instead of velocity-based one? Also, how misguided is to rely on averages in the analysis (e.g. RSS happens 1.1 times per game per player on average) while neglecting distributions and worst case scenarios. Can you please expand more on this along with what might be the worst case scenario for certain position in a game from the data you have? What might RSA training give us in terms of game transfer if there is not much RSS happening in a game?

Chris: Results will always depend on the definition of a repeated sprint sequence, i.e. duration of each individual sprints, how many, over how long etc. RSS determined individually according to metabolic power thresholds for example might be useful though and should be explored to see whether the RSS demands are actually higher than demonstrated in our study. In our data, even the players (fulbacks) who performed the most RSA performed (1.7) about 0.6 actions more per game than CD, the mean & SD across all positions were only 1.1  / 1.1. Specific RSA training has also been shown to help other physiological aspects(VO2max)  so should not be ruled out entirely, but as RSS match data apparently demonstrate that this specific quality is not as important as one might  think then practitioners should reflect on the real world usefulness of implementing such training until we provide power based RSS data.

  

Mladen: A lot of pro clubs track GPS and Acceleration data as a form of evaluating training load. What is common practice is to use absolute velocity zones to evaluate training/game load. Do you believe that using relative intensity zones (for example using individual’s vLT, MAS, v30-15IFT and VMAX) might yield more valid data to keep track of workloads for a given individual? Expressed this way, do you believe that it might help preventing overtraining and/or injury?

Chris: Some practitioners adapt their training data according to personalized sprint speed thresholds for each player which allows a more objective determination of training loading. MAS vLT, one should remember are often determined using continuous linear running protocols that do not really represent the actual physical demands of the game – ie the intermittent running activity profile. These data are definitely useful for monitoring players but require quite a lot of expense (buying enough systems for a squad of players) and human interpretation of vast amounts of information. I recently read an interview with Sir Alex Ferguson who said that he could detect when a player was carrying an injury when the player actually thought he was ok, thus one could say we need some subjective analysis in there too!

Mladen: And for the last question, what is your opinion in using ‘efficiency’ scores? For example instead of only tracking physical performance data one might use both physical and technical/tactical data and combine them: dividing amount of high velocity distance by number of successful passes or some other technical or tactical statistic. Do you think that this might give us more power in evaluating players, clubs, leagues? Also, what about ‘efficiency’ score comparing internal vs. external load: dividing high velocity distance by iTRIMP score or time >90%HRmax? From one source of mine tracking these over time for a given player might give some insights into overtraining and injury potential. What are your thoughts on these? 

Chris: In my club, we use efficiency scores mainly for technical scores ie ratios of shots to goals, possessions to chances created… Problem is the weighting of ratings, do we give equal weightings to physical and technical performance for example or should these be adapted to League position – top teams tend to run less so should we be concentrating on technical parameters whereas lower teams might rely more on physical ability. Teams that are strong in one or the other might simply end up being balanced out and having similar ratings. For HR data, the moment we are somewhat limited by the rules of the game, ie we cannot collect in competition. We can do all the predictions we want using physiological/physical data but many injuries are down to contact situations that the player can do nothing about, also we should not forget that some coaches know their players well enough to detect when there is an issue (see earlier comment). Most managers are clever enough to rotate their team (where possible) to keep players as fresh as possible. Simple, subjective ratings from players (after training and/or match-play) are an easy, cheap and reliable means of keeping track of monitoring players.

Mladen: Thank you very much for sharing these invaluable insights Chris.  My readers and me appreciate your time and energy for doing this interview. I wish you all the best in your future endeavors and I am looking forward to new insights from your research.

Read More