TidyX Episode 108: lists, part 2

Ellis Hughes and I finish up part 2 about lists. In this episode we explain how to do a number of statistical processes across list elements (summary statistics by group, regression models by group, correlation coefficients by group) and then wrap up by showing how to use lists to build a multi-page PDF report where each page represents the contents of one list element.

To watch our screen cast, CLICK HERE.

To access our code, CLICK HERE.

Two Group Comparison – Frequentist vs Bayes – Part 2

In Part 1 of this series we were looking at data from a fake study, which evaluated the improvement in strength scores for two groups — Group 1 was a control group that received a normal training program and Group 2, the experimental group that received a special training program, designed to improve strength. In that first part we used a traditional t-test (frequentist) approach and a Bayesian approach, where we took advantage of a normal conjugate distribution. In order to use the normal-normal conjugate, we needed to make an assumption about a known population standard deviation. By using a known standard deviation it meant that we only needed to perform Bayesian updating for the mean of the distribution, allowing us to compare between group means and their corresponding standard errors. The problem with this approach is that we might not always have a known standard deviation to apply, thus we would want to be able to estimate this along with the mean — we need to estimate both parameters jointly!

Both Part 1 and Part 2 are available in a single file Rmarkdown file on my GitHub page.

Let’s take a look at the first few rows of the data to help remind ourselves what it looked like.

Linear Regression

To estimate the joint distribution of the mean and standard deviation under the Bayesian framework we will work with a regression model where group (control or experimental) is the independent variable and the dependent variable is the change in strength score. We can do this because, recall, t-tests are really just regression underneath.

Let’s look at the output of a frequentist linear regression before trying a Bayesian regression.

fit_lm <- lm(strength_score_change ~ group, data = df)
summary(fit_lm)
confint(fit_lm)

As expected, the results that we get here are the same as those that we previously obtained from our t-test in Part 1. The coefficient for the control group (-0.411) represents the difference in the mean strength score compared to the experimental group, whose mean change in strength score is accounted for in the intercept.

Because the experimental group is now represented as the model intercept, we can instead code the model without an intercept and get a mean change in strength score for both groups. This is accomplished by adding a “0” to the right side of the equation, telling R that we don’t want a model intercept.

fit_lm2 <- lm(strength_score_change ~ 0 + group, data = df)
summary(fit_lm2)
confint(fit_lm2)

Bayesian Regression

Okay, now let’s move into the Bayesian framework. We’ll utilize the help of the brilliant {brms} package, which compiles C++ and runs the Stan language under the hood, allowing you to use the simple and friendly R syntax that you’re used to.

Let’s start with a simple Bayesian Regression Model.

library(brms)

# Set 3 cores for parallel processing (helps with speed)
fit_bayes1 <- brm(strength_score_change ~ group, 
                 data = df,
                 cores = 3,
                 seed = 7849
                 )

summary(fit_bayes1)

The output here is a little more extensive than what we are used to with the normal regression output. Let’s make some notes:

  • The control group coefficient still represents the mean difference in strength score compared to the experimental group.
  • The experimental groups mean strength score is still the intercept
  • The coefficients for the intercept and control group are the same that we obtained with the normal regression.
  • We have a new parameter at the bottom, sigma, which is a value of 0.50. This value represents the shared standard deviation between the two groups. If you recall the output of our frequentist regression model, we had a value called residual standard error, which was 0.48 (pretty similar). The one thing to add with our sigma value here is, like the model coefficients, it has its own error estimate and 95% Credible Intervals (which we do not get from the original regression output).

Before going into posterior simulation, we have to note that we only got one sigma parameter. This is basically saying that the two groups in our model are sharing a standard deviation. This is similar to running a t-test with equal variances (NOTE: the default in R’s t-test() function is “var.equal = FALSE”, which is usually a safe assumption to make). To specify a sigma value for both groups we will wrap the equation in the bf() function, which is a function for specifying {brms} formulas. In there, we will indicate different sigma values for each group to be estimated. Additionally, to get a coefficient for both groups (versus the experimental group being the intercept), we will add a “0″ to the right side of the equation, similar to what we did in our second frequentist regression model above.

