Introduction to Samozino et al. method

# Install vjsim if you haven't already by running the following commands
# install.packages("devtools")
# devtools::install_github("mladenjovanovic/vjsim")

# Install tidyverse and cowplot packages
# install.packages(c("tidyverse", "cowplot", "DT), dependencies = TRUE)

library(vjsim)
library(tidyverse)
library(cowplot)
library(DT)

Before reading this vignette, please read Introduction to vjsim, Simulation in vjsim and Profiling in vjsim vignettes by running:

vignette("introduction-vjsim")
vignette("simulation-vjsim")
vignette("profiling-vjsim")

This vignette is influenced with the works of Samozino, JB Morin et al. (Samozino et al. 2008, 2010, 2012, 2013; Pierre Samozino 2018b, 2018a, 2018a) and Jimenez-Reyes et al. (Jiménez-Reyes et al. 2017; Jiménez-Reyes, Samozino, and Morin 2019). I will provide a short summary and implementation in vjsim as well as some of my opinions regarding Samozino optimization model, but you are urged to check the references for more thorough information.

In the previous vignette Profiling in vjsim I demonstrated the ideas behind Force-Velocity (FV) and Load-Velocity (LV) profiling and their implementation in vjsim, as well as the issues with those methods. This is relatively easy, since we are using simulated data. But collecting data on real athletes is harder, and Samozino et al. (2008) and Samozino et al. (2010) created a simple method for estimating mean force and mean velocity from load~jump height profiling. To demonstrate how this works, and what are the assumptions of this model, I will used the profile from the previous vignette:

# External load in kilograms
external_load <- c(0, 20, 40, 60, 80, 100)

profile_data <- vjsim::vj_profile(
  external_load = external_load,

  # Simulation parameters
  mass = 75,
  push_off_distance = 0.4,
  max_force = 3000,
  max_velocity = 4

  # Other parameters are default in the `vjsim::fgen_get_output`
  # weight = mass * gravity_const,
  # gravity_const = 9.81,
  # decline_rate = 1.05,
  # peak_location = -push_off_distance * 0.15,
  # time_to_max_activation = 0.3
)

As illustrated in the previous vignette, with vjsim we can create various FV profiles. But Samozino et al. (2008) and Samozino et al. (2010) suggested a simple field method that can approximate mean force (\(\bar{F}\)) and mean velocity (\(\bar{v}\)) from a known (1) bodyweight, (2) jump height and (3) push-off distance. Otherwise, more expensive equipment, such as force place, is needed. The goal is to approximate laboratory Force-Velocity profiling using field Load-Height profile.

To represent jump, we will use introduced model in vjsim, where we have object of mass \(m\), that is accelerated over push-off distance (or height) \(h_{po}\) by the action of the Force Generator, eventually reaching jump height \(h\). The jump is performed vertically, on Earth, where the gravity constant \(g\) is equal to \(9.81ms^-2\). At the end of the push-off distance, object reaches the take-off velocity (\(TOV\), or \(v_{to}\)), from which we can calculate height reached. Actually, in this field method, we calculate TOV from jump height reached using the inverse ballistic equation:

\[ v_{to} = \sqrt{2gh} \]

For the object to reach jump height, Force Generator must do some work. Using the work (or energy) equation, we can derive mean force \(\bar{F}\). Please note that this mean force is mean over distance (see Introduction to vjsim). Work done by the Force Generator is equal to weight of the object (i.e. \(mg\)), multiplied by total distance, or sum of push-off distance and height reached (i.e. \(h_{po}+h\)):

\[ W = mg(h_{po}+h) \]

The above equation uses potential energy. We can also use potential plus kinetic energy form, where we sum the potential energy acquired by lifting the object for \(h_{po}\) and kinetic energy using \(v_{to}\). These two equations are identical (actually this is how ballistic equation to get height from \(v_{to}\) and vice versa is derived in the first place).

\[ W = mgh + \frac{mv_{to}^2}{2} \]

Since work is equal \(F \times d\), and Force Generator acts only over the \(h_{po}\), we get the following equation for mean force \(\bar{F}\):

\[ \begin{aligned} W &= \bar{F}\times h_{po} \\ \bar{F} &= \frac{W}{h_{po}} \\ \bar{F} &= \frac{mg(h_{po}+h)}{h_{po}} \\ \bar{F} &= mg(\frac{h}{h_{po}}+1) \end{aligned} \]

Now we have the mean force acting over \(h_{po}\). But, how do we get the mean velocity \(\bar{v}\) over \(h_{po}\) (which is needed for FV profiling)?

The mean velocity over \(h_{po}\) is equal to \(h_{po}\) divided by the time it takes to cover it, or \(t_{po}\):

\[ \bar{v} = \frac{h_{po}}{t_{po}} \]

Unfortunately, we do not know \(t_{po}\) (using this simple field approach). What Samozino et al. (2008) and Samozino et al. (2010) have done (see Samozino et al. (2008) page 2941, equation 6) is to use impulse-momentum equation, which states that impulse (\(I = (\bar{F} - mg) \times t_{po}\)) is equal to the momentum (\(M = m \times v_{to}\)). We thus get the following equation, from which we can extract \(t_{po}\):

\[ \begin{aligned} I &= M \\ (\bar{F} - mg) \times t_{po} &= m \times v_{to} \\ t_{po} &= \frac{m \times v_{to}}{\bar{F} - mg} \end{aligned} \]

If we use the equation that \(v_{to} = \sqrt{2gh}\), we get the following equation for \(\bar{v}\):

\[ \begin{aligned} \bar{v} &= \frac{h_{po}}{t_{po}} \\ \bar{v} &= \frac{h_{po}}{\frac{m \times v_{to}}{\bar{F} - mg}} \\ \bar{v} &= \frac{h_{po}}{\frac{m \times v_{to}}{mg(\frac{h}{h_{po}}+1) - mg}} \\ \bar{v} &= \frac{h_{po}}{\frac{m \times v_{to}}{\frac{mgh}{h_{po}}}} \\ \bar{v} &= \frac{gh}{v_{to}} \\ \end{aligned} \]

We can now use \(v_{to} = \sqrt{2gh}\) to solve this equation:

\[ \begin{aligned} \bar{v} &= \frac{gh}{v_{to}} \\ \bar{v} &= \frac{gh}{\sqrt{2gh}} \\ \bar{v} &= \sqrt{\frac{gh}{2}} \end{aligned} \]

This also implies that \(\frac{v_{to}}{\bar{v}}\) is equal to:

\[ \begin{aligned} \frac{v_{to}}{\bar{v}} &= \\ \frac{\sqrt{2gh}}{\sqrt{\frac{gh}{2}}} &= \\ &= 2 \end{aligned} \]

Why is this? Because, equations used by Samozino et al. (2008) and Samozino et al. (2010) assumes same \(\bar{F}\) over distance (estimated using Work/Energy equations) and \(\bar{F}\) over time (estimated using Impulse - Momentum equations). This assumption is only valid when there is constant force applied and there is uniform acceleration. In that case (assuming starting velocity to be equal to 0 m/s), \(v_{to}\) is twice as big as \(\bar{v}\).

In short, field method proposed by Samozino et al. (2008) and Samozino et al. (2010) approximates \(\bar{F}\) and \(\bar{v}\) assuming constant force and is uniform acceleration across \(h_{po}\). Although wrong, the error introduced by this assumption is estimated in the aforementioned studies when compared with laboratory measurement and it doesn’t seem to be practically significant. But this error (or assumption) simplifies the estimation of \(\bar{F}\) and \(\bar{v}\) in the field settings.

