Looking beyond the mean: quantile regression for comparative physiologists Open Access

A fundamental assumption in biological sciences is that physiological systems are an emergent property of complex interactions between behavioural, morphological and biochemical pathways. Such complexity gives rise to considerable natural variation within and across physiological systems at multiple levels of organisation, which is at the intellectual core of comparative physiology (Somero, 2000). Exploration of variation by experimental biologists has shed light on a wealth of scientific phenomena explaining how life works and has uncovered many new questions that drive fascinating research agendas today. The contributions of comparative physiologists are increasingly useful and important in understanding the impacts of global change from molecular to ecosystem scales (Somero, 2011). Although there is a growing reframing of scientific thought towards the importance of variation and outliers in biological processes (Cook et al., 2021; Kar et al., 2021), most physiological studies address hypotheses about the central tendency of datasets rather than the nature of variation within them (Mykles et al., 2010). By solely focusing on means and averages, researchers can constrain their research questions and the answers they may find.

Least squares regression (LSR), one of the most widely used frequentist statistical approaches in empirical sciences, is a key component of systematic data analysis by physiologists. LSR techniques include ordinary least squares (OLS), multiple regression and ANOVA analyses. LSR has many perks: it is straightforward, parsimonious and describes relationships between dependent and independent variables well when assumptions are satisfied (but see Halsey and Perna, 2019). However, LSR is not without its drawbacks. Physiological datasets are often characterised by outliers, unequal variation (heteroscedasticity) and non-normal distributions – qualities that sometimes limit the utility of LSR. The application of LSR also shapes the way physiologists ask questions and formulate hypotheses. LSR techniques (such as OLS and ANOVA) predict the expected value (average) of the dependent variable for a given value of the independent variable. This is because LSR measures the deviation of data points away from the mean value. However, the relationship between variables based on the expected value of the response can look different to their relationship at other parts of the distribution. The application of standard LSR techniques in the analysis of scientific data is routine and normalised. This disposition of researchers to prioritise analytical methods that focus on central tendency has been called the ‘mean focus fallacy’ (Hohl, 2009), and can shape the generation of research questions to prioritise testing changes in the mean of the response. There is nothing inherently wrong with focusing on means, but not considering alternative tools which explore responses outside the centre of the distribution limits the scope of inquiry by excluding exploration of other areas of the sample distribution.

Quantile regression (QR) is an emerging statistical methodology that goes beyond the limits of traditional LSR. QR estimates models of conditional quantile functions, or the relationship between independent variables and specified percentiles of a response (Fig. 1; Koenker, 2005; Yu et al., 2003), bypassing the need for normally distributed datasets without heteroscedasticity and outliers. QR has received attention in the field of economics (Huang et al., 2017; Koenker and Hallock, 2001) and ecology (Antúnez et al., 2023; Cade and Noon, 2003). However, QR is rarely used in comparative physiology, possibly because it is not well known or understood. QR offers researchers an opportunity to explore the nature of a response over the whole of its distribution (i.e. for a given quantile of the response). Imagine a research project concerned with understanding whether environmental temperature influenced the jumping distance of an adult amphibian. A researcher using LSR (i.e. ANOVA with cool and warm exposure treatments) that detected a statistically significant 2% reduction in mean jump distance for warm amphibians might conclude that temperature decreased jump distance. However, if a proportion of the amphibians in the warm treatment experienced a 15% increase in jump distance (i.e. in the 80th percentile), QR could describe this treatment effect which is otherwise overlooked or described qualitatively.

Here, we argue that QR analyses could be better applied by comparative physiologists to detect responses that other linear models would otherwise miss by prioritising investigation of the central tendency and focusing solely on changes in the mean. To illustrate this, we analysed a simple simulated dataset using traditional LSR and QR techniques. Then, we analysed real physiological data using the same methods. Our goal was not to provide a definitive guide on the mathematics and logic underpinning QR but to point comparative physiologists towards QR as a useful tool to search beyond the mean of physiological metrics when researching natural variation in the systems they study.