group_equation <- bf(strength_score_change ~ 0 + group,
                     sigma ~ 0 + group)

fit_bayes2 <- brm(group_equation, 
                 data = df,
                 cores = 3,
                 seed = 7849
                 )

summary(fit_bayes2)

Now we have an estimate for each group (their mean change in strength score from pre to post testing) and a sigma value for each group (NOTE: To get this value to the normal scale we need to take is exponential as they are on a log scale, as indicated by the links statement at the top of the model output, sigma = log.). Additionally, we have credible intervals around the coefficients and sigmas.

exp(-0.69)
exp(-0.72)

We have not specified any priors yet, so we are just using the default priors. Before we try and specify any priors, let’s get posterior samples from our model (don’t forget to exponentiate the sigma values). We will also calculate a Cohen’s d as a measure of standardized effect.

Cohen’s d = (group_diff) / sqrt((group1_sd^2 + group2_sd^2) / 2)

bayes2_draws <- as_draws_df(fit_bayes2) %>%
  mutate(across(.cols = contains("sigma"), 
                ~exp(.x)),
         group_diff = b_groupexperimental - b_groupcontrol,
         cohens_d = group_diff / sqrt((b_sigma_groupexperimental^2 + b_sigma_groupcontrol^2)/2))

bayes2_draws %>%
  head()

Let’s make a plot of the difference in means and Cohen’s d across our 4000 posterior samples.

par(mfrow = c(1,2))
hist(bayes2_draws$group_diff,
main = "Posterior Draw of Group Differences",
xlab = "Group Differences")
abline(v = 0,
col = "red",
lwd = 3,
lty = 2)
hist(bayes2_draws$cohens_d,
main = "Posterior Draw of Cohen's d",
xlab = "Cohen's d")

Adding Priors

Okay, now let’s add some priors and repeat the process of plotting the posterior samples. We will use the same normal prior for the means that we used in Part 1, Normal(0.1, 0.3) and for the sigma value we will use a Cauchy prior, Cauchy(0, 1).

 

## fit model
fit_bayes3 <- brm(group_equation, 
                 data = df,
                 prior = c(
                       set_prior("normal(0.1, 0.3)", class = "b"),
                       set_prior("cauchy(0, 1)", class = "b", dpar = "sigma")
                 ),
                 cores = 3,
                 seed = 7849
                 )

summary(fit_bayes3)

## exponent of the sigma values
exp(-0.66)
exp(-0.70)

 

## posterior draws
bayes3_draws <- as_draws_df(fit_bayes3) %>%
  mutate(across(.cols = contains("sigma"), 
                ~exp(.x)),
         group_diff = b_groupexperimental - b_groupcontrol,
         cohens_d = group_diff / sqrt((b_sigma_groupexperimental^2 + b_sigma_groupcontrol^2)/2))

bayes3_draws %>%
  head()

 

 

## plot sample of group differences and Cohen's d
par(mfrow = c(1,2))
hist(bayes3_draws$group_diff,
     main = "Posterior Draw of Group Differences",
     xlab = "Group Differences")
abline(v = 0,
       col = "red",
       lwd = 3,
       lty = 2)
hist(bayes3_draws$cohens_d,
     main = "Posterior Draw of Cohen's d",
     xlab = "Cohen's d")

Combine all the outputs together

Combine all of the results together so we can evaluate what has happened.

 

no_prior_sim_control_mu <- mean(bayes2_draws$b_groupcontrol)
no_prior_sim_experimental_mu <- mean(bayes2_draws$b_groupexperimental)
no_prior_sim_diff_mu <- mean(bayes2_draws$group_diff)

no_prior_sim_control_sd <- sd(bayes2_draws$b_groupcontrol)
no_prior_sim_experimental_sd <- sd(bayes2_draws$b_groupexperimental)
no_prior_sim_diff_sd <- sd(bayes2_draws$group_diff)

