https://orcid.org/0000-0002-0767-8565
Abstract
The stress phenotype is multivariate. Recent advances have broadened our understanding beyond characterizing the stress response in a single dimension. Simultaneously, the toolbox available to ecophysiologists has expanded greatly in recent years, allowing the measurement of multiple biomarkers from an individual at a single point in time. Yet these advances—in our conceptual understanding and available methodologies—have not yet been combined in a unifying multivariate statistical framework. Here, we offer a brief review of the multivariate stress phenotype and describe a general statistical approach for analysis using nonparametric multivariate analysis of variance with residual randomization in permutation procedures (RRPP) implemented using the “RRPP” package in R. We also provide an example illustrating the novel insights that can be gained from a holistic multivariate approach to stress and provide a tutorial for how we analyzed these data, including annotated R code and a guide to interpretation of outputs (Online Appendix 1). We hope that this statistical methodology will provide a quantitative framework facilitating the unification of our theoretical understanding and empirical observations into a quantitative, multivariate theory of stress.
“I suppose it is tempting, if the only tool you have is a hammer, to treat everything as if it were a nail.”
—Abraham Maslow (1966), The Psychology of Science.
Introduction
Stress, as a biological concept, has proved notoriously difficult to define (Levine 1985; Koolhaas et al. 2011; Del Giudice et al. 2018). Even so, stress-related studies in wild animal populations have exploded over the past three decades (e.g., this issue of Integrative and Comparative Biology and numerous others). An emerging consensus from this research is that stress physiology is complex, with numerous traits across scales of biological organization interacting to facilitate tolerance or escape from challenging conditions (e.g., Greenberg et al. 2002; Breuner et al. 2013; Jessop et al. 2013; Del Giudice et al. 2018; MacDougall-Shackleton et al. 2019; Telemeco et al. 2019; Wada 2019). Yet a predictive understanding of the stress response and its resultant effects remains elusive (Busch and Hayward 2009; Dickens and Romero 2013; Romero et al. 2015). We propose that slow progress toward a unified, quantitative theory of biological stress is due, in part, to a mismatch between the high dimensionality of stress physiology and the univariate tools typically applied to the problem.
Historically, biological stress was often described as a single phenomenon that varies along a conceptual axis from low “stress” to high “stress.” This univariate focus is especially evident when we consider research examining the hypothalamic–pituitary–adrenal axis (HPA axis) and glucocorticoid hormones (GCs) specifically. For decades, numerous researchers quantified circulating concentrations of the most common GCs, corticosterone or cortisol (together referred to as CORT hereafter), as a primary bioindicator of the physiological stress-state of animals (elevated CORT interpreted to indicate elevated stress, for example, Romero and Wingfield 1999; Romero 2004; Hunt et al. 2006; Sheriff et al. 2011; Telemeco and Addis 2014; Gangloff et al. 2016; Vera et al. 2017). However, the role of GCs in the stress response is complicated and interpretation not straightforward (Sapolsky et al. 2000; Romero 2004; Denver 2009; Dhabhar 2009; Dantzer et al. 2014). For example, receptor type and density, transcriptional effects of bound receptors, and duration that GCs remain in circulation can vary across tissues and contexts, modulating the activity of GCs (reviewed in Breuner and Orchinik 2002; Landys et al. 2006; de Kloet et al. 2008; Groeneweg et al. 2011; Breuner et al. 2013; Desantis et al. 2013; Polman et al. 2013; Telemeco et al. 2019). Thus, measuring one endpoint of HPA activation—circulating GCs—belies the complex interactions that ultimately drive effects on behavior, metabolism, cognition, and immune function. Even the simple presumption that chronic stress results in increased HPA activation is not generally supported across taxa (Dickens and Romero 2013). Such complexities likely explain why CORT phenotype poorly correlates with survival and reproduction in many systems (Breuner et al. 2008; Bonier et al. 2009; Crespi et al. 2013). Increased recognition of such complexities has resulted in growing advocacy against use of CORT, or any single bioindicator, as a univariate measure of stress (Bradshaw 2003; Breuner et al. 2013; Schoech et al. 2013; Romero et al. 2015; Romero and Gormally 2019; Telemeco et al. 2019). While measures of GCs are certainly important (Schoech et al. 2013), it is only within the context of the full multivariate physiological phenotype that we can properly interpret the action of GCs or most other individual measures of stress physiology (Angelier and Wingfield 2013; Breuner et al. 2013; Goessling et al. 2015; Vera et al. 2017).
Even when multiple traits are measured, most analyses of stress still rely on traditional univariate tools and hypothesis tests (e.g., t-tests, analysis of variance [ANOVA], chi-square, etc.); our own prior work being no exception (e.g., Gangloff et al. 2016; Gangloff et al. 2017; Telemeco et al. 2017; Telemeco et al. 2019). Thus, stress researchers commonly measure numerous traits in multiple tissues, but then analyze those traits one-by-one and attempt to discover patterns of covariation post hoc. This approach is suboptimal for multiple reasons. First, a series of univariate analyses does not allow researchers to directly quantify covariation among measured traits, account for this covariation when performing hypothesis tests, or allow predictive descriptions of the multivariate phenotype. Additionally, multiple univariate tests require multiple-comparison corrections which correct Type-1 error at the expense of Type-2 error, and thus reduce power to detect potentially important biological consequences of stressor exposure (Perneger 1998; Cabin and Mitchell 2000; Moran 2003). Finally, drawing simple conclusions from numerous univariate tests can be exceedingly difficult, with hypothesis tests on different variables frequently contradicting one another. Rather than the covariation and context-dependency implied by such contradictions adding to our understanding of stress physiology, this situation often adds confusion. More pernicious, such complexities effectively incentivize researchers to measure fewer traits or selectively report measurements so that conclusions are “cleaner” and easier to publish.
In addition to affecting how we analyze stress-physiology data, reliance on univariate tools might influence how we think about stress. Many prominent models cast stress along a single axis of variation (Wingfield and Kitaysky 2002; McEwen and Wingfield 2003; Romero et al. 2009; McEwen and Wingfield 2010; Angelier and Wingfield 2013; Dickens and Romero 2013; Wingfield 2013; Dantzer et al. 2014). Forcing stress onto a single axis makes sense if researchers must use univariate tools to test models, but contradicts the growing consensus that stress is a complex, multivariate phenotype. While it is possible that a single axis through a physiological hypervolume defines stress (Jessop et al. 2013), multiple points or surfaces throughout this hypervolume could also represent “stressed” or “unstressed” states. This latter possibility is incongruous with the most current theory. Although recent models of stress, such as the control-theory model (Del Giudice et al. 2011; Del Giudice et al. 2018) and damage-fitness model (Wada 2019; Wada and Heidinger 2019), allow a diversity of stress responses, they largely do so by bypassing the underlying physiological mechanisms. Moreover, stress models tend to be verbal rather than mathematical (Romero et al. 2015), potentially due to difficulties converting univariate ideas into quantitative multivariate predictions. Thus, an explicit multivariate toolkit could prove immensely useful for both our analyses and theory of stress.
Recent development of multivariate statistical tools within the R statistical programming environment (R Core Team 2019) makes application of multivariate analyses to physiological datasets relatively simple. Here, we present recommendations for multivariate analyses of stress and other complex physiological traits leveraging these tools and provide an example demonstrating some of the novel insights that can be gained by using these multivariate tools. Additionally, Online Appendix 1 provides a detailed worked example and tutorial demonstrating the steps and considerations needed to analyze stress as a multivariate phenotype.
Recommendations for multivariate analyses of stress
Multivariate analyses are not new (Hotelling 1931; Wilks 1932; Fisher 1936; Bartlett 1939; Pillai 1955), but ecological physiologists tend to receive limited training in these methods. Additionally, traditional parametric multivariate tests can be limiting because they require normally distributed data (Mardia 1971; Olson 1974) and sample sizes greater than the number of response variables measured (Anderson 2001; Legendre and Legendre 2012). Parametric multivariate tests are not possible for high-dimensional datasets where many traits are measured on a few individuals, which can be a common problem for ecological physiologists and other biologists. However, recent advances in the fields of community ecology, geometric morphometrics, and evolutionary biology have removed most of these hurdles (Anderson 2001; Bolnick et al. 2018; Adams and Collyer 2019). To start, Anderson (2001) and McArdle and Anderson (2001) developed a non-parametric multivariate ANOVA (NP-MANOVA) that is robust to deviations from normality, allows a greater number of response variables than samples and is generalizable to a wide array of linear models. This NP-MANOVA relies on estimating distances between replicates within and among treatments and assesses statistical significance via full-randomization permutation tests (Anderson 2001). More recently, Collyer and others demonstrated that assessing significance of NP-MANOVA via residual randomization in permutation procedures (RRPP), rather than full randomization, improves test performance by decreasing Type-1 error rates and improving statistical power with increasing data dimensionality (reviewed in Collyer et al. 2015). Finally, Collyer and Adams (2018) introduced the R package “RRPP”, which implements the above methods in R with functions similar to those for traditional univariate analyses. In addition to models analogous to single-factor ANOVA, “RRPP” allows multivariate models analogous to linear regression, ANCOVA, multi-factor AN(C)OVA, mixed-effects AN(C)OVA, and phylogenetic generalized least squares (PGLS, Collyer and Adams 2018). More generally, the “RRPP” package in R allows both ordinary least squares (OLS) and generalized least squares estimation of model coefficients with hypothesis tests utilizing Type-1, Type-2, or Type-3 sums of squares (Collyer and Adams 2018). Thus, virtually all univariate statistical analyses commonly performed by stress researchers have a multivariate analog that can be readily implemented via NP-MANOVA with RRPP utilizing the “RRPP” package in R. With minimal additional training, stress researchers can now readily perform multivariate hypothesis tests.
The core difference between multivariate and univariate analyses is the number of dimensions in the dependent variable. Whereas univariate analyses use a single response vector (i.e., column in a spreadsheet), multivariate analyses use a response matrix (i.e., multiple columns in a spreadsheet) as the dependent variable (Legendre and Legendre 2012). Essentially, multiple univariate dependent variables are grouped side-by-side into a matrix which becomes the single dependent variable for multivariate analyses. However, care must be taken when grouping variables into a response matrix. For analyses to be valid, three conditions must be met: (1) skew and outliers should be minimized, (2) the response matrix must be complete, and (3) data must be continuous with commensurate units and scale (Legendre and Legendre 2012; Adams and Collyer 2019). As mentioned above, NP-MANOVA with RRPP does not require strict adherence to normality (Collyer et al. 2015; Collyer and Adams 2018), but because distances are used to estimate effects, strong skew and outliers can bias results and should be minimized (Legendre and Legendre 2012). Additionally, transformations necessary to ensure data are in commensurate units and scale (Condition 3) may require that data are approximately normally distributed (Legendre and Legendre 2012).
A complete response matrix implies that every variable was measured in every individual. In practice, missing data are often inevitable due to technical difficulties with an assay or individuals being removed before study completion. Missing values can be dealt with in three ways: (1) remove the entire row (i.e., delete the replicate), (2) remove the entire column (i.e., delete the trait), or (3) impute the value (i.e., fill in the value; Little and Rubin 1987; Legendre and Legendre 2012). We recommend a combination of these strategies that minimizes both data loss and over imputation—remove any individuals or traits missing a large amount of data and then impute remaining, isolated missing values. Unfortunately, we do not know of a validated rule for deciding when to delete versus impute but we suggest that ≥10% missing data are a reasonable initial cutoff. When imputing, we recommend estimating values via multiple regressions that include the variable to be imputed as the dependent variable and remaining response variables as independent variables. These multiple regression models can then be used to predict the missing value for imputation. Importantly, the independent variables for final multivariate analyses should not be included in these multiple-regression models; only the traits that will be part of the response matrix should be included to keep imputation from biasing final results. For additional discussion of imputation strategies, see Legendre and Legendre (2012).
Finally, for NP-MANOVA with RRPP models to be valid, all response variables must be represented as continuous, quantitative data in commensurate units and similar scale (Legendre and Legendre 2012; Adams and Collyer 2019). Continuous, quantitative data are necessary because NP-MANOVA with RRPP relies on Euclidean geometry (Collyer and Adams 2018). If the response matrix is composed of data better represented by a non-Euclidean distance metric (e.g., count data), an appropriate distance matrix can first be estimated and then input into the NP-MANOVA with RRPP analysis (Anderson 2001, Collyer and Adams 2018). Importantly, however, multivariate analyses are never valid if the response matrix is composed of a combination of variables with fundamentally different geometries (Legendre and Legendre 2012; Adams and Collyer 2019). For example, it is not possible to conduct a multivariate analysis with a response matrix composed of a continuous trait such as hormone concentration, a binary trait such as presence/absence of a bioindicator, and a count trait such as clutch size. Similarly, care should be taken before including measurements from multiple tissues or time points within a response matrix. Conceptually, the response matrix should describe a single biological feature, such as the shape or physiological state of an organism.
After confirming all response variables are continuous and quantitative, researchers must next ensure that variables are all in commensurate units and scale (Legendre and Legendre 2012; Adams and Collyer 2019). For some physiological data, such as gene expression or metabolite measures, the data are inherently commensurate and no further processing is necessary. This is generally true for “-omics” datasets, making multivariate analyses via NP-MANOVA with RRPP straightforward. However, other common stress-physiology datasets combine diverse measures such as circulating hormone concentrations, indices of immune function, quantity of reactive oxygen species, etc. Because these traits are not measured in the same units, they are not inherently commensurate. Moreover, even traits measured in the same units may be on drastically different scales, such as is common when multiple hormones are measured. Failure to transform such data to similar scale will result in the variable with the greatest variance having an outsized effect on conclusions. If possible, converting measures to the same units is preferable. However, when this is not possible, traits can be z-standardized by subtracting the mean from the data and then dividing by the standard deviation (SD; Legendre and Legendre 2012). The z-standardization results in all data having a mean of zero and unit variance, and thus commensurate units and scale. Because z-standardization relies on the mean and SD, data need to be approximately normally distributed before standardization (Legendre and Legendre 2012).
Although z-standardization is often necessary to perform a multivariate analysis, it is not without constraints. For example, z-standardized data all have equal weight, which might not be desirable if a trait is known to be more informative than others a priori. That said, if appropriate weights are known, they can be included in NP-MANOVA with RRPP analyses. Similarly, because all variables have similar weight by default, care should be taken to include similar numbers of variables describing features within the response matrix or appropriate weight corrections should be made to models. For example, if a response matrix includes 20 endocrine variables and 2 immune variables, the endocrine variables will by default constitute 10 times more potential variation in the data thereby implying that endocrine variables are 10 times more important for describing the organism’s physiological status. Additionally, because variances are all forced to one, z-standardized data cannot be used for analyses of dispersion (Adams and Collyer 2019).
In addition to performing hypothesis tests via NP-MANOVA with RRPP, the dominant trends within a multivariate dataset can be visualized via ordination analyses, such as a principal component analysis (PCA, Fig. 1). For a complete review of ordination analyses, see Legendre and Legendre (2012) Chapter 9. When response matrices are composed of continuous, quantitative data, we recommend PCA be used for visualization because PCA preserves relationships within the data without distortion, unlike some other methods such as canonical analyses or nonmetric multidimensional scaling (Legendre and Legendre 2012). Conceptually, PCA works by rotating the data such that the maximum variation possible is visible without distorting the relationships within the dataset. First, the analysis identifies the longest line through the high-dimensional hypervolume defined by the data. This line becomes the first principal component (PC1) and represents the axis of greatest variation within the data. Next, the analysis identifies the longest line through the hypervolume that is orthogonal (this can be thought of as perpendicular) to the first line, which becomes PC2. The process repeats until no orthogonal axes of variation remain. In practice, this process is completed via eigendecomposition of the covariance matrix or singular value decomposition (Legendre and Legendre 2012) and can be performed with the prcomp() function in base R or via numerous functions within the “RRPP” package (see Online Appendix 1). Typically, PC plots display PC1 on the x-axis and PC2 on the y-axis. Thus, the x-axis illustrates the maximum amount of variation in the dataset and the y-axis illustrates the most variation orthogonal (perpendicular) to the x-axis. When PC1 and PC2 are plotted, variation on the remaining axes is hidden from view. The remaining PCs can be examined but caution is warranted because they are largely defined by being orthogonal to PC1 and PC2. Rules including the broken-stick model and only interpreting PCs with corresponding eigenvalues greater than the mean eigenvalue of all PCs are available for deciding which of the minor PCs may be appropriate for interpretation (Legendre and Legendre 2012). Data requirements for PCA are similar to those for NP-MANOVA with RRPP (Legendre and Legendre 2012); thus, a single data preparation phase following the guidelines in the previous paragraphs will allow both hypothesis testing (via NP-MANOVA) and visualization (via PCA).
PCA plots generated using “RRPP” tools in R displaying effects of three stressor treatments on the multivariate physiological phenotype of eastern fence lizards (S. undulatus). (A) Least-squares means and 95% confidence ellipses from an NP-MANOVA with RRPP model projected onto PC space. Loadings from this PCA are provided in Table 1. (B) Fitted values from an NP-MANOVA with RRPP model for each individual projected onto PC space. Rotations are different for these two PC plots because they are derived from separate PCAs. PC1 and PC2 total 100% because three treatments can always be cast onto two dimensions. Percentages explained by the first two PCs will be lower when models include more treatment levels or factors. Code for recreating each panel is provided in Online Appendix 1.
Additional insights provided by multivariate analyses: an example with eastern fence lizards (Sceloporus undulatus)
Here, we demonstrate how a holistic multivariate analysis can add useful insights beyond those attainable with traditional univariate tools. For this example, we use data published by Telemeco et al. (2019) examining the effect of acute exposure to two natural stressors (high temperature exposure or attack by red imported fire ants [Solenopsis invicta]) on the physiology of eastern fence lizards (S. undulatus). Telemeco et al. (2019) measured multiple physiological traits across levels of biological organization, including plasma CORT and blood glucose concentrations (indicators of endocrine physiology), complement and natural antibody efficacy (CH50 and NAb, indicators of innate immunity), gene expression of a panel of five heat-shock proteins (HSPs, indicator of the cellular stress response), and qualitative escape behaviors in two life-history stages (juveniles and adult males). They used these data to test whether exposure to the two stressors induced similar or divergent stress responses and used a combination of univariate and multivariate tools, although they did not group response variables into a unified analysis as we do here. Telemeco et al. (2019) concluded that similarities between responses to the two stressors varied across levels of organization. For example, escape behaviors appeared similar, but cellular responses were drastically different and endocrine responses were mixed. For further details on experimental design, analyses, and conclusions, see Telemeco et al. (2019).
Online Appendix 1 provides a tutorial and worked example reanalyzing the data for adult male lizards from Telemeco et al. (2019) via NP-MANOVA with RRPP using the “RRPP” package in R. The tutorial walks readers through the complete multivariate analysis, including all data processing steps needed to ensure the response matrix is suitable for analysis as described in the previous section. This reanalysis does not alter the broad, qualitative conclusions of Telemeco et al. (2019), but adds unique details and insight. First, NP-MANOVA with RRPP indicated a strong effect of treatment on lizard physiology (F2,13 = 4.00, P < 0.0001), similar to the conclusion drawn by Telemeco et al. (2019). However, Telemeco et al. (2019)’s analyses suggested that this effect was driven by the high temperature treatment alone, with little difference between the fire ant treatment and controls. By contrast, post-hoc pairwise tests utilizing “RRPP” tools indicate that the multivariate physiology of animals from both stressor treatments differed from control-treatment animals (P < 0.01 for both), and the stressor treatments marginally differed from one another (P = 0.049). Moreover, Fig. 1 illustrates a PCA plot of least-squares means with 95% confidence ellipses from the NP-MANOVA with RRPP model. This analysis suggests that one axis of variation through the data generally distinguishes control individuals from those exposed to stressors (Fig. 1, PC1) while another axis describes the differences between the stressors (Fig. 1, PC2). Further examination of the PC loadings (Table 2) suggests that stressor-exposed individuals were characterized by increased CORT and HSPs, with little difference in immune function or plasma glucose. The stressor treatments, however, were primarily distinguished by a contrast between plasma glucose, CORT, and CH50 on one side, and HSPA1A (also called HSP70, the canonical HSP) and NAb on the other. Additionally, the NP-MANOVA with RRPP analysis allowed estimation of predicted means and 95% confidence intervals for each measured trait while accounting for covariation (Fig. 2). We display the scaled values in Fig. 2, but these values are readily back-transformed into their original units (see Online Appendix 1) that could be used for simulation studies or as quantitative predictions for future experiments. Examination of these predicted values adds context for interpreting the model and PCA. When taken together, the “RRPP” analyses suggest that regardless of the stressor, exposure induced elevated circulating plasma CORT and HSP-gene expression. However, the stressors also induced unique physiological changes: fire ant attack elevated plasma glucose while having little effect on CH50 or HSPA1A, whereas high-temperature exposure decreased plasma glucose and CH50 while elevating HSPA1A. These quantitative conclusions could not be drawn from the original analyses of Telemeco et al. (2019).
Least-squares means and 95% confidence intervals generated from an NP-MANOVA with RRPP model for each physiological trait included in a multivariate response matrix. Data are from eastern fence lizards (S. undulatus) exposed to three stressor treatments. Values shown are predicted values from the model after accounting for covariation within the response matrix. Values are displayed on a z-standardized scale, but can be back-transformed to original units if desired. Abbreviations for the response variables are provided in Table 1. Code for generating predicted values, plotting, and back-transforming values are provided in Online Appendix 1.
Variable loadings from a PCA on predicted values for the multivariate physiological phenotype of eastern fence lizards (S. undulatus) exposed to three acute stressor treatments
| Response variable | PC1 (63.15%) | PC2 (36.85%) |
|---|---|---|
| Glucose | −0.053 | −0.77 |
| CORT | 0.566 | −0.265 |
| CH50 | −0.241 | −0.318 |
| NAb | −0.094 | 0.24 |
| HSP90AB1 | 0.379 | −0.162 |
| HSP90B1 | 0.259 | −0.135 |
| HSPA1A | 0.494 | 0.36 |
| HSPA9 | 0.395 | −0.066 |
| Response variable | PC1 (63.15%) | PC2 (36.85%) |
|---|---|---|
| Glucose | −0.053 | −0.77 |
| CORT | 0.566 | −0.265 |
| CH50 | −0.241 | −0.318 |
| NAb | −0.094 | 0.24 |
| HSP90AB1 | 0.379 | −0.162 |
| HSP90B1 | 0.259 | −0.135 |
| HSPA1A | 0.494 | 0.36 |
| HSPA9 | 0.395 | −0.066 |
Predicted values were generated using an NP-MANOVA with RRPP model implemented using the “RRPP” package in R. Figure 1A illustrates this PCA and Table 1 provides definitions for response variable abbreviations. All variance is explained by two PCs because three treatments can always be explained by two axes (i.e., a plane through the hypervolume).
Variable loadings from a PCA on predicted values for the multivariate physiological phenotype of eastern fence lizards (S. undulatus) exposed to three acute stressor treatments
| Response variable | PC1 (63.15%) | PC2 (36.85%) |
|---|---|---|
| Glucose | −0.053 | −0.77 |
| CORT | 0.566 | −0.265 |
| CH50 | −0.241 | −0.318 |
| NAb | −0.094 | 0.24 |
| HSP90AB1 | 0.379 | −0.162 |
| HSP90B1 | 0.259 | −0.135 |
| HSPA1A | 0.494 | 0.36 |
| HSPA9 | 0.395 | −0.066 |
| Response variable | PC1 (63.15%) | PC2 (36.85%) |
|---|---|---|
| Glucose | −0.053 | −0.77 |
| CORT | 0.566 | −0.265 |
| CH50 | −0.241 | −0.318 |
| NAb | −0.094 | 0.24 |
| HSP90AB1 | 0.379 | −0.162 |
| HSP90B1 | 0.259 | −0.135 |
| HSPA1A | 0.494 | 0.36 |
| HSPA9 | 0.395 | −0.066 |
Predicted values were generated using an NP-MANOVA with RRPP model implemented using the “RRPP” package in R. Figure 1A illustrates this PCA and Table 1 provides definitions for response variable abbreviations. All variance is explained by two PCs because three treatments can always be explained by two axes (i.e., a plane through the hypervolume).
Conclusion
A growing body of evidence suggests that stress is a multivariate phenotype (e.g., Jessop et al. 2013; MacDougall-Shackleton et al. 2019; Telemeco et al. 2019; Wada and Heidinger 2019) and recent tools available through the “RRPP” package in R make multivariate analyses of stress (or any physiological phenotype) readily accessible to researchers (Collyer and Adams 2018). We further hope that our tutorial and worked example (Online Appendix 1) reduce barriers to utilizing these tools. In the short-term, we think that adoption of multivariate techniques will allow researchers to gain greater quantitative insights from their data while explicitly accounting for physiological covariation, similar to the additional insights gained from our worked example. The complexity of models accessible via NP-MANOVA with RRPP should also facilitate analyses of temporal changes in the stress response and phylogenetic comparative examinations of the stress response. Additionally, the predicted values and confidence intervals readily estimated from these multivariate models should facilitate design of future experiments and simulation studies that could advance our theoretical understanding of the stress response. In the longer-term, production of a large body of literature explicitly examining stress as a multivariate phenotype should facilitate meta-analyses searching for general patterns in the stress response across species and contexts. Identification of the key networks and regions of physiological state space associated with stress could allow us to finally move toward a predictive, quantitative theory of biological stress.
Acknowledgments
We thank Dean Adams and Michael Collyer for helpful discussions about multivariate analyses and the “RRPP” package in R.
Supplementary data
Supplementary data available at ICB online.
List of abbreviations and acronyms
| Abbreviation | Definition |
|---|---|
| AN(C)OVA | Analysis of variance or analysis of covariance |
| CH50 | Measure of plasma complement ability to lyse heterologous red blood cells independent of antibodies |
| CORT | Corticosterone (rodents, avian, and non-avian reptiles) or cortisol (non-rodent mammals) |
| GC | Glucocorticoid hormones such as corticosterone/cortisol |
| GLS | Generalized least squares |
| HPA | Hypothalamic–pituitary–adrenal axis |
| HSP | Heat-shock proteins. Number and letters associated with “HSP” denote a specific protein with names following Telemeco et al. (2019) |
| MAN(C)OVA | Multivariate analysis of variance or covariance. When used in isolation, refers to traditional parametric methods |
| NAb | Measure of natural antibody ability to agglutinate heterologous red blood cells |
| NP-MANOVA | Non-parametric multivariate analysis of variance sensuAnderson (2001) |
| OLS | Ordinary least squares |
| PCA | Principal components analysis |
| PGLS | Phylogenetic generalized least squares |
| RRPP | Residual randomization in permutation procedures. When in quotes (“RRPP”) refers to R package of the same name |
| SS | Sum of squares |
| Abbreviation | Definition |
|---|---|
| AN(C)OVA | Analysis of variance or analysis of covariance |
| CH50 | Measure of plasma complement ability to lyse heterologous red blood cells independent of antibodies |
| CORT | Corticosterone (rodents, avian, and non-avian reptiles) or cortisol (non-rodent mammals) |
| GC | Glucocorticoid hormones such as corticosterone/cortisol |
| GLS | Generalized least squares |
| HPA | Hypothalamic–pituitary–adrenal axis |
| HSP | Heat-shock proteins. Number and letters associated with “HSP” denote a specific protein with names following Telemeco et al. (2019) |
| MAN(C)OVA | Multivariate analysis of variance or covariance. When used in isolation, refers to traditional parametric methods |
| NAb | Measure of natural antibody ability to agglutinate heterologous red blood cells |
| NP-MANOVA | Non-parametric multivariate analysis of variance sensuAnderson (2001) |
| OLS | Ordinary least squares |
| PCA | Principal components analysis |
| PGLS | Phylogenetic generalized least squares |
| RRPP | Residual randomization in permutation procedures. When in quotes (“RRPP”) refers to R package of the same name |
| SS | Sum of squares |
List of abbreviations and acronyms
| Abbreviation | Definition |
|---|---|
| AN(C)OVA | Analysis of variance or analysis of covariance |
| CH50 | Measure of plasma complement ability to lyse heterologous red blood cells independent of antibodies |
| CORT | Corticosterone (rodents, avian, and non-avian reptiles) or cortisol (non-rodent mammals) |
| GC | Glucocorticoid hormones such as corticosterone/cortisol |
| GLS | Generalized least squares |
| HPA | Hypothalamic–pituitary–adrenal axis |
| HSP | Heat-shock proteins. Number and letters associated with “HSP” denote a specific protein with names following Telemeco et al. (2019) |
| MAN(C)OVA | Multivariate analysis of variance or covariance. When used in isolation, refers to traditional parametric methods |
| NAb | Measure of natural antibody ability to agglutinate heterologous red blood cells |
| NP-MANOVA | Non-parametric multivariate analysis of variance sensuAnderson (2001) |
| OLS | Ordinary least squares |
| PCA | Principal components analysis |
| PGLS | Phylogenetic generalized least squares |
| RRPP | Residual randomization in permutation procedures. When in quotes (“RRPP”) refers to R package of the same name |
| SS | Sum of squares |
| Abbreviation | Definition |
|---|---|
| AN(C)OVA | Analysis of variance or analysis of covariance |
| CH50 | Measure of plasma complement ability to lyse heterologous red blood cells independent of antibodies |
| CORT | Corticosterone (rodents, avian, and non-avian reptiles) or cortisol (non-rodent mammals) |
| GC | Glucocorticoid hormones such as corticosterone/cortisol |
| GLS | Generalized least squares |
| HPA | Hypothalamic–pituitary–adrenal axis |
| HSP | Heat-shock proteins. Number and letters associated with “HSP” denote a specific protein with names following Telemeco et al. (2019) |
| MAN(C)OVA | Multivariate analysis of variance or covariance. When used in isolation, refers to traditional parametric methods |
| NAb | Measure of natural antibody ability to agglutinate heterologous red blood cells |
| NP-MANOVA | Non-parametric multivariate analysis of variance sensuAnderson (2001) |
| OLS | Ordinary least squares |
| PCA | Principal components analysis |
| PGLS | Phylogenetic generalized least squares |
| RRPP | Residual randomization in permutation procedures. When in quotes (“RRPP”) refers to R package of the same name |
| SS | Sum of squares |
References
R Core Team.
