Scrpts and data for the manuscript "The influence of playing level on season-long development of physical fitness test performance in male youth soccer".

YouthFootballProgression Introduction In this analysis we will investigate whether there are differences in performance tests between different performance levels of youth footballers (Grassroots (GR), ProYouth (PY), Performance School (PS)) across a season. We are interested in the size of any differences and the probability of differences between the groups. The data consist preand post-season measures of 20m sprint time, squat jump, change of direction, and Yo-Yo test. There is also data on biological age and maturity offset1& we will account for these in our analysis. In this document we will only analyse 20m sprint data but analysis scripts for each test are available in as supplementary material (see here). This analysis uses the R programming language. Setup First we load libraries to be used in the analysis and set other options. library(tidyverse) Warning: package 'ggplot2' was built under R version 4.4.1 Warning: package 'tibble' was built under R version 4.4.1 Warning: package 'purrr' was built under R version 4.4.1 Warning: package 'stringr' was built under R version 4.4.1 Warning: package 'lubridate' was built under R version 4.4.1 1 -- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 -- v dplyr 1.1.4 v readr 2.1.5 v forcats 1.0.0 v stringr 1.5.2 v ggplot2 4.0.0 v tibble 3.3.0 v lubridate 1.9.4 v tidyr 1.3.1 v purrr 1.1.0 -- Conflicts ------------------------------------------ tidyverse_conflicts() -- x dplyr::filter() masks stats::filter() x dplyr::lag() masks stats::lag() i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors library(ggplot2) theme_set(theme_bw()) library(readxl) Warning: package 'readxl' was built under R version 4.4.1 library(rstanarm) Loading required package: Rcpp Warning: package 'Rcpp' was built under R version 4.4.1 This is rstanarm version 2.32.1 - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors! - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults. - For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores()) library(modelbased) Warning: package 'modelbased' was built under R version 4.4.1 library(parameters) Warning: package 'parameters' was built under R version 4.4.1 2 Attaching package: 'parameters' The following object is masked from 'package:rstanarm': compare_models library(marginaleffects) # for model interpretation Warning: package 'marginaleffects' was built under R version 4.4.1 library(ggdist) # stat_halfeye Warning: package 'ggdist' was built under R version 4.4.1 library(ggrain) Registered S3 methods overwritten by 'ggpp': method from heightDetails.titleGrob ggplot2 widthDetails.titleGrob ggplot2 library(corrplot) Warning: package 'corrplot' was built under R version 4.4.1 corrplot 0.95 loaded library(gt) Warning: package 'gt' was built under R version 4.4.1 set.seed(1234) Next we read in the data, set the Group variable as a factor and define the order of that factor. 3 # get data data_in <- read_xlsx("reduced_data_Feb25.xlsx",sheet = "Master") # set Group as factor and reorder data_in <- data_in |> mutate(Group = as_factor(Group)) |> mutate(Group = fct_relevel(Group, c("GR","PY","PS"))) Exploratory analysis The data in this file include: • subj: Subject ID • Group: Player group id (GR: Grassroots, PY: ProYouth, PS: Performance School) • Age: Player age (yrs) • Height1 & Height2: Height at start & end of season respectively (cm) • Weight1 & Weight2: Weight at at start & end of season respectively (kg) • MO2 & MO2: Maturity offset at start & end of season respectively • Twenty1 & Twenty2: 20m sprint time at start & end of season respectively (s) • SJ1 & SJ2: Squat jump performance at start & end of season respectively (cm) • COD1 & COD2: Change of direction perfomance at start & end of season respectively (s) • YY1 & YY2: Yoyo test performance at start & end of season respectively (s) Before any formal modeling we will plot the change in several measures across the season. Raincloud plots2provide useful visualisation of raw data and distributional characteristics of that data. The overall pre- & post data are shown in Figure 1 pre_post_data <- data_in |> select(subj, Group, Age:YY2) |> pivot_longer(cols = -c(subj,Group,Age), names_to = "variable",values_to = "values") # add a timepoint variable pre_post_data <- mutate(separate_wider_regex(data = pre_post_data, cols = variable, c(variable = ".+",timepoint = "\\d"))) # make plots clrs = c("chocolate","dodgerblue") pre_post_data |> ggplot(aes(timepoint, values, fill = timepoint, colour = timepoint)) + geom_rain(alpha = 0.7)+ scale_fill_manual(values = clrs, name = "",# for legend 4 labels = c("Pre","Post")) + scale_colour_manual(values = clrs, guide = "none")+ # remove legend in colour aesthetic theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + facet_wrap(~variable, scales = "free_y") YY SJ Twenty Weight COD Height MO −2 0 2 20 40 60 80 130 140 150 160 170 180 3.0 3.3 3.6 3.9 5.5 6.0 6.5 7.0 7.5 30 40 50 60 1000 2000 3000 timepoint values Pre Post Figure 1: Overall pre & post data. The data all appear reasonable & there are no extreme values. There is some skew in some variables (e.g. MO,YY) but generally the data appear approximately normally distributed. Next we create a summary table detailing the pre and post measures for each performance metric for each group. vars <- c("Group","MO1","MO2","Twenty1","Twenty2", "SJ1","SJ2","COD1","COD2","SJ1","SJ2","YY1","YY2") 5 tbl_data <- data_in |> select(all_of(vars)) summary_data <- tbl_data |> group_by(Group) |> summarise_at(vars(MO1:YY2), list(mean, sd)) # get n per group group_n <- tbl_data |> summarise(n(), .by = Group) |> arrange(Group) summary_data <- summary_data |> add_column(group_n["n()"], .after = "Group") summary_data |> gt() |> cols_merge(columns = c(MO1_fn1, MO1_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(MO2_fn1, MO2_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(Twenty1_fn1, Twenty1_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(Twenty2_fn1, Twenty2_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(COD1_fn1, COD1_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(COD2_fn1, COD2_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(SJ1_fn1, SJ1_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(SJ2_fn1, SJ2_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(YY1_fn1, YY1_fn2), pattern = "{1} ({2})")|> cols_merge(columns = c(YY2_fn1, YY2_fn2), pattern = "{1} ({2})")|> cols_label(`n()`="N", MO1_fn1 = "MO Pre", MO2_fn1 = "MO Post", Twenty1_fn1 = "20m Pre", Twenty2_fn1 = "20m Post", COD1_fn1 = "COD Pre", 6 Table 1 Group N MO Pre MO Post 20m Pre 20m Post SJ Pre SJ Post COD Pre COD Post YYIRTL1 Pre YYIRTL1 Post GR 35.00 -0.91 (1.32) -0.20 (1.37) 3.67 (0.23) 3.50 (0.24) 36.48 (6.18) 40.50 (5.35) 6.53 (0.35) 6.32 (0.24) 1,632.00 (494.25) 1,778.86 (640.88) PY 52.00 -0.52 (1.30) 0.29 (1.32) 3.45 (0.20) 3.34 (0.22) 39.23 (5.54) 41.11 (5.84) 6.33 (0.31) 6.18 (0.33) 1,752.69 (535.98) 2,271.15 (717.27) PS 88.00 -1.04 (1.10) -0.34 (1.19) 3.48 (0.18) 3.43 (0.20) 39.95 (5.36) 41.04 (5.31) 6.27 (0.32) 6.18 (0.34) 2,375.45 (758.45) 2,731.14 (704.60) COD2_fn1 = "COD Post", SJ1_fn1 = "SJ Pre", SJ2_fn1 = "SJ Post", YY1_fn1 = "YYIRTL1 Pre", YY2_fn1 = "YYIRTL1 Post")|> fmt_number(decimals = 2) It would be useful to plot the differences in performance metrics we are interested in on a group by group basis. # make group level plots clrs = c("chocolate","dodgerblue") # create negative %in% function for variable extraction `%nin%`<- negate(`%in%`) plt_data <- pre_post_data |> filter(variable %nin% c("MO","Height","Weight")) |> # reorder variable names mutate(variable = fct_relevel(variable, "Twenty","COD","SJ","YY")) |> mutate(Group = fct_relevel(Group, "GR","PY","PS")) # define facet labels tests <- c(Twenty = "20m sprint (s)",COD = "COD (s)",SJ = "SJ (cm)",YY = "YYIRT1 (m)") # plot plt_data |> ggplot(aes(timepoint, values, fill = timepoint, colour = timepoint)) + geom_rain(alpha = 0.7)+ scale_fill_manual(values = clrs, name = "",# for legend labels = c("Pre","Post")) + scale_colour_manual(values = clrs, guide = "none")+ # remove legend in colour aesthetic theme(axis.text.x=element_blank(), axis.ticks.x=element_blank()) + 7 facet_grid(variable~Group, scales = "free", labeller = labeller( variable = tests ) )+ labs(x = "Timepoint",y = NULL) # ggsave("plots/groupwise_diffs.pdf") GR PY PS 20m sprint (s) COD (s) SJ (cm) YYIRT1 (m) 3.0 3.3 3.6 3.9 5.5 6.0 6.5 7.0 7.5 30 40 50 60 1000 2000 3000 Timepoint Pre Post Figure 2: The preto post differences across the season for each performance group. This figure is the same as figure 2 in the paper. Collinearity Confounding effects of age & maturity status could arise because these variables are correlated with each other and potentially with the performance metrics i.e. collinearity3. In the code below we generate and plot a correlation matrix to examine the potential for collinearity. 8 Here we want to examine collinearity between Age, maturity status (i.e. MO1 &MO2). data_in |> select(Age, MO1, MO2:YY2) |> cor() |> corrplot(method = "number",number.cex=0.9) 1.00 0.87 0.87 −0.54 −0.57 0.46 0.52 −0.48 −0.49 0.35 0.28 0.87 1.00 0.98 −0.57 −0.63 0.43 0.49 −0.46 −0.50 0.22 0.20 0.87 0.98 1.00 −0.57 −0.64 0.42 0.53 −0.46 −0.53 0.22 0.19 −0.54 −0.57 −0.57 1.00 0.82 −0.65 −0.61 0.69 0.69 −0.44 −0.50 −0.57 −0.63 −0.64 0.82 1.00 −0.57 −0.66 0.63 0.71 −0.33 −0.44 0.46 0.43 0.42 −0.65 −0.57 1.00 0.70 −0.46 −0.42 0.35 0.35 0.52 0.49 0.53 −0.61 −0.66 0.70 1.00 −0.47 −0.60 0.35 0.33 −0.48 −0.46 −0.46 0.69 0.63 −0.46 −0.47 1.00 0.65 −0.47 −0.50 −0.49 −0.50 −0.53 0.69 0.71 −0.42 −0.60 0.65 1.00 −0.45 −0.51 0.35 0.22 0.22 −0.44 −0.33 0.35 0.35 −0.47 −0.45 1.00 0.75 0.28 0.20 0.19 −0.50 −0.44 0.35 0.33 −0.50 −0.51 0.75 1.00 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Age MO1 MO2 Twenty1 Twenty2 SJ1 SJ2 COD1 COD2 YY1 YY2 Age MO1 MO2 Twenty1 Twenty2 SJ1 SJ2 COD1 COD2 YY1 YY2 Figure 3: The correlation matrix for the preto post performance metrics, maturity offset and age for the data used in this study. From Figure 3we see that age and MO are highly correlated (as expected) and there is some high correlation between some pre & post measures (also expected). None of the performance metrics are highly correlated with age or MO though. 9 data = model_data, refresh = 0)# suppress output We can review the priors using prior_summary(). prior_summary(bayes_model) Priors for model 'bayes_model' ------ Intercept (after predictors centered) Specified prior: ~ normal(location = 3.4, scale = 0.22) Adjusted prior: ~ normal(location = 3.4, scale = 0.049) Coefficients Specified prior: ~ normal(location = [0,0,0,...], scale = [2.5,2.5,2.5,...]) Adjusted prior: ~ normal(location = [0,0,0,...], scale = [2.59,1.20,1.10,...]) Auxiliary (sigma) Specified prior: ~ exponential(rate = 1) Adjusted prior: ~ exponential(rate = 4.5) ------ See help('prior_summary.stanreg') for more details Next we examine the model summary. parameters(bayes_model) Parameter | Median | 95% CI | pd | Rhat | ESS | Prior ------------------------------------------------------------------------------------- (Intercept) | 0.67 | [ 0.25, 1.07] | 99.90% | 0.999 | 2078 | Normal (3.42 +- 0.05) Twenty1 | 0.77 | [ 0.66, 0.88] | 100% | 0.999 | 2202 | Normal (0.00 +- 2.59) GroupPY | 0.02 | [-0.03, 0.08] | 78.90% | 1.000 | 2176 | Normal (0.00 +- 1.20) GroupPS | 0.08 | [ 0.02, 0.13] | 99.80% | 0.999 | 2079 | Normal (0.00 +- 1.10) age_mat | -0.04 | [-0.06, -0.02] | 99.98% | 1.000 | 2282 | Normal (0.00 +- 0.55) 16 Uncertainty intervals (equal-tailed) computed using a MCMC distribution approximation. The coefficients are almost identical to those from the frequentist model above. As in the frequentist analysis above we can begin by plotting model effects before examining the pairwise contrasts. plot_predictions(bayes_model, condition = c("Twenty1","Group")) + theme(legend.position = "bottom")+ labs(title = "Bayesian Model predictions") Ignoring unknown labels: * linetype : "Group" 3.2 3.6 4.0 3.0 3.3 3.6 3.9 Twenty1 Twenty2 Group GR PY PS Bayesian Model predictions Figure 5: The Bayesian model based prediction for post 20m sprint time with the age/maturity variable held at its mean. 17 Unsurprisingly Figure 5is almost identical to ?@fig-freq-mod-preds. Contrasts between the groups The contrasts were are interested in are the differences in the change across the season comparing each group to the others. Again we can predict these from our model at the mean of age_mat &Twenty1. # get contrasts bayes_cmp <- avg_comparisons(bayes_model, variables = list(Group = "pairwise")) bayes_cmp Contrast Estimate 2.5 % 97.5 % PS - GR 0.0752 0.0250 0.1293 PS - PY 0.0528 0.0110 0.0949 PY - GR 0.0219 -0.0331 0.0793 Term: Group Type: response Lack of p-values in Bayesian analysis can be uncomfortable for some. Bayesian inference depends on the entire posterior distribution and not on a maximum likelihood point estimate & a sampling distribution for an assumed null value. Use of an entire distribution leads to several possible probabilistic indices of “effect”. The probability of direction (pd) is defined as the proportion of the posterior distribution that is of the posterior median sign8and is useful for assessing the direction of effect in exploratory analysis9. We can extract the MCMC samples from the Bayesian model and create plots to illustrate the probability of direction for effects in the model. # get draws cmp_draws <- bayes_cmp |> get_draws() Next we create plots for each contrast. # create plot clrs <- c("cornflowerblue","chocolate2") cmp_draws |> 18 ggplot(aes(x = draw, y = contrast, fill = after_stat(x <0))) + geom_vline(xintercept = 0, linetype = "dashed", linewidth = 1)+ # useful for appearance of stat_halfeye stat_halfeye(.width = 0.95,slab_linewidth = 0.5, slab_color = "black")+ scale_fill_manual(name = "Probability of direction: ", labels = c("Positive","Negative"), values = clrs) + labs(title = "20m Sprint",x = "Plausible Differences (s)",y = "")+ theme(legend.position = "bottom", axis.text.x = element_text(size = 18), axis.text.y = element_text(size = 18), axis.title.x = element_text(size = 20), legend.text = element_text(size = 20), legend.title = element_text(size = 20), plot.title = element_text(hjust=0.5,size = 24)) 19 PS − GR PS − PY PY − GR 0.0 0.1 Plausible Differences (s) Probability of direction: Positive Negative 20m Sprint Figure 6: Probability of direction plots for the 20m sprint across each groupwise difference. This figure is the same as figure 3 in the paper. The density plots in Figure 6are posterior distributions for the groupwise differences in the line heights in Figure 5. As per Figure 5we see that GR has a higher probability of improvement over the season than PS &PY.PY also has a higher probability of improvement over the season compared to PS. Calculate probability of direction We can use the MCMC draws for each contrast to calculate probability of direction using the estimate_contrasts() function from the modelbased package. modelbased::estimate_contrasts(bayes_model, contrast = "Group", test = "pd") 20 Marginal Contrasts Analysis Level1 | Level2 | Median | 95% CI | pd ------------------------------------------------- PY | GR | 0.02 | [-0.03, 0.08] | 78.90% PS | GR | 0.08 | [ 0.02, 0.13] | 99.80% PS | PY | 0.05 | [ 0.01, 0.09] | 99.28% Variable predicted: Twenty2 Predictors contrasted: Group Predictors averaged: Twenty1 (3.5), age_mat (-6.7e-17) These pd values reflect the visual impression of groupwise differences we see in the plot above. Specifically there is a >99% probability that GR improves 20m sprint time more than PS & a ~79% probability that GR improves 20m sprint time more than PY. There is a >99% probability that PY improves sprint time more than PS. Summary The frequentist & Bayesian ANCOVA analyses suggest that the GR group has a greater improvement in 20m sprint time across the season than the PS group. The Bayesian analysis additionally tells us there is a >99% probability for this effect. The GR group also outperforms the PY group although we are less confident in this effect (~79% probability of a larger effect in GR). References 1. Mirwald, R. L., Baxter-Jones, A. D. G., Bailey, D. A. & Beunen, G. P. An assessment of maturity from anthropometric measurements.Med Sci Sports Exerc 34, 689–694 (2002). 2. Allen, M. et al. Raincloud plots: A multi-platform tool for robust data visualization. Wellcome Open Res 4, 63 (2019). 3. Slinker, B. K. & Glantz, S. A. Multiple regression for physiological data analysis: The problem of multicollinearity.The American Journal of Physiology 249, R1–12 (1985). 4. York, R. Residualization is not the answer: Rethinking how to address multicollinearity. Social Science Research 41, 1379–1386 (2012). 5. Vickers, A. J. & Altman, D. G. Statistics notes: Analysing controlled trials with baseline and follow up measurements.BMJ (Clinical research ed.) 323, 1123–1124 (2001). 21 6. Vickers, A. J. The use of percentage change from baseline as an outcome in a controlled trial is statistically inefficient: A simulation study.BMC Medical Research Methodology 1, 6 (2001). 7. Rafi, Z. & Greenland, S. Semantic and cognitive tools to aid statistical science: Replace confidence and significance by compatibility and surprise.BMC Medical Research Methodology 20, 244 (2020). 8. Makowski, D., Ben-Shachar, M. S., Chen, S. H. A. & Lüdecke, D. Indices of Effect Existence and Significance in the Bayesian Framework.Frontiers in Psychology 10, (2019). 9. Kelter, R. How to Choose between Different Bayesian Posterior Indices for Hypothesis Testing in Practice.Multivariate Behavioral Research 58, 160–188 (2023). 22