with_prior_sim_control_mu <- mean(bayes3_draws$b_groupcontrol)
with_prior_sim_experimental_mu <- mean(bayes3_draws$b_groupexperimental)
with_prior_sim_diff_mu <- mean(bayes3_draws$group_diff)

with_prior_sim_control_sd <- sd(bayes3_draws$b_groupcontrol)
with_prior_sim_experimental_sd <- sd(bayes3_draws$b_groupexperimental)
with_prior_sim_diff_sd &<- sd(bayes3_draws$group_diff) 

data.frame(group = c("control", "experimental", "difference"), observed_avg = c(control_mu, experimental_mu, t_test_diff), posterior_sim_avg = c(posterior_mu_control, posterior_mu_experimental, mu_diff), no_prior_sim_avg = c(no_prior_sim_control_mu, no_prior_sim_experimental_mu, no_prior_sim_diff_mu), with_prior_sim_avg = c(with_prior_sim_control_mu, with_prior_sim_experimental_mu, with_prior_sim_diff_mu), observed_standard_error = c(control_sd / sqrt(control_N), experimental_sd / sqrt(experimental_N), se_diff), posterior_sim_standard_error = c(posterior_sd_control, posterior_sd_experimental, sd_diff), no_prior_sim_standard_error = c(no_prior_sim_control_sd, no_prior_sim_experimental_sd, no_prior_sim_diff_sd), with_prior_sim_standard_error = c(with_prior_sim_control_sd, with_prior_sim_experimental_sd, with_prior_sim_diff_sd) ) %>%
  mutate(across(.cols = -group,
                ~round(.x, 3))) %>%
  t() %>%
  as.data.frame() %>%
  setNames(., c("Control", "Experimental", "Difference")) %>%
  slice(-1) %>%
  mutate(models = rownames(.),
    group = c("Average", "Average", "Average", "Average", "Standard Error", "Standard Error", "Standard Error", "Standard Error")) %>%
  relocate(models, .before = Control) %>%
  group_by(group) %>%
  gt::gt()

Let’s make some notes:

  • First, observed refers to the actual observed data, posterior_sim is our normal-normal conjugate (using a known standard deviation), no_prior_sim is our Bayesian regression with default priors and with_prior_sim is our Bayesian regression with pre-specified priors.
  • In the normal-normal conjugate (posterior_sim) analysis (Part 1), both the control and experimental groups saw their mean values get pulled closer to the prior leading to a smaller between group difference than we saw in the observed data.
  • The Bayesian regression with no priors specified (no_prior_sim) resulted in a mean difference that is pretty much identical to the outcome we saw with our t-test on the observed data.
  • The Bayesian Regression with specified priors (with_prior_sim) ends up being somewhere in the middle of the observe data/Bayes Regression with no priors and the normal-normal conjugate. The means for both groups are pulled close to the prior but not as much as the normal-normal conjugate means (posterior_sim). Therefore, the mean difference between groups is higher than the posterior_sim output but not as large as the observed data (because it is influenced by our prior). Additionally, the group standard errors are more similar to the observed data with the Bayesian regression with priors than the Bayesian regression without priors and the normal-normal Bayesian analysis.

Wrapping Up

We’ve gone over a few approaches to comparing two groups using both Frequentist and Bayesian frameworks. Hopefully working through the analysis in this way provides an appreciation for both frameworks. If we have prior knowledge, which we often do, it may help to code it directly into our analysis and utilize a Bayesian approach that helps us update our present beliefs about a phenomenon or treatment effect.

Both Part 1 and Part 2 are in a single file on my GitHub page.

If you notice any errors or issues feel free to reach out via email!

Visualizing Group Changes

CJ Mayes recently posted some really nice plots for visualizing group differences to Twitter.

My personal favorite was the bottom right plot, which can be a nice way of visualizing pre and post changes in a study. I believe the original plots were done in Tableau, so I’ve gone ahead and reproduced that bottom right plot in R.

You can grab the full code, all in one piece, from my GitHub page.

Simulate Some Data

library(tidyverse)
library(gghalves)
theme_set(theme_light())