Two datasets were used for LSR and QR analyses: a simulated dataset and a real physiological dataset. The simulated dataset (n=200) was generated from the normal distribution for an arbitrary control (n=100; mean=50, s.d.=2) and treatment (n=100; mean=50, s.d.=7) group in R (v4.3.0; http://www.R-project.org/). The physiological dataset was obtained from an existing work that explored the effects of temperature on the amount of DNA damage incurred by amphibian larvae when exposed to ultraviolet radiation (UVR; see Hird et al., 2023). Only a subset of the physiological data (n=100) was used in the analyses presented here. In the current study, we focused on the concentration of DNA damage in whole tadpoles acclimated to 25°C (dependent variable) when exposed to UVR at 5 temperature (independent variable) levels (10, 15, 20, 25 and 30°C; n=20 larvae per acute temperature treatment).

All analyses were conducted in the R statistical environment (v4.3.0; http://www.R-project.org/). Regression diagnostics were analysed using the performance package (v0.10.3; Lüdecke et al., 2021). Models were two-tailed and assumed a Gaussian error structure. α was at 0.05 for all tests.

For the simulated dataset, the response was compared between the arbitrary control and treatment group in an OLS regression that modelled the response as a function. The model met most assumptions of LSR except for equality of variance. Nine quantile regressions were fitted to the same data for the 10th to 90th percentile in 10% increments using the quantreg package (v5.95; https://CRAN.R-project.org/package=quantreg). Responses for OLS and QR were modelled using the categorical group as an independent dummy variable (with two levels: control or treatment). While QR is more often used with continuous response variables, we demonstrate here that it is reasonable to use QR for tests of comparison between categorical predictors because of the mathematically identical approaches of ANOVA and OLS. In both OLS and QR, standard errors for each quantile regression were constructed by bootstrapping (10,000 iterations).

For the DNA damage dataset, cyclobutene pyrimidine dimers (CPDs) formed from exposure to UVR were modelled using OLS regression with acute exposure temperature as the independent variable. CPDs were distributed non-normally so were log transformed prior to modelling with OLS regression. Assumptions were satisfied for the model. Nine quantile regressions with bootstrapped standard errors were computed as achieved using the simulated dataset above. DNA damage responses for OLS and QR were modelled using temperature as a continuous independent variable. Slope coefficients different from 0 represent a treatment effect at the conditional quantile. All bootstrapping was achieved using the xy-pair method.

For the simulated dataset, there was no change in mean response between the two arbitrary groups using LSR (t_1,198=0.54, P=0.59). QR showed no change in the response between treatment groups for the middle quantiles (0.4–0.6), an increase in the median treatment response for the upper quantiles (0.7–1.0) and a decrease in the median treatment response for the lower quantiles (0.1–0.3; Fig. 2).

For the physiological dataset, LSR showed that for every degree increase in temperature, the average (log) level of whole-animal DNA damage (UVC-irradiated equivalent) decreased by approximately −0.04 (t_1,30=−5.48, P<0.001). Quantile regression revealed that the median CPD concentration in amphibian larvae decreased as temperature increased in the lower and middle of the distribution (0.1–0.7) but did not change in the upper quantiles (0.8–0.9; Fig. 3).

The results from the simulated data and the physiological dataset highlight two opposing possibilities of QR that tell different stories compared with the LSR outputs. For the simulated dataset, QR revealed a significant effect of treatment on the tails of the distribution of the response that was not considered when comparing the near-identical means. In the physiological dataset, QR regression showed significant decreases in DNA damage with increasing temperature for the 10th to 70th quantile, consistent with the findings of LSR that DNA damage from UVR reduces as temperature increases in amphibian larvae on average. However, the 80th to 90th percentile slopes were not significant. This suggests that temperature did not influence the level of DNA damage for larvae that incurred the most damage, because some larvae across all temperatures experienced extremely high levels of DNA damage. In other words, although OLS revealed that, on average, increased temperature decreased DNA damage, QR revealed that the highest levels of DNA damage were not associated with animals being in cooler temperatures. This finding of variation in the response itself raises interesting questions about what aspects of the biology of the high temperature, high DNA damage larvae were preventing them from undergoing a reduction in DNA damage that most of the larvae in this group experienced. Perhaps more importantly, approaching a statistical investigation with a bias toward solely analysing the mean may constrain researchers from generating interesting research questions and hypotheses that address questions and hypotheses concerned with differences around the upper and/or lower response. If there was evidence that a higher DNA damage level of 7.5 represented a lethal cutoff point, questions around whether temperature influences the capacity for maximal DNA damage following UVR exposure would be of critical importance.

For hypotheses concerning the central tendency of data, QR can be considered as a complementary or even alternative statistical tool to LSR. However, QR has important limitations that must be considered prior to its application. QR is unable to be used with binary dependent variables (Staffa et al., 2019). One of the biggest concerns with QR for experimental biology and comparative physiology is obtaining adequate statistical power to detect biologically relevant effects with low sample sizes. Low achieved statistical power due to sample size constraints is a pervasive issue across scientific disciplines and statistical methodologies (Button et al., 2013; Burgess et al., 2022). Many experimental designs in comparative physiology, for example, could not logistically hope to achieve a sample size of 100 as we used in this data simulation. QR estimates for data in the tails of the distribution, where there are generally fewer samples, create higher uncertainty than at central quantiles (Koenker and d'Orey, 1987; Staffa et al., 2019). While small sample size can also lead to underpowered LSR tests, it may disproportionately limit the analysis and interpretation of extreme quantiles in QR. When dealing with low sample sizes, it is recommended to use resampling techniques such as bootstrapping. It is possible that if the sample size is too small to bootstrap, both QR and LSR would not be adequately powered to detect significant effects. However, bootstrapped QR has shown an impressive capacity to handle the methodological difficulties of small samples (Koenker and Zhao, 1994; Cade and Richards, 2006; Tarr, 2012; Nikitina et al., 2019). Bootstrapping confidence around QR estimates also provides an important alternative to focusing on the P-value (Halsey, 2019; Halsey et al., 2015). While statistical software such as the quantreg package will provide bootstrapped estimates and standard errors for extreme percentiles (i.e. 0 and 100th percentiles), extreme care should be taken in trusting these estimates when the sample size is low.

While the examples modelled here were kept conceptually simple, QR is applicable to more complex models including non-linear forms with generalised additive models and splines (Koenker, 2005; Wei et al., 2019; Chabot et al., 2016). This is particularly useful in in vivo physiological datasets where many factors may contribute to variation in the response variable. QR has proven a useful modelling technique in studies of gene expression (Deng and Chowdhury, 2022), metabolism (Chabot et al., 2016; Sparling and Fedak, 2004), toxicity (Villain et al., 2016) and more (Baqué et al., 2015; Fleeger et al., 2010; McClain and Rex, 2001; Steluti et al., 2013; Yaniv et al., 2014). These and our own analyses illustrate critical and interrelated points around the drawbacks of focusing on the mean in physiological datasets. First, the variance of the response variable can change across an independent variable (homoscedasticity), which is not captured or explored using LSR. Second, there are important reasons that the tail ends of the response distribution can behave differently to the centre of the distribution in response to independent variables. Third, research questions of interest may include or directly relate to quantiles, rather than means of physiological data. And fourth, QR is a useful statistical method to describe relationships within physiological datasets instead of or in addition to LSR. Ultimately, the inclusion of QR in the comparative physiologist's statistical toolbox can help to ask and answer more questions about variation in the natural world.

The authors gratefully acknowledge all volunteers involved with data collection for the physiological data used in this article. We acknowledge support from The University in Queensland in providing core facilities for the research. We acknowledge the cultural significance of amphibians globally and respectfully acknowledge Yuggera Country where the amphibians involved in this research were collected. We thank two anonymous reviewers for their valued critical reading and comments on earlier drafts of the manuscript that helped to improve the clarity of the article.

This article has an associated ECR Spotlight interview with Coen Hird.