In vjsim this method of profiling is implemented using true \(\bar{v}\) and mean GRF over distance (\(\bar{F}\)), which I called theoretical (and thus represent Mean Force ~ Mean Velocity profile) and using \(\frac{1}{2}v_{to}\) as \(\bar{v}\), which I called practical. We will explore their similarities and differences in this and the next vignette. Here is how these are calculated in vjsim (focus on \(F_0\) and \(V_0\) and disregard optimal metrics for now - they will soon be explained):

samozino_profiles <- vjsim::get_all_samozino_profiles(profile_data)


datatable(samozino_profiles$data_frame, rownames = FALSE) %>%
  formatRound(columns = 1:ncol(samozino_profiles$data_frame), digits = 2)

Here is the plot of these two profiles:

plot_data <- data.frame(
  profile = rep(c("Theoretical", "Practical"), each = 2),
  group = rep(1:2, each = 2),
  x = c(samozino_profiles$list$samozino_theoretical_profile$F0, 0, samozino_profiles$list$samozino_practical_profile$F0, 0),
  y = c(0, samozino_profiles$list$samozino_theoretical_profile$V0, 0, samozino_profiles$list$samozino_practical_profile$V0)
)

plot_data_points <- rbind(
  data.frame(
    profile = "Theoretical",
    group = 1,
    x = profile_data$mean_GRF_over_distance,
    y = profile_data$mean_velocity
  ),
  data.frame(
    profile = "Practical",
    group = 2,
    x = profile_data$mean_GRF_over_distance,
    y = profile_data$mean_velocity_as_TOV_half
  )
)

ggplot(plot_data, aes(x = x, y = y, group = group, color = profile)) +
  theme_cowplot(8) +
  geom_line(alpha = 0.8) +
  geom_point(data = plot_data_points, alpha = 0.8) +
  labs(x = "Mean GRF over distance (N)", y = "Mean Velocity (m/s)")

Effect of \(\bar{F}\) on jump height

If mean force \(\bar{F}\) acts on the body of mass \(m\), weight \(w=mg\) over a push-off distance \(h_{po}\), starting with initial velocity \(v=0\), then the final take-off velocity reached \(v_{to}\) is equal to:

\[ \begin{aligned} (\bar{F} - mg)\times h_{po} &= \frac{v_{to}^2}{2} \\ v_{to} &= \sqrt{\frac{2h_{po}(\bar{F} - mg)}{m}} \end{aligned} \]

Assuming athlete weights 75kg, we get the following graph for varying \(h_{po}\) and \(\bar{F}\)

# Need this package to label the lines
# install.packages("directlabels")
require(directlabels)

# Function to calculate TOV
get_tov <- function(mass, gravity_const, mean_force, push_off_distance) {
  sqrt((2 * push_off_distance * (mean_force - mass * gravity_const)) / mass)
}

params <- expand_grid(
  mass = 75,
  gravity_const = 9.81,
  push_off_distance = c(0.2, 0.3, 0.4, 0.5, 0.6),
  mean_GRF_over_distance = seq(0, 3000)
)

params$take_off_velocity <- get_tov(
  mass = params$mass,
  gravity_const = params$gravity_const,
  mean_force = params$mean_GRF_over_distance,
  push_off_distance = params$push_off_distance
)

gg <- ggplot(
  params,
  aes(x = mean_GRF_over_distance, y = take_off_velocity, group = push_off_distance, color = push_off_distance)
) +
  theme_cowplot(8) +
  geom_line(color = "black") +
  scale_color_gradient(low = "black", high = "black")

direct.label(gg, list("last.bumpup", cex = 0.6))

As can be seen from this graph, higher the \(\bar{F}\), the higher the \(v_{to}\) reached. The same thing with \(h_{po}\). If the \(\bar{F}\) is below \(mg\) (in this case 735.75N), then there will be no jump at all, since force acting on the object is lower than object’s weight. We can also do this for the jump height, rather than \(v_{to}\):

params$height <- vjsim::get_height(params$take_off_velocity, params$gravity_const)

gg <- ggplot(params, aes(x = mean_GRF_over_distance, y = height, group = push_off_distance, color = push_off_distance)) +
  theme_cowplot(8) +
  geom_line(color = "black") +
  scale_color_gradient(low = "black", high = "black")

direct.label(gg, list("last.bumpup", cex = 0.6))

Long story short, assuming \(h_{po}\) is fixed, the more force one applies (particularly relative force, i.e. \(\frac{\bar{F}}{m}\)), the higher the jump. Rocket science!

Adding Force Generator ‘constraints’

