Technical writeup: Reference range meta-analysis for LT1 and LT2 as %HRmax, %VO2max, and %HRR

This is an extremely technical writeup detailing the analysis behind my recent article on individual variation in heart rates at LT1 and LT2 in runners. You probably want to read that article instead.

This write-up is mostly to (1) “show my work” and (2) remind future me how I did this analysis in case I ever want to formally publish it. If you want to see some cutting-edge biostatistics in action, read on.

Studies

I did a non-systematic search for studies on runners that measured LT1 and/or LT2 using blood lactate data, and expressed LT and/or LT2 as percentages of VO2max, HRmax, or HR reserve (HRR). Studies needed to report group means and standard deviations (SD), or provide raw data, to be included.

I did not include studies that measured MLSS or critical speed, nor did I include studies that estimated LT1 or LT2 using gas exchange data (e.g. VT1, VT2, RCP).

At a later date, I might update this meta-analysis to (1) use a true systematic search on PubMed and Web of Science, and (2) include studies that estimated LT1 with GET/VT1.

For each study, I recorded the mean and standard deviation at the group level, at the smallest “level of resolution” possible. For example, Maldonado-Martin 2004 reported mean/SD values separately for men and women; I recorded these as separate groups.

When reported, I also recorded the reported age, VO2max, HRmax, and male/female distribution of the subjects. I also made a subjective classification of the subjects into recreational, trained, highly trained, and elite categories based on my “expert judgement.” If needed, I converted standard error of the mean (SEM) back to SD using the reported number of subjects.

In some cases I had to make a judgement call about which precise metric to include. Weltman 1990, for example, provided different numbers for different fixed levels of lactate. In these cases I chose the metrics that corresponded most closely to “LT1 = baseline + 0.5 mM” and “LT2 = baseline + 1.5 mM” based on Jamnick 2020.

Special cases

Denis 2021 provided a table of raw subject-level data. For consistency, I calculated the mean and SD of the subject-level data for inclusion alongside the other studies.

Fleckenstein 2023 provided a plot with raw subject-level data; I digitized the data points using a web plot digitizer and calculated the mean and SD of the subject-level data.

Weltmann 1990 did not provide a mean and SD, but did provide count-level data of the number of subjects above LT1 and LT2 at different percentages of VO2max, HRmax, and HRR. For these data, I exploited the Gaussian latent variable properties of the probit model by fitting a linear probit regression model to the percentage-over-threshold data as follows:

mod <- glm(pct_above_LT1 ~ metric_pct, family = binomial(link= "probit"), data = metric_df)

From the $\beta_0$ and $\beta_1$ coefficients of this model, you can reconstruct the mean and SD of the latent Gaussian distribution that would have generated this percentage-above-threshold data.

Theoretically it might be better to assume log-normality, as I do later on, but it is not obvious how to define a log-normal latent binomial model in any statistically meaningful way, and the normal distribution model fit the data quite well.

Disappointingly, Weltman 1990 was the only study to report LT1 or LT2 as a percentage of HRR. In the main article, I report two estimates for HRR: one using Weltman’s data only, and another using %VO2max, relying on the close HRR-VO2max correspondence described in Swain 1997 (for runners, VO2max is high enough that VO2max and VO2reserve are practically equivalent at >50% VO2max, which is the range our analysis takes place in).

Data processing

The primary outcome of interest is the individual range of variation in LT1 and LT2, expressed as a %VO2max or %HRmax. In other words, we are interested in a meta-analytic reference range. Siegel et al. 2020 and Siegel et al. 2022 provide a method for calculating just such a statistic using a REML-based random-effects meta-analysis.

The raw mean and SD values suggest some physiologically implausible results: the subjects in Santos-Concejero 2013, for example, report that, in well-trained runners, LT2 occurs at 90 +/- 7.2% VO2max. Naively interpreted, this range suggests that ~10% of runners have an LT2 that occurs above 100% VO2max, which is clearly impossible. As such, a data transformation that prevents these sorts of implausible values is necessary.

Fortunately, Siegel et al. provide an analysis method both for “raw” mean and SD values, and also a transformation for log-normally-transformed data.

Careful consideration of our data suggests that a $x^* = log(100 - x)$ transformation would have the desired properties, while still being amenable to meta-analytic determination of a reference range.

The transformation would be left-skewed and unable to exceed 100, which is exactly what we want. Additionally, adding or subtracting values to the mean before transformation does not affect the standard deviation.

The figure below shows the mean and +/- 2 SD values for each group estimate, both on the standard scale (red bars) and the back-transformed log(100-x) estimate (black bars).

The result is exactly what we want: mean estimates very close to the original raw estimates, and physiologically plausible data distributions that do not exceed 100.

