It is generaly hard to do randomized trial in performance environment, especially with higher-level athletes. It is a gold standard method that tries to minimize confounders (although they do happen as I will demonstrate) between control and treatment groups, so the effects could be assigned to the treatment and not something else (confounders). This allows for certainty in making causal claims (Treatment A cause effects B).

But can we do same with observational studies and data? Less likely, since we are not controlling (by the design) all the potential confounders that might influence the effects and hence reduce the confidence in causal claims. They are easier to do though. Much easier. But can we do something about it? Reading Data Analysis Using Regression and Multilevel/Hierarchical Models, which is must have book btw, gave me couple of ideas/solutions that I want to explore in this post. So, yes this is a playbook, and I am hoping it will give you some food for thought and practical tools for data analysis.

Data set

For the purpose of this playbook, I created data set that represents hypothetical randomized trial of novel periodization method. We are interested in the effects of new periodization method to the increase in back squat 1RM.

We have 100 subjects, evenly split into control and treatment group, whose 1RM back squat is measured before and after treatment. Subject characteristics include training age, body weight and pre-test 1RM. By randomizing we want to make sure that the control and treatment groups are not-significantly different in those characteristic.

Here is the data generation script:

# Set random number generator seed so you can reproduce the scores

athletesNumber = 100

# Generate group (Control/Treatment)
group <- gl(n = 2, length = athletesNumber,
          k = 1, labels = c("Control", "Treatment"))

# Generate Age
ageTraining <- round(rnorm(n = athletesNumber, mean= 10, sd = 2), 0)

# Generate BW
bodyWeight <- round(rnorm(n = athletesNumber, mean= 80, sd = 10), 0)

# Pretest squat 1RM
preTestSquat1RM <- round(1.5 * ageTraining + 1.5 * bodyWeight +
                        rnorm(n = athletesNumber, mean= 0, sd = 10), 0)

# Intervention training volume (Arbitrary number)
trainingVolume <- round(100 + ifelse(group == "Treatment",
                                     15 + rnorm(n = athletesNumber, mean= 0, sd = 15),
                                     0 +  rnorm(n = athletesNumber, mean= 0, sd = 8)), 0)

# Posttest squat 1RM
postTestSquat1RM <- round(1.1 * preTestSquat1RM + 
                          0.15 * trainingVolume +
                          ifelse(group == "Treatment", 4, 0) +
                          rnorm(n = athletesNumber, mean= 0, sd = 8), 0)

# Difference
difference <- postTestSquat1RM - preTestSquat1RM

# Percent change
change <- round(difference / preTestSquat1RM, 2)

# Save the data in data frame
studyData <- data.frame(ageTraining,

Let’s plot the data. Each dot represent a subject (they are “jittered” for easier viewing). The vertical line and a dot represent mean and 1 SD.


# Convert to long format for plotting
studyDataLong <- melt(studyData, id.vars = "group")

gg <- ggplot(studyDataLong, aes(x = group, y = value, color = group)) +
      geom_jitter(alpha = 0.3, width = 0.4) +
      stat_summary(geom = "pointrange", = "mean_sdl",
                   size = 0.4, color = "black") +
      facet_wrap(~variable, scales = "free") +
      theme_bw() +