One of the premise of Samozino’s model is that the previous mechanical constraints (\(v_{to}\) increasing with \(\bar{F}\)) can be combined with FV relationship of the Force Generator (i.e. \(F_0\) and \(V_0\)). As previously explained, Samozino’s model assumes that \(\bar{v}\) is equal to \(\frac{v_{to}}{2}\). Thus we get the following figure (using 75kg bodyweight and \(h_{po}\) of 0.4m:

params$mean_velocity_as_TOV_half <- params$take_off_velocity / 2

ggplot(
  filter(plot_data, profile == "Practical"),
  aes(x = x, y = y, group = group, color = profile)
) +
  theme_cowplot(8) +
  geom_line(alpha = 0.8) +
  labs(x = "Mean GRF over distance (N)", y = "Mean Velocity (m/s)") +
  geom_line(
    data = filter(params, push_off_distance == 0.4),
    aes(x = mean_GRF_over_distance, y = mean_velocity_as_TOV_half), group = 1, color = "black", alpha = 0.8
  )

The predicted \(\hat{v_{to}}\) (referenced papers refer to this as maximal \(v_{to}\), but I do not like the name, so I just refer to it as take-off velocity) is the point where two profiles intersect. This point can mathematically be expressed as:

\[ \hat{v_{to}} = \frac{h_{po}\left(\sqrt{\left(\frac{F_0^2}{4 V_0^2}-\frac{2 m (g m-F_0)}{h_{po}}\right)}-\frac{F_0}{2 V_0}\right)}{m} \]

This equation is implemented in vjsim::get_samozino_take_off_velocity function, that we will shortly use.

Look at \(\hat{v_{to}}\) this way - given the FV characteristic of the individual (as well as his/her bodyweight, \(h_{po}\) and of course gravity constant) we can predict take-off velocity and height. Thus, this is a model prediction.

As already explained, in vjsim we can implement practical Samozino profile (approximating mean velocity with half of TOV) and theoretical profile (using true mean velocity over \(h_{po}\)). We can plot these and see that the model \(v_{to}\) predictions differ, but they are both quite close to the first data point (i.e. bodyweight jump performance):

ggplot(plot_data, aes(x = x, y = y, group = group, color = profile)) +
  theme_cowplot(8) +
  geom_line(alpha = 0.8) +
  geom_point(data = plot_data_points, alpha = 0.8) +
  labs(x = "Mean GRF over distance (N)", y = "Mean Velocity (m/s)") +
  geom_line(
    data = filter(params, push_off_distance == 0.4),
    aes(x = mean_GRF_over_distance, y = mean_velocity_as_TOV_half), group = 1, color = "black", alpha = 0.8
  )

The differences between true jump height or true take-off velocity and model predicted can be traced back to the ratio between take-off velocity and mean velocity (since Samozino’s model assumes uniform acceleration and 2:1 ratio between these two). This is something we will explore in the next vignette.

This model might seems a bit circular, and that is one of my arguments against it. We use load and jump height performance to infer \(F_0\) and \(V_0\), from which we predict jump height. Besides, we estimate Force Generator constraints or determinants of performance (in this case jump height) from the performance itself. Then we claim that \(F_0\) and \(V_0\) determine jump performance (which we estimated from performance itself in the first place)?!

As seen in this and previous vignette, these parameters doesn’t recover true Force Generator characteristic. Here I highly agree with the comment by Cleather (2019) that we need to refer to this as “Load-jump height profile”. I will come back to further critique in the next vignette, so let’s continue exploring this model.

Optimization

“What should I be working on to increase my jump height?” is the question we as coaches often get. From the ballistic equation, the answer is simple - increase the force applied to the ground (or force impulse). But how should this be done? Should athlete focus on high-loads and low-velocity activities or low-load high-velocity activities?

Samozino’s model approach answering this question by making few assumptions. First, model assumes that \(F_0\) and \(V_0\) represent true Force Generator profile (mechanical) characteristics (but as seen they represent manifested and extrapolated characteristics given \(h_{po}\) and bodyweight, rather than ontological truth). Second, they make assumption that max power \(P_{max}\), which is equal to \(\frac{F_0 \times V_0}{4}\) represent some ontological characteristic of the Force Generator (again, rather than manifested performance metric for a given conditions). Third, model assumes that there are “force” and “velocity” capabilities of the entire lower extremities, and thus that there are training methods to hit those qualities. For example, one might improve force characteristics or velocity characteristics of the Force Generator.

Philosophically, these represent realist positions, where these metrics (i.e., \(F_0\), \(V_0\), and \(P_{max}\)) are treated as real mechanical characteristics of the Force Generator and thus represent determinants of performance that can be directly affected with specific training interventions. These implications will be explained and critiqued thorough this and next vignette. For now, let’s maintain this assumption.

According to Samozino et al., \(P_{max}\) is the determinant mechanical characteristic of the Force Generator. Additional, and assumed to be independent determinant, is the slope of the FV curve \(S_{FV}\):

\[ S_{FV} = -\frac{F_0}{V_0} \]

Thus, one can either use \(F_0\) and \(V_0\), or \(P_{max}\) and \(S_{FV}\) to represent mechanical characteristics of the Force Generator.

Similarly to the probing we did in the previous vignettes (e.g., how does increasing max force, max velocity or time to max activation affects the jump performance or the profile, while keeping everything else the same), we can ask ourselves the following question: “Assuming the \(P_{max}\) stays the same, but we change the \(S_{FV}\), can we increase the predicted jump height or take-off velocity?”

What does this mean? It is easier to understand it graphically. Let’s use the hypothetical profile where \(F_0\) is equal to 3,000N and \(V_0\) is equal to 2m/s. Push-off distance \(h_{po}\) is equal to 0.4m and bodyweight is equal to 75kg (gravity constants is of course, \(9.81ms^-2\)).

plot_data <- data.frame(
  x = c(3000, 0),
  y = c(0, 2)
)

gg <- ggplot(plot_data, aes(x = x, y = y)) +
  theme_cowplot(8) +
  geom_line(color = "blue", alpha = 0.8) +
  geom_line(
    data = filter(params, push_off_distance == 0.4),
    aes(x = mean_GRF_over_distance, y = mean_velocity_as_TOV_half),
    alpha = 0.8
  ) +
  labs(x = "Mean GRF over distance (N)", y = "Mean Velocity (m/s)")

gg

Max power and slope of this profile are equal to:

vjsim::get_max_power(3000, 2)
#> [1] 1500
vjsim::get_slope(3000, 2)
#> [1] -1500

Predicted \(v_{to}\) is equal to:

vjsim::get_samozino_take_off_velocity(
  F0 = 3000, V0 = 2,
  bodyweight = 75,
  push_off_distance = 0.4,
  gravity_const = 9.81
)
#> [1] 2.336561

Here comes the probing question: “If we assume \(P_{max}\) stays the same, can we tilt the FV profile so that take-off velocity is maximized (for a given \(P_{max}\))?”.

But how do we tilt the FV curve without affecting \(P_{max}\)? If you look at the equation for \(P_{max}\), which is \(\frac{F_0 \times V_0}{4}\), we tilt the slope by multiplying \(F_0\) by a specific factor \(k\), while simultaneously dividing \(V_0\) by the same factor. This way \(P_{max}\) stays the same, while we tilt the FV curve.

Here is the code and graphical representation using \(k\) of 0.8, 0.9, 1, 1.1, 1.2:

k <- c(0.8, 0.9, 1, 1.1, 1.2)

plot_data <- tibble(
  k = k,
  F0 = 3000 * k,
  V0 = 2 / k,
  max_power = vjsim::get_max_power(F0, V0),
  slope = vjsim::get_slope(F0, V0),
  predicted_TOV = vjsim::get_samozino_take_off_velocity(F0, V0, 75, 0.4, 9.81),
  group = 1,
  x = F0,
  y = 0,
)

plot_data_2 <- plot_data
plot_data_2$group <- 2
plot_data_2$x <- 0
plot_data_2$y <- plot_data_2$V0

plot_data <- rbind(
  plot_data,
  plot_data_2
)

plot_data <- plot_data %>%
  mutate(
    k = factor(k)
  )

gg <- ggplot(plot_data, aes(x = x, y = y, group = k, color = slope)) +
  theme_cowplot(8) +
  geom_line(alpha = 0.8, color = "blue") +
  geom_line(
    data = filter(params, push_off_distance == 0.4),
    aes(x = mean_GRF_over_distance, y = mean_velocity_as_TOV_half),
    alpha = 0.8, group = 1, color = "black"
  ) +
  labs(x = "Mean GRF over distance (N)", y = "Mean Velocity (m/s)") +
  scale_color_gradient(low = "black", high = "black")

direct.label(gg, list("last.points", cex = 0.5))

Each profile line on the previous figure has the same \(P_{max}\), but different slope. The original slope is equal to -1500, but as can be seen from the figure, according to this model, take-off velocity would be maximized when slope is -960. According to this model, this implies improving “velocity” characteristics of the Force Generator (whatever that might be). In coaches language, this individual should improve his “velocity” characteristics to improve take-off velocity and jump height.

We can represent these FV profiles with different slopes (but same \(P_{max}\)) with this figure:

plot_data <- tibble(
  k = seq(0.5, 1.5, length.out = 100),
  F0 = 3000 * k,
  V0 = 2 / k,
  slope = vjsim::get_slope(F0, V0),
  predicted_TOV = vjsim::get_samozino_take_off_velocity(F0, V0, 75, 0.4, 9.81),
  predicted_height = vjsim::get_height(predicted_TOV)
)

ggplot(plot_data, aes(x = slope, y = predicted_TOV)) +
  theme_cowplot(8) +
  geom_line()

There is algebraic solution to this problem of finding \(S_{fv}{opt}\) (that maximizes \(v_{to}\)), but in vjsim::get_samozino_optimal_profile function I have used the optim function to find the slope that maximizes the take-off velocity while maintaining the \(P_{max}\):

vjsim::get_samozino_optimal_profile(
  F0 = 3000,
  V0 = 2,
  bodyweight = 75,
  push_off_distance = 0.4,
  gravity_const = 9.81
)
#> $F0
#> [1] 3000
#> 
#> $F0_rel
#> [1] 40
#> 
#> $V0
#> [1] 2
#> 
#> $Pmax
#> [1] 1500
#> 
#> $Pmax_rel
#> [1] 20
#> 
#> $Sfv
#> [1] -1500
#> 
#> $Sfv_rel
#> [1] -20
#> 
#> $take_off_velocity
#> [1] 2.336561
#> 
#> $height
#> [1] 0.2782628
#> 
#> $optimal_F0
#> [1] 2527.841
#> 
#> $optimal_F0_rel
#> [1] 33.70455
#> 
#> $optimal_V0
#> [1] 2.373567
#> 
#> $optimal_height
#> [1] 0.2871468
#> 
#> $optimal_height_diff
#> [1] 0.008884008
#> 
#> $optimal_height_ratio
#> [1] 1.031927
#> 
#> $optimal_Pmax
#> [1] 1500
#> 
#> $optimal_Pmax_rel
#> [1] 20
#> 
#> $optimal_take_off_velocity
#> [1] 2.373567
#> 
#> $optimal_take_off_velocity_diff
#> [1] 0.03700627
#> 
#> $optimal_take_off_velocity_ratio
#> [1] 1.015838
#> 
#> $optimal_Sfv
#> [1] -1064.997
#> 
#> $optimal_Sfv_rel
#> [1] -14.19996
#> 
#> $Sfv_perc
#> [1] 140.8455
#> 
#> $FV_imbalance
#> [1] 40.84549
#> 
#> $probe_IMB
#> [1] 143.3826

vjsim::get_samozino_optimal_profile spits out multiple summaries, both original and optimized. To compare \(S_{fv}\) and \(S_{fv}^{optimal}\), Samozino suggested used two metrics, \(S_{FV}\%\) and \(FV_{IMB}\):

\[ \begin{aligned} S_{FV}\% &= 100 \times \frac{S_{FV}}{S_{fv}{opt}} \\ FV_{IMB} &= 100\times |1 - \frac{S_{FV}}{S_{fv}{opt}}| \end{aligned} \]

According to Samozino et al., \(S_{FV}\%\) over 100% represents FV profile too much oriented towards “force”, which implies “velocity” deficit (and the individual should spend more time fixing “velocity” deficit by performing low-load high-speed activities). When \(S_{FV}\%\) is under 100%, it represents FV profile that is too much oriented towards “velocity”, which implies “force” deficit (and the individual should spend more time fixing “force” deficit by performing high-load low-speed activities).

\(FV_{IMB}\) represents an absolute measure of this imbalance without a direction, where 0 represents optimal profile.

These summaries are automatically computed using vjsim::get_all_samozino_profiles for both practical and theoretical profiles.

samozino_profiles <- vjsim::get_all_samozino_profiles(profile_data)


datatable(samozino_profiles$data_frame, rownames = FALSE) %>%
  formatRound(columns = 1:ncol(samozino_profiles$data_frame), digits = 2)

Rather than maintaining \(P_{max}\) we can probe this by alternating increasing \(F_0\) or \(V_0\) by 10% to check which one provides more increasing take-off velocity:

tov <- vjsim::get_samozino_take_off_velocity(
  F0 = 3000,
  V0 = 2,
  bodyweight = 75,
  push_off_distance = 0.4,
  gravity_const = 9.81
)
tov
#> [1] 2.336561

tov_velocity <- vjsim::get_samozino_take_off_velocity(
  F0 = 3000,
  V0 = 2 * 1.1,
  bodyweight = 75,
  push_off_distance = 0.4,
  gravity_const = 9.81
)
tov_velocity
#> [1] 2.477157

tov_force <- vjsim::get_samozino_take_off_velocity(
  F0 = 3000 * 1.1,
  V0 = 2,
  bodyweight = 75,
  push_off_distance = 0.4,
  gravity_const = 9.81
)
tov_force
#> [1] 2.434618

If \(V_0\) increases for 10% (to 2.2m/s), take-off velocity will increase for 6.02%. If \(F0\) increases for 10% (to 3300N), take-off velocity will increase for 4.2%. Thus, increasing “velocity” seems better investment (assuming same effort to improve parameters for 10%), which is same conclusion with the previous optimal slope approach. We can use another metric using these two probing values:

\[ probe_{IMB} = \frac{v_{to}^{10\%velocity} - v_{to}}{v_{to}^{10\%force} - v_{to}} \times 100 \]

Using the above values, \(probe_{IMB}\) is equal to 143.38%, which is similar to \(S_{FV}\%\), just created different way.

We can also plot this probe (use the code below or vjsim::probe_samozino_take_off_velocity function):

probe_data <- vjsim::get_probing_data(
  args_list = list(F0 = 3000, V0 = 2),
  probe_func = function(...) list(take_off_velocity = vjsim::get_samozino_take_off_velocity(...)),
  change_ratio = seq(0.5, 1.5, length.out = 100),
  bodyweight = 75,
  push_off_distance = 0.4,
  gravity_const = 9.81
)


gg <- ggplot(
  probe_data,
  aes(x = change_ratio, y = take_off_velocity, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  xlab("Normalized parameter change") +
  ylab("Take-off velocity") +
  scale_color_manual(values = c(
    "F0" = "#4D4D4D",
    "V0" = "#5DA5DA"
  )) +
  xlim(c(0.5, 1.6)) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "grey") +
  geom_hline(yintercept = tov, linetype = "dashed", color = "grey")


direct.label(gg, list("last.bumpup", cex = 0.6))

As can be seen from the above probing experiment, increasing \(V_0\) in this profile creates higher improvements in take-off velocity, reaching the same conclusion as with optimal slope concept.

Bootstrapping optimization estimators

Like we did in the previous vignette on profiling, we can use bootstrap and bmbstats1 package to get estimates of uncertainty around optimization metrics (i.e. \(S_{FV}\%\), \(FV_{IMB}\), and \(probe_{IMB}\)).

bmbstats needs to be installed from the GitHub repository:

# Install bmbstats if you haven't already by running the following commands
# install.packages("devtools")
# devtools::install_github("mladenjovanovic/bmbstats")

library(bmbstats)

The following code (actually a wrapper code), puts optimization estimators inside the bootstrap loop (1000 resamples) and provides 90% confidence intervals (CIs) estimated using the percentile method.

boot_optimization <- function(profile_data) {
  optimization_estimators <- function(data, SESOI_lower, SESOI_upper, na.rm, init_boot) {

    # Get Samozino profile
    samozino_profile <- vjsim::get_all_samozino_profiles(
      profile_data = data
    )

    # Return chosen optimization estimators

    return(list(
      practical.Sfv_perc = samozino_profile$list$samozino_practical_profile$Sfv_perc,
      practical.FV_imbalance = samozino_profile$list$samozino_practical_profile$FV_imbalance,
      practical.probe_IMB = samozino_profile$list$samozino_practical_profile$probe_IMB,

      theoretical.Sfv_perc = samozino_profile$list$samozino_theoretical_profile$Sfv_perc,
      theoretical.FV_imbalance = samozino_profile$list$samozino_theoretical_profile$FV_imbalance,
      theoretical.probe_IMB = samozino_profile$list$samozino_theoretical_profile$probe_IMB
    ))
  }

  # Perform bootstrap
  boot_data <- bmbstats::bmbstats(
    data = profile_data,
    estimator_function = optimization_estimators,
    control = model_control(
      boot_samples = 1000,
      boot_type = "perc",
      seed = 1667
    )
  )

  # Add plot
  boot_data$graphs <- plot(boot_data)

  return(boot_data)
}

boot_data <- boot_optimization(profile_data)

boot_data$estimators
#>                  estimator       value       lower      upper
#> 1       practical.Sfv_perc 100.8027997  93.6468281 106.347263
#> 2   practical.FV_imbalance   0.8027997   0.1381622   8.225631
#> 3      practical.probe_IMB 101.2062250  93.6949472 107.002227
#> 4     theoretical.Sfv_perc 152.2394581 127.3506639 175.046929
#> 5 theoretical.FV_imbalance  52.2394581  27.3506639  75.046929
#> 6    theoretical.probe_IMB 155.0620464 128.8847947 178.903669

From the 90% CIs above we can see some randomness involved in estimating optimization estimators. Using plot we can plot the bootstrap distribution (I have created this object inside the above function):

boot_data$graphs

One of the main take-away points is the difference between theoretical and practical profiles (and their optimization estimators), as well as uncertainty involved when estimating them with 6 observations. And this is all without any measurement error involved.

Probing Force Generator characteristics

What happens to Samozino profile estimators when we probe the Force Generator? Let’s find out, but let’s use more force dominant profile:

probe_data <- vjsim::probe_profile(
  mass = 75,
  max_force = 3500,
  max_velocity = 3,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80, 100),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = vjsim::get_all_samozino_profiles
)

datatable(probe_data, rownames = FALSE) %>%
  formatRound(columns = 1:ncol(probe_data), digits = 2)

We can also plot the result of probing. I will split them in two figures of the most interesting estimators.

# Convert to long
probe_data_1 <- gather(probe_data, key = "variable", value = "value", -(1:8)) %>%
  filter(variable %in% c(
    "samozino_practical_profile.F0",
    "samozino_practical_profile.V0",
    "samozino_practical_profile.optimal_F0",
    "samozino_practical_profile.optimal_V0",
    "samozino_practical_profile.Pmax",

    "samozino_theoretical_profile.F0",
    "samozino_theoretical_profile.V0",
    "samozino_theoretical_profile.optimal_F0",
    "samozino_theoretical_profile.optimal_V0",
    "samozino_theoretical_profile.Pmax"
  )) %>%
  mutate(variable = factor(variable, levels = c(
    "samozino_practical_profile.F0",
    "samozino_practical_profile.V0",
    "samozino_practical_profile.optimal_F0",
    "samozino_practical_profile.optimal_V0",
    "samozino_practical_profile.Pmax",

    "samozino_theoretical_profile.F0",
    "samozino_theoretical_profile.V0",
    "samozino_theoretical_profile.optimal_F0",
    "samozino_theoretical_profile.optimal_V0",
    "samozino_theoretical_profile.Pmax"
  )))

gg <- ggplot(
  probe_data_1,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  scale_color_manual(values = c(
    "mass" = "#4D4D4D",
    "max_force" = "#5DA5DA",
    "max_velocity" =  "#FAA43A",
    "push_off_distance" = "#60BD68",
    "time_to_max_activation" = "#B276B2"
  )) +
  xlim(c(0.8, 1.35))

direct.label(gg, list("last.bumpup", cex = 0.4))

The above figure depicts effects of probing Force Generator characteristics on selected profile estimators. Let’s see what happens to selected optimization parameters including predicted jump height:

# Convert to long
probe_data_2 <- gather(probe_data, key = "variable", value = "value", -(1:8)) %>%
  filter(variable %in% c(
    "samozino_practical_profile.Sfv_perc",
    "samozino_practical_profile.optimal_height",
    "samozino_practical_profile.height",

    "samozino_theoretical_profile.Sfv_perc",
    "samozino_theoretical_profile.optimal_height",
    "samozino_theoretical_profile.height"
  )) %>%
  mutate(variable = factor(variable, levels = c(
    "samozino_practical_profile.Sfv_perc",
    "samozino_practical_profile.optimal_height",
    "samozino_practical_profile.height",

    "samozino_theoretical_profile.Sfv_perc",
    "samozino_theoretical_profile.optimal_height",
    "samozino_theoretical_profile.height"
  )))

gg <- ggplot(
  probe_data_2,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  scale_color_manual(values = c(
    "mass" = "#4D4D4D",
    "max_force" = "#5DA5DA",
    "max_velocity" =  "#FAA43A",
    "push_off_distance" = "#60BD68",
    "time_to_max_activation" = "#B276B2"
  )) +
  xlim(c(0.8, 1.35))

direct.label(gg, list("last.bumpup", cex = 0.4))

What is interesting from the above figure, is that both practical and theoretical \(S_{FV}\%\) showed predictable behavior when increasing Force Generator max force and max velocity. In other words, increasing max force increases \(S_{FV}\%\) and decreasing max force decreases \(S_{FV}\%\). When it comes to max velocity, increasing it decreases \(S_{FV}\%\) and decreasing max velocity increases \(S_{FV}\%\).

Let’s perform this same analysis, but now using more velocity dominant profile:

probe_data_velocity <- probe_profile(
  mass = 75,
  max_force = 2500,
  max_velocity = 8,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = vjsim::get_all_samozino_profiles
)
# Convert to long
probe_data_1 <- gather(probe_data_velocity, key = "variable", value = "value", -(1:8)) %>%
  filter(variable %in% c(
    "samozino_practical_profile.F0",
    "samozino_practical_profile.V0",
    "samozino_practical_profile.optimal_F0",
    "samozino_practical_profile.optimal_V0",
    "samozino_practical_profile.Pmax",

    "samozino_theoretical_profile.F0",
    "samozino_theoretical_profile.V0",
    "samozino_theoretical_profile.optimal_F0",
    "samozino_theoretical_profile.optimal_V0",
    "samozino_theoretical_profile.Pmax"
  )) %>%
  mutate(variable = factor(variable, levels = c(
    "samozino_practical_profile.F0",
    "samozino_practical_profile.V0",
    "samozino_practical_profile.optimal_F0",
    "samozino_practical_profile.optimal_V0",
    "samozino_practical_profile.Pmax",

    "samozino_theoretical_profile.F0",
    "samozino_theoretical_profile.V0",
    "samozino_theoretical_profile.optimal_F0",
    "samozino_theoretical_profile.optimal_V0",
    "samozino_theoretical_profile.Pmax"
  )))

gg <- ggplot(
  probe_data_1,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  scale_color_manual(values = c(
    "mass" = "#4D4D4D",
    "max_force" = "#5DA5DA",
    "max_velocity" =  "#FAA43A",
    "push_off_distance" = "#60BD68",
    "time_to_max_activation" = "#B276B2"
  )) +
  xlim(c(0.8, 1.35))

direct.label(gg, list("last.bumpup", cex = 0.4))

# Convert to long
probe_data_2 <- gather(probe_data_velocity, key = "variable", value = "value", -(1:8)) %>%
  filter(variable %in% c(
    "samozino_practical_profile.Sfv_perc",
    "samozino_practical_profile.optimal_height",
    "samozino_practical_profile.height",

    "samozino_theoretical_profile.Sfv_perc",
    "samozino_theoretical_profile.optimal_height",
    "samozino_theoretical_profile.height"
  )) %>%
  mutate(variable = factor(variable, levels = c(
    "samozino_practical_profile.Sfv_perc",
    "samozino_practical_profile.optimal_height",
    "samozino_practical_profile.height",

    "samozino_theoretical_profile.Sfv_perc",
    "samozino_theoretical_profile.optimal_height",
    "samozino_theoretical_profile.height"
  )))

gg <- ggplot(
  probe_data_2,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  scale_color_manual(values = c(
    "mass" = "#4D4D4D",
    "max_force" = "#5DA5DA",
    "max_velocity" =  "#FAA43A",
    "push_off_distance" = "#60BD68",
    "time_to_max_activation" = "#B276B2"
  )) +
  xlim(c(0.8, 1.35))

direct.label(gg, list("last.bumpup", cex = 0.4))

Both theoretical and practical profiles (using predicted height) showed higher improvements in predicted jump height when max force improves. I will further explore these topics in the next vignette (and at the end of this one), but so far, this model seems to predict intervention effects of increasing either max force or max velocity based on \(S_{FV}\%\) metric.

Alternatives?

The goal of optimization profile is to help answer intervention question: “What should athlete do?”. This is a complex question, but for the sake of this vignette, we are answering it from the vertical jump perspective: “What should athlete do to yield biggest effect in jump performance?”. In other way, should athlete focus on high-load and low-velocity activities (e.g. heavy squatting, or deadlifting), or low-load and high-velocity activities (e.g. unloaded jumping or even assisted jumping). Samozino optimization profile seems to help in answering this question. But let’s explore another approach and see if it works.

In the previous vignette I have demonstrated multitude of profiles that we can make. The simplest profile is Load~TOV, or mass and take-off velocity. This is actually THE profile we get from jumping, from which the Samozino profile is further developed using assumptions of constants acceleration and equal mean force over distance and time. Figure below depicts Load~TOV profile, termed simple profile which is estimated using vjsim::get_simple_profile function.

simple_profile <- vjsim::get_simple_profile(profile_data)

plot_data <- data.frame(
  x = c(simple_profile$L0, 0),
  y = c(0, simple_profile$TOV0)
)

plot_data_points <- profile_data %>%
  mutate(x = mass, y = take_off_velocity)

ggplot(plot_data, aes(x = x, y = y)) +
  theme_cowplot(8) +
  geom_vline(xintercept = profile_data$bodyweight[1], color = "grey", linetype = "dashed") +
  geom_hline(yintercept = simple_profile$take_off_velocity[1], color = "grey", linetype = "dashed") +
  geom_line(alpha = 0.8) +
  geom_point(data = plot_data_points, alpha = 0.8) +
  labs(x = "Load (kg)", y = "Take-off Velocity (m/s)")

Vertical dashed line indicates bodyweight (75kg in this case), and horizontal dashed line indicates predicted take-off velocity. As can be seen it is very close to true take-off velocity during the bodyweight jump.

Similar to Samozino model, we can check how changing slope of the profile affects the predicted vertical jump height. The logic is same as with Samozino profile, but simpler: we check where the profile line cuts the vertical dashed line:

k <- c(0.8, 0.9, 1, 1.1, 1.2)

plot_data <- tibble(
  k = k,
  L0 = simple_profile$L0 * k,
  TOV0 = simple_profile$TOV0 / k,
  slope = vjsim::get_slope(L0, TOV0),
  predicted_TOV = vjsim::get_simple_take_off_velocity(L0, TOV0, 75),
  group = 1,
  x = L0,
  y = 0,
)

plot_data_2 <- plot_data
plot_data_2$group <- 2
plot_data_2$x <- 0
plot_data_2$y <- plot_data_2$TOV0

plot_data <- rbind(
  plot_data,
  plot_data_2
)

plot_data <- plot_data %>%
  mutate(
    k = factor(k),
    slope = round(slope, 2)
  )

gg <- ggplot(plot_data, aes(x = x, y = y, group = k, color = slope)) +
  theme_cowplot(8) +
  geom_line(alpha = 0.8, color = "blue") +
  geom_vline(xintercept = profile_data$bodyweight[1], color = "grey", linetype = "dashed") +
  labs(x = "Mean GRF over distance (N)", y = "Mean Velocity (m/s)") +
  scale_color_gradient(low = "black", high = "black")

direct.label(gg, list("last.points", cex = 0.5))

The results of this search can be represented with the following figure:

plot_data <- tibble(
  k = seq(0.3, 2.5, length.out = 100),
  L0 = simple_profile$L0 * k,
  TOV0 = simple_profile$TOV0 / k,
  slope = vjsim::get_slope(L0, TOV0),
  predicted_TOV = vjsim::get_simple_take_off_velocity(L0, TOV0, 75)
)

ggplot(plot_data, aes(x = slope, y = predicted_TOV)) +
  theme_cowplot(8) +
  geom_line() +
  geom_vline(xintercept = simple_profile$Sfv, color = "grey") +
  geom_vline(xintercept = simple_profile$optimal_Sfv, color = "grey", linetype = "dashed")

Vertical dashed line represents optimal slope estimated using vjsim::get_simple_profile function. Here is the output of the same function:

simple_profile
#> $L0
#> [1] 265.1471
#> 
#> $TOV0
#> [1] 3.471573
#> 
#> $Sfv
#> [1] -76.37664
#> 
#> $Sfv_rel
#> [1] -1.018355
#> 
#> $take_off_velocity
#> [1] 2.489597
#> 
#> $height
#> [1] 0.3159069
#> 
#> $optimal_L0
#> [1] 150
#> 
#> $optimal_TOV0
#> [1] 6.136516
#> 
#> $optimal_Sfv
#> [1] -24.44384
#> 
#> $optimal_Sfv_rel
#> [1] -0.3259178
#> 
#> $optimal_take_off_velocity
#> [1] 3.068258
#> 
#> $optimal_height
#> [1] 0.479827
#> 
#> $Sfv_perc
#> [1] 312.4576
#> 
#> $FV_imbalance
#> [1] 212.4576
#> 
#> $probe_IMB
#> [1] 278.8824

Using vjsim::probe_simple_take_off_velocity we can further probe the behavior of the predicted take-off velocity when, in this case \(L_0\) and \(TOV_0\) change.

simple_probe_data <- probe_simple_take_off_velocity(
  L0 = simple_profile$L0,
  TOV0 = simple_profile$TOV0,
  bodyweight = 75,
  change_ratio = seq(0.8, 1.2, length.out = 1001)
)

ggplot(
  simple_probe_data,
  aes(
    x = change_ratio,
    y = take_off_velocity,
    color = probing
  )
) +
  theme_cowplot(8) +
  geom_line() +
  geom_vline(xintercept = 1, color = "grey", linetype = "dashed") +
  geom_hline(yintercept = simple_profile$take_off_velocity, color = "grey", linetype = "dashed")

Let’s probe this model when the Force Generator characteristics change. We can use the“force-dominant” and “velocity-dominant” profiles and compare it to Samozino model behavior

probe_data_velocity <- vjsim::probe_profile(
  mass = 75,
  max_force = 2500,
  max_velocity = 8,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = function(...) {
    list(list = vjsim::get_simple_profile(...))
  }
)


probe_data_force <- vjsim::probe_profile(
  mass = 75,
  max_force = 3500,
  max_velocity = 3,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = function(...) {
    list(list = vjsim::get_simple_profile(...))
  }
)

probe_data <- rbind(
  data.frame(profile = "Velocity dominant", probe_data_velocity),
  data.frame(profile = "Force dominant", probe_data_force)
)

probe_data <- gather(probe_data, key = "variable", value = "value", -(1:9)) %>%
  filter(variable %in% c("L0", "TOV0", "height", "optimal_Sfv", "Sfv_perc", "probe_IMB")) %>%
  mutate(variable = factor(variable, levels = c("L0", "TOV0", "height", "optimal_Sfv", "Sfv_perc", "probe_IMB")))


gg <- ggplot(
  probe_data,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(profile ~ variable, scales = "free", nrow = 2, ncol = 6) +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  xlim(c(0.8, 1.4))

direct.label(gg, list("last.bumpup", cex = 0.4))

As can be seen from the above graph, this simple model is not able to “predict” effects of increasing max force or max velocity of the force generator. If we check the height (actually predicted height), we can see that in velocity dominant profile, increasing max force parameter yields biggest effects. In the force dominant profile, increasing max velocity parameter of the Force generator yields biggest effect. This is not picked up with \(S_{FV}\%\) nor \(probe_{IMB}\).

Let’s repeat the same visual (one more time) using practical Samozino method (which also uses take-off velocity):

probe_data_velocity <- vjsim::probe_profile(
  mass = 75,
  max_force = 2500,
  max_velocity = 8,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = vjsim::get_all_samozino_profiles
)


probe_data_force <- vjsim::probe_profile(
  mass = 75,
  max_force = 3500,
  max_velocity = 3,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = vjsim::get_all_samozino_profiles
)

probe_data <- rbind(
  data.frame(profile = "Velocity dominant", probe_data_velocity),
  data.frame(profile = "Force dominant", probe_data_force)
)

probe_data <- gather(probe_data, key = "variable", value = "value", -(1:9)) %>%
  filter(variable %in% c(
    "samozino_practical_profile.F0",
    "samozino_practical_profile.V0",
    "samozino_practical_profile.height",
    "samozino_practical_profile.optimal_Sfv",
    "samozino_practical_profile.Sfv_perc",
    "samozino_practical_profile.probe_IMB"
  )) %>%
  mutate(variable = factor(
    variable,
    levels = c(
      "samozino_practical_profile.F0",
      "samozino_practical_profile.V0",
      "samozino_practical_profile.height",
      "samozino_practical_profile.optimal_Sfv",
      "samozino_practical_profile.Sfv_perc",
      "samozino_practical_profile.probe_IMB"
    ),
    labels = c("F0", "V0", "height", "optimal_Sfv", "Sfv_perc", "probe_IMB")
  ))


gg <- ggplot(
  probe_data,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(profile ~ variable, scales = "free", nrow = 2, ncol = 6) +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  xlim(c(0.8, 1.4))

direct.label(gg, list("last.bumpup", cex = 0.4))

And finally using theoretical Samozino model:

probe_data_velocity <- vjsim::probe_profile(
  mass = 75,
  max_force = 2500,
  max_velocity = 8,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = vjsim::get_all_samozino_profiles
)


probe_data_force <- vjsim::probe_profile(
  mass = 75,
  max_force = 3500,
  max_velocity = 3,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = vjsim::get_all_samozino_profiles
)

probe_data <- rbind(
  data.frame(profile = "Velocity dominant", probe_data_velocity),
  data.frame(profile = "Force dominant", probe_data_force)
)

probe_data <- gather(probe_data, key = "variable", value = "value", -(1:9)) %>%
  filter(variable %in% c(
    "samozino_theoretical_profile.F0",
    "samozino_theoretical_profile.V0",
    "samozino_theoretical_profile.height",
    "samozino_theoretical_profile.optimal_Sfv",
    "samozino_theoretical_profile.Sfv_perc",
    "samozino_theoretical_profile.probe_IMB"
  )) %>%
  mutate(variable = factor(
    variable,
    levels = c(
      "samozino_theoretical_profile.F0",
      "samozino_theoretical_profile.V0",
      "samozino_theoretical_profile.height",
      "samozino_theoretical_profile.optimal_Sfv",
      "samozino_theoretical_profile.Sfv_perc",
      "samozino_theoretical_profile.probe_IMB"
    ),
    labels = c("F0", "V0", "height", "optimal_Sfv", "Sfv_perc", "probe_IMB")
  ))


gg <- ggplot(
  probe_data,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(profile ~ variable, scales = "free", nrow = 2, ncol = 6) +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  xlim(c(0.8, 1.4))

direct.label(gg, list("last.bumpup", cex = 0.4))

Both Samozino profiles (practical and theoretical) correctly identified the effects of increasing Force Generator max force or max velocity characteristics for two different profiles.

Load~Peak Force profile?

Another profile that is commonly used in practice is the Load~Peak Force profile (LPF). This profile is useful once someone has force platform, but this also makes it less practical compared to Samozino model/profile (one just needs to calculate push-off distance and use jump height which can be estimated with a simple jump mat).

In the example below, I will provide LPF profile for velocity-dominant and force-dominant examples:

profile_data_velocity <- vjsim::vj_profile(
  external_load = c(0, 20, 40, 60, 80),

  # Simulation parameters
  mass = 75,
  push_off_distance = 0.4,
  max_force = 2500,
  max_velocity = 8
)

profile_data_force <- vjsim::vj_profile(
  external_load = c(0, 20, 40, 60, 80),

  # Simulation parameters
  mass = 75,
  push_off_distance = 0.4,
  max_force = 3500,
  max_velocity = 3
)

merged_profiles <- rbind(
  data.frame(profile = "Velocity dominant", profile_data_velocity),
  data.frame(profile = "Force dominant", profile_data_force)
)

plot_merged_profiles <- function(merged_profile_data, x_var, y_var, profile = "profile") {
  df <- data.frame(
    x_var = merged_profile_data[[x_var]],
    y_var = merged_profile_data[[y_var]],
    profile = merged_profile_data[[profile]]
  )

  gg <- ggplot(df, aes(x = x_var, y = y_var, group = profile, color = profile)) +
    theme_cowplot(8) +
    geom_line(alpha = 0.8) +
    geom_point(alpha = 0.8) +
    labs(x = x_var, y = y_var)

  return(gg)
}

plot_merged_profiles(merged_profiles, "mass", "peak_GRF")

What can be noticed from this profile is that the more “force-dominant” athlete has more steep profile, while “velocity-dominant” athlete has more horizontal profile. We can say that “velocity-dominant” athlete is able to squeeze out ability to generate maximal force at lower loads (mass), and can thus benefit from developing max force capabilities (whatever they are in ontological sense). The “force-dominant” athlete, on the other hand, is unable to manifest maximal forces at lighter loads and would benefit from training that allows him to do that.

This LPF profile can be summarized with a slope metric:

LPF_velocity <- lm(peak_GRF ~ mass, profile_data_velocity)
coef(LPF_velocity)[[2]]
#> [1] 2.674102
LPF_force <- lm(peak_GRF ~ mass, profile_data_force)
coef(LPF_force)[[2]]
#> [1] 9.277718

Let’s check how this slope metric behaves when probed in for both profiles (and compare it to Samozino’s practical \(probe_{IMB}\)):

LPF_profile <- function(profile_data) {
  LPF_profile_model <- lm(peak_GRF ~ mass, profile_data)
  samozino_profile <- get_all_samozino_profiles(profile_data)

  return(list(list = list(
    slope = coef(LPF_profile_model)[[2]],
    probe_IMB = samozino_profile$list$samozino_practical_profile$probe_IMB
  )))
}

probe_data_velocity <- vjsim::probe_profile(
  mass = 75,
  max_force = 2500,
  max_velocity = 8,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = LPF_profile
)


probe_data_force <- vjsim::probe_profile(
  mass = 75,
  max_force = 3500,
  max_velocity = 3,
  time_to_max_activation = 0.3,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = LPF_profile
)

probe_data <- rbind(
  data.frame(profile = "Velocity dominant", probe_data_velocity),
  data.frame(profile = "Force dominant", probe_data_force)
)

probe_data <- gather(probe_data, key = "variable", value = "value", -(1:9))

gg <- ggplot(
  probe_data,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(profile ~ variable, scales = "free", nrow = 2, ncol = 6) +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  xlim(c(0.8, 1.4))

direct.label(gg, list("last.bumpup", cex = 0.4))

As can be seen from the figure above, increasing max force always increase \(probe_{IMB}\), while increasing max velocity always decrease \(probe_{IMB}\). We can also conclude that \(probe_{IMB}\) is below 100 for velocity dominant profile and over 100 for force dominant profile. When it comes to LPF slope, the situation is similar to \(probe_{IMB}\).

Let’s explore this little further using two different profiles: explosive profile and slow profile:

profile_data_explosive <- vjsim::vj_profile(
  external_load = c(0, 20, 40, 60, 80),

  # Simulation parameters
  mass = 75,
  push_off_distance = 0.4,
  max_force = 3000,
  max_velocity = 4,
  time_to_max_activation = 0.1
)

profile_data_slow <- vjsim::vj_profile(
  external_load = c(0, 20, 40, 60, 80),

  # Simulation parameters
  mass = 75,
  push_off_distance = 0.4,
  max_force = 3000,
  max_velocity = 4,
  time_to_max_activation = 0.5
)

merged_profiles <- rbind(
  data.frame(profile = "Explosive profile", profile_data_explosive),
  data.frame(profile = "Slow profile", profile_data_slow)
)

plot_merged_profiles(merged_profiles, "mass", "peak_GRF")

From the figure above, the explosive profile has more horizontal slope, while slow profile has steeper slope. When it comes to mean GRF-TOV profile (the basis behind Samozino profile), we get the following:

plot_merged_profiles(merged_profiles, "mean_GRF_over_distance", "take_off_velocity")

This profile is not much affected. It seems that LPF slope reflects the time to max activation and max force characteristics of the Force Generator.

When we probe LPF slope and \(probe_{IMB}\) using these two profile, we get the following:

# --------------------------------------------
probe_data_explosive <- vjsim::probe_profile(
  mass = 75,
  max_force = 3000,
  max_velocity = 4,
  time_to_max_activation = 0.1,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = LPF_profile
)


probe_data_slow <- vjsim::probe_profile(
  mass = 75,
  max_force = 3000,
  max_velocity = 4,
  time_to_max_activation = 0.5,
  time_step = 0.001,
  external_load = c(0, 20, 40, 60, 80),
  change_ratio = seq(0.8, 1.2, length.out = 7),

  # Profile variables
  profile_func = LPF_profile
)

probe_data <- rbind(
  data.frame(profile = "Explosive profile", probe_data_explosive),
  data.frame(profile = "Slow profile", probe_data_slow)
)

probe_data <- gather(probe_data, key = "variable", value = "value", -(1:9))

gg <- ggplot(
  probe_data,
  aes(x = change_ratio, y = value, color = probing)
) +
  theme_cowplot(6) +
  geom_line() +
  facet_wrap(profile ~ variable, scales = "free", nrow = 2, ncol = 6) +
  xlab("Normalized parameter change") +
  ylab("Estimator value") +
  xlim(c(0.8, 1.4))

direct.label(gg, list("last.bumpup", cex = 0.4))

\(probe_{IMB}\) shows similar values and similar behavior in both profiles, while LPF slope shows different values in each profile, but similar behavior. But can it tell us what we need to improve to reach higher vertical jump height? This will be further discussed in the exploring vignette.

Similar idea to LPF slope concept is utilized with Bosco Index where the ratio of 2x BW jump height to BW jump height is calculated. Unfortunately, it seems that neither of these two (Bosco Index nor LPF slope) gives us a practical summary such as Samozino profile (i.e. using \(S_{FV}\%\) or \(probe_{IMB}\)) in terms what part of the “curve” one needs to improve to reach higher vertical jump height. It might be that these approaches evaluate optimality of different qualities of the Force Generator.

In the next exploring vignette I will continue with the discussion and exploration of these model by using simulated data for a large amount of different combinations of Force Generator parameters. I will follow up with discussion of intervention studies using this model (Jiménez-Reyes et al. 2017; Jiménez-Reyes, Samozino, and Morin 2019) and provide final conclusion about it’s application.

Shiny App

Now start playing with the code or with the Shiny App by click on the previous link or by running the following code:

vjsim::run_simulator()

The Shiny App will allow you much more interactive environment for exploring the vjsim.

Want to learn more?

Please continue by reading Exploring vignette:

vignette("exploring-vjsim")

References

Cleather, Daniel. 2019. “Force-Velocity Profiling Is a Misnomer: Load-Jump Height Profiling Would Be a Better Name.” Preprint. SportRxiv. https://doi.org/10.31236/osf.io/58fu3.

Jiménez-Reyes, Pedro, Pierre Samozino, Matt Brughelli, and Jean-Benoît Morin. 2017. “Effectiveness of an Individualized Training Based on Force-Velocity Profiling During Jumping.” Frontiers in Physiology 7 (January). https://doi.org/10.3389/fphys.2016.00677.

Jiménez-Reyes, Pedro, Pierre Samozino, and Jean-Benoît Morin. 2019. “Optimized Training for Jumping Performance Using the Force-Velocity Imbalance: Individual Adaptation Kinetics.” Edited by Daniel Boullosa. PLOS ONE 14 (5): e0216681. https://doi.org/10.1371/journal.pone.0216681.

Samozino, P., P. Edouard, S. Sangnier, M. Brughelli, P. Gimenez, and J.-B. Morin. 2013. “Force-Velocity Profile: Imbalance Determination and Effect on Lower Limb Ballistic Performance.” International Journal of Sports Medicine 35 (06): 505–10. https://doi.org/10.1055/s-0033-1354382.

Samozino, Pierre. 2018a. “A Simple Method for Measuring Lower Limb Force, Velocity and Power Capabilities During Jumping.” In Biomechanics of Training and Testing, edited by Jean-Benoit Morin and Pierre Samozino, 65–96. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-05633-3_4.

———. 2018b. “Optimal Force-Velocity Profile in Ballistic Push-Off: Measurement and Relationship with Performance.” In Biomechanics of Training and Testing, edited by Jean-Benoit Morin and Pierre Samozino, 97–119. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-05633-3_5.

Samozino, Pierre, Jean-Benoît Morin, Frédérique Hintzy, and Alain Belli. 2008. “A Simple Method for Measuring Force, Velocity and Power Output During Squat Jump.” Journal of Biomechanics 41 (14): 2940–5. https://doi.org/10.1016/j.jbiomech.2008.07.028.

Samozino, Pierre, Jean-Benoît Morin, Frédérique Hintzy, and Alain Belli. 2010. “Jumping Ability: A Theoretical Integrative Approach.” Journal of Theoretical Biology 264 (1): 11–18. https://doi.org/10.1016/j.jtbi.2010.01.021.

Samozino, Pierre, Enrico Rejc, Pietro Enrico Di Prampero, Alain Belli, and Jean-Benoît Morin. 2012. “Optimal Force–Velocity Profile in Ballistic Movements—Altius: Citius or Fortius?” Medicine & Science in Sports & Exercise 44 (2): 313–22. https://doi.org/10.1249/MSS.0b013e31822d757a.


  1. bmbstats package and accompanying book are currently in development, but you can access them on GitHub: https://github.com/mladenjovanovic/bmbstats and https://github.com/mladenjovanovic/bmbstats-book. The bmbstats book can be read on the GitHub pages: https://mladenjovanovic.github.io/bmbstats-book/↩︎