set.seed(1234)
dat <- tibble( subject = LETTERS[1:26], pre = rnorm(n = 26, mean = 10, sd = 3) ) %>%
  mutate(post = pre + rnorm(n = 26, mean = 0, sd = 2))

dat_long <- dat %>%
  pivot_longer(cols = -subject) %>%
  mutate(name = factor(name, levels = c("pre", "post"))) 

Create the plot

dat_long %>%
  ggplot(aes(x = name, y = value)) +
  geom_line(aes(group = subject),
            color = "light grey",
            size = 0.7) +
  geom_point(aes(group = subject,
                 color = name),
             alpha = 0.7,
             size = 2) +
  geom_half_violin(aes(x = name),
                   data = dat_long %&gt;% filter(name == 'pre'),
                   fill = "light grey",
                   color = "light grey",
                   side = 'l',
                   alpha = 0.7) +
  geom_half_violin(aes(x = name),
                   data = dat_long %&gt;% filter(name == 'post'),
                   fill = "palegreen",
                   color = "palegreen",
                   side = 'r',
                   alpha = 0.7) +
  scale_color_manual(values = c("pre" = "light grey", "post" = "palegreen")) +
  labs(x = "Test Time",
       y = NULL,
       title = "Changes from Pre to Post") +
  theme(axis.text = element_text(face = "bold", size = 12),
        axis.title = element_text(face = "bold", size = 15),
        plot.title = element_text(size = 18),
        legend.position = "none")

 

Visualizing Vald Force Frame Data

Recently, in a sport science chat group on Facebook someone asked for an example of how other practitioners are visualizing their Vald Force Frame data. For those that are unfamiliar, Vald is a company that originally pioneered the Nordbord, for eccentric hamstring strength testing. The Force Frame is their latest technological offering, designed to help practitioners test abduction and adduction of the hips for performance and return-to-play purposes.

The Data

I never used the Force Frame, personally, so the author provided a screen shot of what the data looks like. Using that example, I created a small simulation of data to try and create a visual that might be useful for practitioners. Briefly, the data is structured as two rows per athlete, a row representing the force output squeeze (adduction) and a row representing pull (abduction). My simulated data looks like this:

### Vald Force Frame Visual
library(tidyverse)

set.seed(678)
dat <- tibble( player = rep(c("Tom", "Alan", "Karl", "Sam"), each = 2), test = rep(c("Pull", "Squeeze"), times = 4) ) %>%
  mutate(l_force = ifelse(test == "Pull", rnorm(n = nrow(.), mean = 305, sd = 30),
                          rnorm(n = nrow(.), mean = 360, sd = 30)),
         r_force = ifelse(test == "Pull", rnorm(n = nrow(.), mean = 305, sd = 30),
                          rnorm(n = nrow(.), mean = 360, sd = 30)),
         pct_imbalance = ifelse(l_force &gt; r_force, ((l_force - r_force) / l_force) * -1,
                            (r_force - l_force) / r_force))

In this simplified data frame we see the two rows per athlete with left and right force outputs. I’ve also calculated the Bilateral Strength Asymmetry (BSA), indicated in the percent imbalance column. This measure, as well as several other measures of asymmetry, was reported by Bishop et al (2016) in their paper on Calculating Asymmetries. The equation is as follows:

BSA = (Stronger Limb – Weaker Limb) / Stronger Limb

Additionally, if the left leg was stronger than the right I multiplied the BSA by -1, so that the direction of asymmetry (the stronger limb) can be reflected in the plot.

The Plot

Prior to plotting, I set some theme elements that allow me to style the axis labels, axis text, plot title and subtitle, and the headers of the two facets of tests that we have (one for pull and one for squeeze). Additionally, in order to make the plot look cleaner, I get rid of the legend since it didn’t offer any new information that couldn’t be directly interpreted by looking at the visual.