The log-normal assumption is not perfect—a beta-distribution-based analysis would be better—but reference range methods do not exist yet for any other distributions! So, we will do all our analysis on the log-normal-transformed scale, then invert back for presenting our results.

Statistical analysis

From here, the analysis for each combination of (LT1, LT2) and (HRmax, VO2max) is a relatively straightforward implementation of Siegel et al’s meta-analytic procedure, with a few modifications.

First, as noted earlier, some studies report data for multiple different groups. I wanted to disentangle group-level effects (e.g. age, sex, VO2max, level of experience) as much as possible from lab-level effects (e.g. different equipment, definitions of LT1/LT2, etc.), so for analyses that included studies with these “within-study” groups, I added an additional level of random effects to the meta-analysis model. The rma.mv() function in R package “metafor” provides support for this analysis:

meta_mod <- rma.mv(yi = mean, V = vi, random = ~ 1 | study_cluster/study,  data = use_df, method = "REML")

For all analyses using this method, I checked the profile likelihood plots to ensure the model was not over-parameterized.

This multi-level analysis provides two parameters for $\tau^2$ instead of one; when necessary I assumed $\tau_1^2$ and $\tau_2^2$ were independent of one another, so the total study-to-study variance is $\tau_T^2 = \tau_1^2 + \tau_2^2$ .

Parameter interpretation

I interpreted the different variance components from the model as follows:

The pooled subject-level variance ( $\sigma^2$ ) represents the person-to-person variation in (e.g.) LT2 as a percentage of HRmax, if you knew the true, underlying biological value of LT2 with certainty.

$SE(\mu)^2$ , the squared standard error of the meta-analytic mean, represents the uncertainty in the average value of LT2 as a percentage of HRmax, across all studies. With infinitely many studies, this would eventually converge on zero.

$\tau^2$ (or $\tau_1^2 + \tau_2^2$ ) or represents the variation in how LT2 would be estimated from one lab to another, and in one group of subjects to another. It accounts for things like differences in equipment, treadmill protocol, and definition of LT2, but also differences in the subject pools from one study to the next.

In multi-level analyses, and in meta-regression analyses with moderators like VO2max or age, $\tau_1^2$ in theory represents mostly the lab-to-lab variation, with some or most of the group-to-group variation (independent of the lab) “regressed out” by the group-level random effects ( $\tau_2^2$ ) and moderators. In practice, with a small number of studies with multiple groups, and no significant moderators, it is not likely that much of the group-level variation is “absorbed” by this parameter; much of it leaks into $\tau_1^2$ .

In the full blog post, I visualize these different sources of variation as stacked-up log-normal plots, but the reference range estimate, and the final summary plot, comes from the total variance ( $\tau^2 + \sigma^2$ , or $\tau_1^2 + \tau_2^2 + \sigma^2$ ).

That’s because, for any given runner, what you want to know is essentially the following: “If I went to many different physiology labs, what range of heart rate values could my LT1/LT2 occur at?”

When you pick any random lab, you don’t know where on the distribution of lab-level variation its estimates occur, so you need to incorporate that uncertainty as well (which is exactly what Siegel et al.’s total variance captures).

Interestingly, if you do have a real estimate from a physiology lab, you should in theory still incorporate this uncertainty ( $\tau_1^2$ ) into your LT1 and LT2 heart rates, for the same reason.

Moderators were tested using the standard "mods" argument to rma() or rma.mv(). They were initially modeled as simple linear fixed effects, and plotted. There was no obvious evidence of nonlinearity (or really, any evidence of anything at all) so I did no further moderator analysis.

One thing that confused me for quite a while was why the pooled variance ( $\tau^2 + \sigma^2$ ) didn't incorporate the uncertainty in the estimated mean across studies, $\mu_{RE}$ . I eventually realized it does incorporate it: the variance in $\mu_{RE}$ is dictated purely by $tau^2$ and the various individual $\sigma_i^2$ values for each study, so pooling the individual study $\sigma_i$ values "bakes in" that contribution. So, the final reference range incorporates both the individual-level variation, and the study-to-study variation in estimated mean that's captured by $\tau^2$ .

Visualizing the distributions: Don't forget the Jacobian correction

One critical component of plotting the log-normal distribution from the meta-analytic mean and reference (total) variance is incorporating the Jacobian correction to the normal PDF values. Otherwise the density plot will be wrong.

The following code illustrates how to correctly transform the density values:

f_x <- f_x_star * (1 / (100 - x_vals))

Where f_x_star is the PDF values from the dnorm() function, when supplied with the meta-analytic mean and standard deviation (which are on the log(100-x) scale). The scaling factor comes from the Jacobian (really, just the absolute value of derivative) of the log(100-x) transformation.

Full code + data are here: Github link