Before creating the visual, I first add a few variables that will help me give more context to the numbers. The goal of this visual is to help practitioners look at the tests results for a larger group of athletes and quickly identify those athletes that may require specific consideration. Therefore, I create a text string that captures the left and right force outputs so that they can be plotted directly onto the plot and a corresponding “flag” variable that indicates when an athlete may be below some normalized benchmark of strength in the respective test. Finally, I created an asymmetry flag to indicate when an athlete has a BSA that exceeds 10% in either direction. This threshold can (and should) be whatever is meaningful and important with your athletes and in your sport.

For the plot itself, I decided that plotting the BSA values for both tests would be something that practitioners would find valuable and comprehending the direction of asymmetry in a plot is also very easy. Remember, the direction that the bar is pointed represents the stronger limb. To provide context for the asymmetry direction, I created a shaded normative range and whenever the bar is outside of this range, it changes to red. When it is inside the range it remains green. To provide the raw value force numbers, I add those to the plot as labels, in the middle of each plot region for each athlete. If the athlete is flagged as having a force output for either leg that is below the predetermined threshold the text will turn red.

 

theme_set(theme_classic() +
          theme(strip.background = element_rect(fill = "black"),
                strip.text = element_text(size = 13, face = "bold", color = "white"),
                axis.text = element_text(size = 13, face = "bold"),
                axis.title.x = element_text(size = 14, face = "bold"),
                plot.title = element_text(size = 18),
                plot.subtitle = element_text(size = 14),
                legend.position = "none"))

dat %>%
  mutate(left = paste("Left =", round(l_force, 1), sep = " "),
         right = paste("Right =", round(r_force, 1), sep = " "),
         l_r = paste(left, right, sep = "\n"),
         asym_flag = ifelse(abs(pct_imbalance) > 0.1, "flag", "no flag"),
         weakness_flag = ifelse((test == "Pull" & (l_force < 250 | r_force < 250)) |
                                 (test == "Squeeze" & (l_force < 330 | r_force < 330)), "flag", "no flag")) %>%
  ggplot(aes(x = pct_imbalance, y = player)) +
  geom_rect(aes(xmin = -0.1, xmax = 0.1),
            ymin = 0,
            ymax = Inf,
            fill = "light grey",
            alpha = 0.3) +
  geom_col(aes(fill = asym_flag),
           alpha = 0.6,
           color = "black") +
  geom_vline(xintercept = 0, 
             size = 1.3) +
  annotate(geom = "text",
           x = -0.2, 
           y = 4.5,
           size = 6,
           label = "Left") +
  annotate(geom = "text",
           x = 0.2, 
           y = 4.5,
           size = 6,
           label = "Right") +
  geom_label(aes(x = 0, y = player, label = l_r, color = weakness_flag)) +
  scale_color_manual(values = c("flag" = "red", "no flag" = "black")) +
  scale_fill_manual(values = c("flag" = "red", "no flag" = "palegreen")) +
  labs(x = "% Imbalance",
       y = NULL,
       title = "Force Frame Team Testing",
       subtitle = "Pull Weakness < 250 | Squeeze Weakness < 330") +
  facet_grid(~test) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 0.1),
                     limits = c(-0.3, 0.3))

 

At a quick glance we can notice that 3 of the athletes exhibit a strength asymmetry for Pull and two exhibit a strength asymmetry for Squeeze. Additionally, one of the athletes, Karl, is also below the strength threshold for both Pull and Squeeze while Sam is exhibiting strength below the threshold for Squeeze only.

Wrapping Up

There are a lot of ways to visualize this sort of single day testing data. This is just one example and would be different if we were trying to visualize serial measurements, where we are tracking changes over time. Hopefully this short example provides some ideas. If you’d like to play around with the code and adapt it to your own athletes, you can access it on my GitHub page.

TidyX Episode 107: lists, part 1

This week, Ellis Hughes and I discuss a very important and useful object in R, lists. Lists are critical to know about because of their versatility for storing any type of data formats and their ability to handle lots of data tasks very quickly.

In this first part on lists we discuss how to create a list, access the contents within a list, show their versatility, and walk through examples of how to iterate over lists with for loops and base R functions such as lapply.

To watch out screen cast, CLICK HERE.

To access our code, CLICK HERE.