Testing hypotheses in ecoimmunology using mixed models: disentangling hierarchical correlations.

Integrative and Comparative Biology Advance Access published May 26, 2014

Integrative and Comparative Biology Integrative and Comparative Biology, pp. 1–12 doi:10.1093/icb/icu035

Society for Integrative and Comparative Biology

SYMPOSIUM

Testing Hypotheses in Ecoimmunology Using Mixed Models: Disentangling Hierarchical Correlations C. J. Downs*,1 and N. A. Dochtermann† *Department of Natural Resources and Environmental Sciences, University of Nevada, Reno, NV 89557, USA; † Department of Biological Sciences, North Dakota State University, Fargo, ND 58102, USA Both authors contributed equally to this work.

1

E-mail: [email protected]

Synopsis Considerable research in ecoimmunology focuses on investigating variation in immune responses and linking this variation to physiological trade-offs, ecological traits, and environmental conditions. Variation in immune responses exists within and among individuals, among populations, and among taxonomic groupings. Understanding how variation and covariation are distributed and how they differ across these levels is necessary for drawing appropriate ecological and evolutionary inferences. Moreover, variation at the among-individual level directly connects to underlying quantitative genetic parameters. In order to fully understand immune responses in evolutionary and ecological contexts and to reveal phylogenetic constraints on evolution, statistical approaches must allow (co)variance to be partitioned among levels of individual, population, and phylogenetic organization (e.g., population, species, genera, and various higher taxa). Herein, we describe how multi-response mixed-effects models can be used to partition variation in immune responses among different hierarchical levels, specifically within-individuals, among-individuals, and among-species. We use simulated data to demonstrate that mixed models allow for proper partitioning of (co)variances. Importantly, these simulations also demonstrate that conventional statistical tools grossly misestimate relevant parameters, which urges caution in relating ecoimmunological hypotheses to existing empirical research. We conclude by discussing the advantages and caveats of a mixed-effects modeling approach.

Introduction Life histories tend to fall on a gradient of slow paces to fast paces (Stearns 1992; Ghalambor and Martin 2001). Slow-paced life histories are characterized by longer lifespan, higher adult survival, longer development time, and lower reproductive rate; conversely fast-paced life histories are characterized by shorter lifespans, shorter development times, and higher reproductive rates (Saether 1988; Promislow and Harvey 1990; Ricklefs 2000). Recent attempts to explain this pattern have focused on how natural selection acts on life histories, as well as on the physiological and genetic mechanisms underlying tradeoffs among life-history characteristics (Zera and Harshman 2001; Ricklefs and Wikelski 2002). Physiology, particularly immune function, endocrine function, and metabolic rates, elucidates mechanisms

that underlie life-history trade-offs (Sibly and Calow 1986; Ricklefs and Wikelski 2002). The pace-of-life hypothesis integrates environmental influences on life-history traits with the physiological mechanisms that underlie life-history trade-offs, to explain why most life histories fall along a continuum of slowto-fast pace (Ricklefs and Wikelski 2002). The pace-of-life hypothesis posits that investment in immune function should vary in a predictable manner across the pace-of-life continuum (Ricklefs and Wikelski 2002; Lee, 2006). Predictions related to the pace-of-life hypothesis generally divide immune responses along two axes: specific to non-specific responses and constitutive to induced responses (Schmid-Hempel and Ebert 2003; Lee 2006). Nonspecific immunity is characterized by rapid, nonantigen-specific responses including the complement

ß The Author 2014. Published by Oxford University Press on behalf of the Society for Integrative and Comparative Biology. All rights reserved. For permissions please email: [email protected].

Downloaded from http://icb.oxfordjournals.org/ at UQ Library on May 31, 2014

From the symposium ‘‘Methods and Mechanisms in Ecoimmunology’’ presented at the annual meeting of the Society for Integrative and Comparative Biology, January 3–7, 2014 at Austin, Texas.

2

Why might we expect different relationships among traits at different levels? Allocation versus acquisition as an example An underlying assumption of most hypotheses related to life-history theory, and thus ecoimmunology, is the principle of allocation. The principle of allocation posits that resources are limited so individuals must choose how to allocate resources among functions essential for survival and reproduction (Stearns 1992). Decisions about allocation depend on varying

selective pressures and are expressed through variation in investment in traits among individuals (Ricklefs and Wikelski 2002). The general prediction of the allocation hypothesis is that traits that share a common resource should be negatively correlated with each other across individuals because resources invested in one trait are unavailable for another trait (Williams 1966; Stearns 1992). The principle of allocation assumes that the quantity and types of resources available to each individual are the same. This assumption is invalid, however, because acquisition of resources is affected by numerous factors including habitat, rank, genetic background, maternal effects, and state (Flatt and Heyland 2011). Differential acquisition posits that one might expect any type of correlation between life-history traits among individuals, including a counterintuitive positive correlation, because individuals acquire different amounts of resources. A positive correlation might arise among individuals because individuals with more resources overall can allocate more resources to all functions (Van Noordwijk and De Jong 1986). A classic illustration is the big house, big car paradox (Van Noordwijk and De Jong 1986; Reznick et al. 2000). The allocation hypothesis would predict that if money were limited, then families that allocate money to buy a big house should not have enough money to also buy a big car. Instead, a positive correlation is often seen between cars and houses across families because families differ in income, that is, the total amount of money they have available to spend (Van Noordwijk and De Jong 1986; Reznick et al. 2000). Despite individuals with big houses also having big cars, it still holds that when resource levels are limited within an individual, that individual must reduce resources allocated to one or more traits to increase allocation to another trait. When comparing across individuals that have access to different amounts of resources, however, a positive correlation between costly traits is possible with differential acquisition or other differences in state. Thus, the combination of the principle of allocation and differential acquisition suggest that within an individual, a negative relationship is expected between traits, whereas among individuals a positive relationship is expected among traits (Fig. 1A and B). As a result, when correlations between traits are estimated without properly considering within-individual and among-individual components, inferences about patterns that should be manifested at either of these levels can be incorrect. Herein we demonstrate how multi-response mixed-effects models also being discussed in other areas of evolutionary ecology (e.g.,


pathway and the acute-phase response. Specific immunity is characterized by slower, antigen-specific responses, including those from B-cell and mammalian toll-like receptors (Schmid-Hempel and Ebert 2003; Lee 2006). The second axis forms a continuum from constitutive to induced immune responses. Constitutive responses are always present and capable of immediate physiological defense, whereas induced responses are deployed after the pathogen has been recognized (Schmid-Hempel and Ebert 2003; Lee 2006). The energetic, resource, and self-damage costs associated with an immune response, and therefore the selective pressures and predictions about investment into immune function, depend on the type and strength of the immune response. Lee (2006) developed testable hypotheses about how investment in different types of immune function should vary across the pace-of-life continuum at the interspecific level. For example, slow-paced species should invest more in specific, induced immune responses, such as antibody responses, than do fastpaced species, because these types of immune response require substantial investment and are most beneficial against repeated infections that are more likely in long-lived, slow-paced species (Klasing and Leshchinsky 1999; Lee 2006). Hypotheses also extend to the individual level to develop predictions about how investment in immune function should vary among individuals and populations depending on other energetic requirements, including reproductive status, sex, latitude, and age (Lee 2006). At this level of comparison, individuals that invest more intensively in activities such as reproduction should invest less in specific, induced immune responses than do individuals that invest less effort because reducing investment in immune function will free up resources for other activities (Lochmiller and Deerenberg 2000; Lee 2006). These hypotheses are the basis for much research in ecoimmunology, and a statistical framework is needed to test predictions simultaneously among species, among individuals, and within individuals.

C. J. Downs and N. A. Dochtermann

3

Mixed models for ecoimmunology

Objectives and an example

Fig. 1 At the level of within individuals, the principle of allocation suggests a trade-off between investments in basal metabolic rate and antibody response; thus, we predict a negative correlation for within-individuals (A). At the among-individual level, differential acquisition would suggest that individuals which acquired more resources have more to invest both in immune function and in BMR; thus, we predict a positive relationship between these traits for among-individuals (B). Finally, based on the pace-of-life hypothesis we predict a negative correlation between these traits at the interspecific level (i.e., among species; C).

Dingemanse et al. 2012; Dingemanse and Dochtermann 2013) can be applied to test ecoimmunological hypothesis that make predictions about the correlations between traits within individuals and among individuals. We will also test other statistical

For simplicity, in this article we will discuss how to statistically model a single trade-off predicted by the pace-of-life hypotheses: energy used for basal metabolic rate (BMR) and investment in antibody production. Please note that we could have chosen for our example any two traits that are expected to be correlated. The predictions developed below for BMR and antibody production are intended to provide a concrete example to illustrate the statistical models presented herein. Although the hypotheses are based in theory, other rationale could lead to different predictions. We do not test the predictions below with empirical data (necessary data have not yet been collected by ecoimmunologists) but rather we simulated data to test the effectiveness of different statistical models to properly estimate within-individual, among-individual, and among-species correlations. Antibody production is an induced, specific response. As such, investment in antibody response is expected to decrease as pace-of-life becomes faster when compared across species, and increase with intensity of effort when compared across individuals (Lee 2006). BMR is the minimum energy required


frameworks for estimating the correlation of traits among individuals, including raw phenotypic correlations, and a correction factor developed by Adolph and Hardin (2007). A similar problem arises when trying to quantify how traits are correlated across species. If numerous individuals from a species are measured to develop a mean value for that species, correlations between two traits that are developed from these means confound interspecific and intraspecific variation. Instead, variation in traits needs to be partitioned between these two types of variation to accurately estimate the correlations. Note that in this discussion, intraspecific variation is the same as among-individual variation; that is, both terms refer to the variation in a trait assigned to the differences among individuals of the same species. It should be noted, that if each individual is not measured multiple times, then the intraspecific term will include variation that arises from variation within and among individuals. Herein we extend the multi-response mixed-effects models to partition correlations between traits into among-individual correlations and among-species correlations. We also present a statistical model to partition correlations simultaneously to the component parts of the within-individual correlation, the among-individual correlation, and the among-species correlation.

4


The problem When theory predicts correlations of different signs (or magnitudes) across levels of organization, how might these predictions be evaluated? This issue can be addressed based on knowledge of how different factors contribute to phenotypic correlations derived from measuring numerous individuals. Although potentially much more complex (Fig. 2), phenotypic correlations stem from the contribution of among-individual correlations and withinindividual correlations (Dingemanse et al. 2012; Dingemanse and Dochtermann 2013): vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ! u u Vindy Vind z rPy ,Pz ¼ rindy ,indz t VPy V Pz vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ! u u VPy Vindy VP Vind z z þ rey ,ez t V Py VPz

ð1Þ

where rPy ,Pz is the phenotypic correlation between traits y and z, rindy ,indz is the among-individual correlation for y and z, rey ,ez is the within-individual correlations for y and z, VPy and VPz are the phenotypic variances of y and z, respectively; Vindy and Vindz are the among-individual variances of y and z, respectively (Dingemanse et al. 2012; Dingemanse and Dochtermann 2013). Conceptually, we can think about among-individual estimates of variances as how much individuals differ in their average

phenotypes while within-individual variances estimate how much an individual differs in its phenotypic expression from one measurement to the next. Moreover, the among-individual variance for a trait divided by the total phenotypic variance for that trait Vind (e.g., VP y ) is repeatability. Likewise, rindy ,indz concepy tually would be the correlation between individuals’ mean y’s and z’s. Conceptually the within-individual correlation, rey ,ez , is then the average correlation across multiple individuals when y and z are measured multiple times in each individual and the correlation for that individual is estimated. Importantly, ecoimmunological theory makes specific predictions about trade-offs corresponding to both correlations in Equation (1). Specifically, allocation theory estimates that trade-offs will be manifested at the rey ,ez level. However, differences in states, differences in genetic propensities in allocation, and differences in permanent environmental effects will be manifested at the rindy ,indz level. As these correlations are conflated at the phenotypic level, they cannot be reliably estimated by simple phenotypic correlations. For example, the phenotypic correlation will be zero if y and z are each measured multiple times, have repeatabilities of 0.5, and have a correlation of 0.5 at the withinindividual level and a correlation of 0.5 at the among-individual level. Fortunately this conflation can be addressed by separately estimating the within-individual and among-individual correlations, which can be accomplished using multi-response mixed-effects models whenever multiple measurements at each level are estimated. Moreover, other ecoimmunological and general evolutionary ecological theories make predictions at various hierarchical levels such as amongindividuals and among species. When there is replication within levels, mixed models can be used to estimate correlations at these levels as well. Here we provide a simple demonstration of disentangling correlations for the among-individual and within-individual levels, for the among-individual and among-species levels, and for all three levels simultaneously.

A potential solution: multi-response mixed-effects models Disentangling among-individual and within-individual correlations Recall the example described previously: for two traits, such as BMR (y) and antibody production (z), a negative allocation trade-off is predicted within individuals (Fig. 1A) and the positive


for a resting, post-absorptive animal that is not growing, reproducing, or thermoregulating (McNab 1997). At the level of within-individuals, the principle of allocation suggests that if an individual has a given amount of resources and uses more energy for BMR, it will have fewer resources available to mount an antibody response; thus a negative correlation is predicted within-individuals (Fig. 1A). The amongindividual prediction is less clear, but differential acquisition would suggest that individuals which acquired more energy have more available both for immune function and for BMR. Thus, a positive relationship between these traits is predicted amongindividuals (Fig. 1B). Finally, slow-paced species have lower BMRs than do fast-paced species (Ricklefs and Wikelski 2002; Wiersma et al. 2007), though slowpaced species should also invest more in antibody responses than do fast-paced species for reasons explained previously (Klasing and Leshchinsky 1999; Lee 2006). Thus, a negative correlation between these traits is predicted at the interspecific level (i.e., among species; Fig. 1C).

5


relationship is expected among individuals because of differences in acquisition of resources (Fig. 1B). In such a case, phenotypic correlations can be misleading and these correlations must be estimated separately. When both traits are measured at the same time and each is measured repeatedly, this can be done with a multi-response mixed-effects model (Dingemanse and Dochtermann 2013). Such a model for traits y and z has the form: yij ¼ 0y þ ind0yj þ e0yij zij ¼ 0z þ ind0zj þ e0zij

ð2Þ

where yij and zij are BMR (y) and immune response (z) measured for individual j at sampling time i. This model is fit by estimating a random intercept (ind0j ) for each individual, that is, ‘‘Individual’’ is a random factor. These random intercepts, for a model like Equation (2), can be thought of as an individual’s deviation (with a mean of zero) from the population mean, 0y and 0z . These population means are typically treated as independent (e.g., Matsuyama and Ohashi 1997). e0yij and e0zij are the residual or within-individual deviations from this model for individual j at sampling time i and are estimated as having averages of zero. What differentiates a multi-response mixed-effects model of Equation (2) from more common mixed

models is that the random intercepts and the withinindividual variances for each trait are not estimated as independent. Instead, both are fit to a single multivariate normal distribution. These estimated distributions have—for two traits—two variances and a covariance estimated from the data. The variances and covariances estimated among the random intercepts are among-individual variances and the among-individual covariance. The variances and covariances among the within-individual deviations are the within-individual variances and covariance: ind0yj MVN ð0,ind Þ : ind0zj " # ð3Þ Vind0y ind ¼ Covind0y ,ind0z Vind0z

e0yj e0zj

"

MVN ð0,e Þ :

e ¼

#

Ve0y Cove0y ,e0z

Ve0z

From the among-individual and within-individual variances and covariances we can calculate rind0y ,ind0z and re0y ,e0z as [Equations (4) and (5)]: COVind0y ,ind0z rind0y ,ind0z ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Vind0y Vind0z

ð4Þ

COVe0y ,e0z re0y ,e0z ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ve0y Ve0z

ð5Þ


Fig. 2 Hierarchical diagram illustrating how genetic effects, environmental effects, and measurement error underpin between-individual and within-individual correlations that underpin raw phenotypic correlations (originally published by Dingemanse and Dochtermann 2013).

6


An example of jointly estimating among-individual and within-individual correlations To illustrate the utility of this approach, recall our predicted relationships in Fig. 1A and B. How do different approaches to estimating correlations, at different sample sizes, perform in estimating the hypothesized relationships? To address this question we simulated data according to the following: yij ¼ 0 þ ind0yj þ e0yij zij ¼ 0 þ ind0zj þ e0zij

ind0yj ind0zj

The mean within-individual correlation for this model was 0.5, and the mean among-individual correlation was 0.5. Based on these distributional relationships we generated datasets with two measures per trait, per individual. Datasets were generated with population sizes of 25, 50, 75, 100, 125, and 150 individuals and 200 datasets for each size of population. This approach to the generation of data (used throughout) conflates sampling error with estimation error. However, we consider this conflation to accurately reflect how real-world data are collected. For each of these 1200 datasets we calculated the correlation between y and z in the following ways: (1) based on a single measure of each trait for each individual; (2) based on the average of each trait by individual; (3) based on Adolph and Hardin’s (2007) repeatability-based correction of individual mean-based correlations; and (4) based on the described multi-response mixed-effects approach to modeling. Multi-response mixed-effects models were fit using the MCMCglmm package (Hadfield 2010) for the R programing language (R Core Development Team 2013). ‘‘Individual’’ was included in the model as a random factor. Models were fit with a prior that was flat for the correlation and using 13,000 iterations, a 3000 burn-in, and a thinning interval of 10. For the various estimates of correlation, we evaluated how well each represented underlying parameters of interest.

Estimates of the correlation between y and z were generally poor when not made specifically at the within-individual or among-individual levels (Fig. 3). As correlations based both on single measures of individuals and individual-level means conflated within-individual and among-individual correlations, they performed very poorly at identifying trade-offs (Fig. 3). If inferences were drawn from single estimates per individual, then researchers would conclude from our simulated dataset that no trade-offs are present. Although correlations across individual means performing better at identifying among-individual correlations, they nonetheless underestimate values used to generate the populations (i.e., generating values) (Fig. 3). Of non-mixed model approaches, Adolph and Hardin’s (2007) weighting of correlations based on repeatability most closely approximated the generating values. However, this value still dramatically underestimated the generating values (Fig. 3). The inability of Adolph and Hardin’s (2007) measure to properly estimate among-individual correlations is likely due to its calculation necessitating the assumption that withinindividual correlations are zero, an assumption we know to be violated for this simulated dataset. A multi-response mixed-effects modeling approach most accurately estimated among-individual correlations. This approach was also the only approach to identify the pattern of within-individual correlations that underlaid observed phenotypes (Fig. 3). These results strongly suggest that mixed-model approaches may be necessary to ask appropriate questions in ecoimmunology. An alternative approach, known as within-individual centering (van de Pol and Wright 2009), would also require the collection of repeated measures and the use of mixed models. The approach we describe here is preferable, as within-individual centering requires researchers to assume that one of their two measures was estimated with no error. Discerning patterns among species Multi-response mixed models can also be used to partition variation across other hierarchical levels, making the approach useful for testing hypotheses derived from ecoimmunological theory that describe patterns of correlations of traits across species. This includes questions from the pace-of-life hypothesis, such as how does investment in specific versus nonspecific immune responses or induced versus constitutive immune responses differ among species along the pace-of-life continuum (Schmid-Hempel and


MVN ð0,ind Þ : " # Vind0y ¼ 1 ind ¼ Covind0y ,ind0z ¼ 0:5 Vind0z ¼ 1 e0yj MVN ð0,e Þ : e0zj " # Ve0y ¼ 1 e ¼ Cove0y ,e0z ¼ 0:5 Ve0z ¼ 1

Results

7


Ebert 2003; Lee 2006)? For such questions researchers are asking among-species (i.e. interspecific) questions under which among-individual and withinindividual patterns will be nested. As was the case among individuals, this nesting of patterns results in a conflation of various contributors to an observed correlation (rObsy ,Obsz ): vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u u Vspecy Vspec z t rObsy ,Obsz ¼ rspecy ,specz VObsy V Pz vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u u Vindy Vind z t þ rindy ,indz VObsy VObsz

ð6Þ

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ! u u VObs Vspec Vind y y y u u VObsy u þ rey ,ez u u VObsz Vspecz Vindz t VObsz

where all values are the same as in Equation (1), and where rspecy ,specz is the among-species correlation (species: spec) and Vspecy and Vspecz are the among-species variances for y and z, respectively.

Above we considered a example where the withinindividual and among-individual predictions differ (Fig. 1A and B); we can add to these expected relationships a prediction that BMR (y) and production of antibodies (z) should be negatively correlated at the among-species level for reasons described previously (Fig. 1C). Hypotheses in comparative physiology and ecoimmunology often make predictions about patterns among species begging the question: how effective are different methods, including multiresponse mixed-effects models, in estimating the correlations of interest? To address this question, we again simulated data incorporating correlations at different levels. The structure of the model was: yhij ¼ 0 þ spec0yh þ ind0yj þ e0yhij zhij ¼ 0 þ spec0zh þ ind0zj þ e0zhij

spec0yh

spec0zh spec

MVN 0,spec : " Vspec0y ¼ 1 ¼ Covspec0y ,spec0z ¼ 0:5

# Vspec0z ¼ 1


Figure 3. Among-individual correlations (mean 95% CI) estimated by single measures per individual (black circles), individual means (white circles), Adolph and Hardin’s (2007) corrected estimate (black squares), and mixed models (black diamonds). The mixed-model framework provided the most accurate estimate of the among-individual correlation, and it was the only method that provided an estimate of the within-individual correlation (mean 95% CI; white triangles). Correlations were estimated from simulated data (see text for details) for different sample sizes of individuals (n ¼ 25, 50, 75, 100, 125, 150). Each individual had two simulated measurements for each trait.

8


ind0yj

MVN ð0,ind Þ : " # Vind0y ¼ 1 ind ¼ Covind0y ,ind0z ¼ 0:5 Vind0z ¼ 1 e0yij MVN ð0,e Þ : e0zij " # Ve0y ¼ 1 e ¼ Cove0y ,e0z ¼ 0:5 Ve0z ¼ 1 ind0zj

Results Correlations based either on species’ means or single measures per species produced correlation estimates that were substantially lower than the actual populations-among-species correlation (Fig. 4). A mixed modeling approach generally provided the best

Integrating within-individual, among-individual, and among-species questions In the previous two examples we estimated correlations at levels that are likely to interest ecoimmunological researchers: (1) within-individual and betweenindividual relationships (e.g., Demas et al. 1997; Tieleman et al. 2005) or (2) among population or interspecific/among-species relationships (e.g., Addison et al. 2009; Downs et al. 2012, 2013; Tieleman et al. 2005, Versteegh et al. 2012). However, the more complete representation of contributors to within-individual, among-individual, and among-species patterns represented by Equation (6) and Fig. 1A–C (in total) has not been tested. Although researchers studying patterns for a species might collect multiple measurements per individual, researchers studying patterns across species instead typically will try to maximize the number of species they examine and, perhaps, the number of individuals per species. When asking questions across species, repeated measurements per individual are likely to be logistically prohibitive. However, if individuals are measured repeatedly and multiple individuals are measured per species, then the whole suite of correlations in Equation (6) can be estimated separately. To illustrate this scenario, we used the same model structure we used for examining among-species patterns, but generated data for 150 species with


where values are as above and spec0yh and spec0zh correspond to a random intercept for each species and each trait. The mean within-individual correlation for this model was 0.5, the mean amongindividuals correlation was 0.5, and the mean among-species correlation was 0.5. Based on this model we generated datasets with two individuals per species and one measure per individual. Datasets were generated with 25, 50, 75, 100, 125, and 150 species and 200 datasets per level of number of species. For each of these 1200 datasets we calculated the correlation between y and z in the following ways: (1) based on a single measure of each trait for each species; (2) based on the average of each trait by species; (3) based on the described multi-response mixed-effects modeling approach. Multi-response mixed-effects models were fit using the MCMCglmm package (Hadfield 2010) for the R programing language (R Core Development Team 2013). ‘‘Species’’ but not ‘‘individual’’ was included in the model as a random factor. Both ‘‘species’’ and ‘‘individual’’ were used in the model to generate the datasets used to test the different statistical models, but we included only ‘‘species’’ as a random effect within analyses because we are analyzing datasets without repeated measurements of individuals. To model within-individual correlations, we would need more than one measurement for each individual and to include ‘‘individual’’ as a random effect (see Section ‘‘Integrating within-individual, among-individual, and among-species questions’’ for details). Models were fit with a prior that was flat for the correlation and using 13,000 iterations, a 3000 burn-in, and a thinning interval of 10. For the various correlation estimates, we evaluated how well each represented underlying parameters of interest.

estimate of the among-species correlation, although 95% confidence intervals (CIs) did not quite overlie the generating distributions of known correlation (Fig. 4). This is likely due to sampling error, that is, differences between the simulated sample population and the generating distributions. Importantly, as we modeled the data based on a single measure of each trait for each individual sampled in a species, we were not able to separately estimate the withinindividual and among-individual correlations. As a result, the within-species correlation conflated these two effects and was therefore estimated as 0 (Fig. 4). As in the case of patterns within and across individuals, our simulation suggests that multi-response mixed-effects models are strongly recommended for asking questions at the among-species level. It should also be noted that this framework could be used to compare populations-level differences if ‘‘population’’ replaced ‘‘species’’ in Equation (6). Indeed the flexibility of the multi-response mixedeffects models is illustrated by the fact that if all of the terms for the species level are removed, Equation (6) becomes Equation (1).

9


three individuals measured twice each for both traits. Three individuals per species were chosen in accordance with current recommendations for the minimum for obtaining good estimates of BMRs (McKechnie and Wolf 2004). In general, sampling more individuals should improve the accuracy of the estimates of the correlation (Dingemanse and Dochtermann 2013). Exploring how the number of individuals sampled per species and how the characteristics of the individuals sampled (e.g., sex, state, and age) changes the estimates is beyond the scope of this article. We fit a multi-response mixed-effects model to the data, including ‘‘species’’ and ‘‘individual’’ both as random factors (individual nested within species). The data were fit with a prior that was flat for the correlation and using 13,000 iterations, a 3000 iteration burn-in, and a thinning interval of 10. We did not explore additional combinations of parameters because varying number of species, number of individuals, and number of repeated measures per individual would lead to an overwhelming number of combinations of parameters. The joint estimation of all three correlations using a multi-response mixed-effects model led to quite accurate estimates (Fig. 5). All three correlations

were quite close to the generating distribution’s parameters. This suggests that the described modeling approach could be very effective in asking integrative questions about correlations of traits across levels of organization. However, these results and this particular simulation come with very strong caveats. Our analysis assumes that variances and covariances at the among-individual and within-individual levels do not differ across species. In other words, patterns and trade-offs for individuals and state-trait relationships among individuals are being estimated as the same for all species. This is a highly unlikely assumption and strongly informs the estimates of variance and covariance. Statistical models can be formulated that would allow unique estimation of variance and covariance on a species basis. However, this approach would require a large number of measurements of both traits for each individual and a large number of unique individuals. Although these conditions are unlikely for free-ranging animals, it might be possible to meet these conditions to conduct a cross-species comparison using captive animals.

Conclusions Mixed models provide better estimates of relevant correlations than do alternative approaches.


Figure 4. Between-species correlations (mean 95% CI) estimated by single measures per species (black circles), species’ means (white circles), and mixed models (black diamonds). The mixed-model framework provided the most accurate estimate of the amongspecies correlation, and it was the only method that provided an estimate of the within-species correlation (white triangles). Note that the within-species correlation includes variation from within-individual and among-individual correlations. Correlations were estimated from simulated data (see text for details) for different samples sizes of species (n ¼ 25, 50, 75, 100, 125, 150). Each species had simulated measurements of two individuals for both traits.

10


Multi-response mixed-effects models can be very data-hungry. If sample sizes are not sufficiently large, results and the conclusions drawn from those results are unreliable. Figs 3, 4, and those within Dingemanse and Dochtermann (2013) provide additional guidance for sampling schemes.

Fig. 5 Within-individual, among-individual, and among-species correlations (mean 95 CI) estimated from simulated data for 150 species with three individuals sampled twice for both traits for each species.

When estimating phenotypic correlations across individuals, the mixed model is the only approach that provides an estimate of the within-individual correlation. This same modeling approach also provides the most accurate estimates of the among-individual correlation. When estimating correlations across species, a mixed modeling framework provides more accurate estimates of species’ correlation because it is accounting for the variation caused by individual variation. These models do not estimate the among-individual or within-individual correlations when only single measures per individual are available. Among-individual or withinindividual correlations are conflated as part of the within-species correlation. If repeated measures per individual and repeated individuals per species are both available, all components of variance – within-individual, among-individual, among-species variances and correlations – can be estimated. However, the assumptions underlying the example and the structure of the model provided herein are unlikely and there are currently no empirical estimates to support these assumptions. To properly account for differences in trade-offs among species requires much greater sample sizes and more complicated mixed models.

Many hypotheses about ecoimmunological tradeoffs assume that immune responses are nutritionally or energetically costly and invoke the principle of allocation (e.g., Gustafsson et al. 1994; Lochmiller and Deerenberg 2000; Downs et al. 2013). Hypotheses regarding trade-offs caused by decisions about allocations will be manifested within the individual and require repeated measures to be tested. In contrast, hypotheses about the relationship between traits that are caused by differences in states, differences in genetic propensities in allocation, and differences in permanent environmental effects will be manifested among individuals. The overall phenotypic correlation must be partitioned into within-individual and among-individual components to assess these hypotheses. Similarly, there is substantial variation among individuals of the same species in physiological traits, including immune function (Hayes and Jenkins 1997; Schmid-Hempel 2003; Speakman et al. 2004). To properly test ecoimmunological hypotheses that make predictions about patterns across species, such as those from the pace-oflife hypothesis (Ricklefs and Wikelski 2002; Lee 2006), statistical models must properly account for individual-level correlation (i.e., the combination of within-individual correlation and among-individual correlation). The approach we have described applies not only to ecoimmunological trade-offs; rather it applies to the evaluation of the relationships of traits regardless of trait domain. Examples herein have focused on predictions for which the signs of the correlation are conflicting


If researchers are investigating particular ecoimmunological trade-offs, simple correlations often may be misleading. This may be contributing to current debates regarding the actual role of tradeoffs in natural populations. Specifically, simple correlations are not necessarily testing for the presence of trade-offs, an issue that has also been discussed extensively in life-history and quantitative genetics research (e.g., Roff and Fairbairn 2007).


among levels. If the sign of the correlations are the same, then estimates of the correlations are still conflated and the relative contributions of effects cannot be determined unless a mixed-model approach is used. Models for the comparisons of species assume star-phylogenies. To correct for phylogenetic relationships, phylogenetically independent contrasts can be used or information about branch lengths in phylogenetic trees can be incorporated into the mixed-model framework (Hadfield and Nakagawa 2010; Nakagawa and Santos 2012).

Acknowledgments We thank the Society for Integrative and Comparative Biology’s Divisions of Animal Behavior, Comparative Endocrinology, and Comparative Physiology and Biochemistry and the Research Coordination Network for Ecoimmunology funded by the National Science Foundation [IOS 094177] for supporting the symposium ‘‘Methods and Mechanisms in Ecoimmunology’’ financially. We thank H. Heatwole and two anonymous reviewers for helpful comments on an early draft of this manuscript.

Funding This work was supported by Hatch grant awarded to K.M. Stewart by the University of Nevada Agriculture Experimental Station & the United States Department of Agriculture (to C.J.D.); the ND-EPSCoR program, the NDSU Department of Biological Sciences, and the NDSU College of Science and Mathematics (to N.A.D.).

References Addison B, Klasing KC, Robinson WD, Austin SH, Ricklefs RE. 2009. Ecological and life-history factors influencing the evolution of maternal antibody allocation: a phylogenetic comparison. P Roy Soc B Bio 276:3979–87. Adolph SC, Hardin JS. 2007. Estimating phenotypic correlations: correcting for bias due to intraindividual variability. Funct Ecol 21:178–84. Demas GE, Chefer V, Talan MI, Nelson RJ. 1997. Metabolic costs of mounting an antigen-stimulated immune response

in adult and aged C57BL/6J mice. Am J Physiol Regul Integr Comp Physiol 273:R1631–7. Dingemanse NJ, Dochtermann NA. 2013. Quantifying individual variation in behaviour: mixed-effect modelling approaches. J Anim Ecol 82:39–54. Dingemanse NJ, Dochtermann NA, Nakagawa S. 2012. Defining behavioural syndromes and the role of ‘syndrome deviation’ in understanding their evolution. Behav Ecol Sociobiol 66:1543–8. Downs CJ, Brown JL, Wone B, Donovan ER, Hunter K, Hayes JP. 2013. Selection for increased mass-independent maximal metabolic rate suppresses innate but not adaptive immune function. P Roy Soc B Bio 280:20122636. Downs CJ, Schutz H, Meek TH, Dlugosz EM, Acosta W, de Wolski KS, Malisch JL, Hayes JP, Garland T. 2012. Withinlifetime trade-offs but evoutionary freedom for hormonal and immunological traits: evidence from mice bred for high voluntary exercise. J Exp Biol 215:1651–61. Flatt T, Heyland A. 2011. Mechanisms of life history evolution: the genetics and physiology of life history traits and trade-offs. New York, NY: Oxford University Press. Ghalambor CK, Martin TE. 2001. Fecundity-survival trade-offs and parental risk-taking in birds. Science 292:494–7. Gustafsson L, Nordling D, Andersson MS, Sheldon BC, Qvarnstrom A. 1994. Infectious diseases, reproductive effort and the cost of reproduction in birds. Philos T Roy Soc B 346:323–31. Hadfield JD. 2010. MCMC methods for multi-response generalized linear mixed models: the MCMCglmm R package. J Stat Softw 33:1–22. Hadfield JD, Nakagawa S. 2010. General quantitative genetic methods for comparative biology: phylogenies, taxonomies and multi-trait models for continuous and categorical characters. J Evol Biol 23:494–508. Hayes JP, Jenkins SH. 1997. Individual variation in mammals. J Mammal 78:274–93. Klasing KC, Leshchinsky TV. 1999. Functions, costs, and benefits of the immune system during development and growth. Johannesburg: BirdLife South Africa. Lee KA. 2006. Linking immune defenses and life history at the levels of the individual and the species. Integr Comp Biol 46:1000–15. Lochmiller RL, Deerenberg C. 2000. Trade-offs in evolutionary immunology: just what is the cost of immunity? Oikos 88:87–98. Matsuyama Y, Ohashi Y. 1997. Mixed models for bivariate response repeated measures data using Gibbs sampling. Stat Med 16:1587–601. McKechnie AE, Wolf BO. 2004. The allometry of avian basal metabolic rate: good predictions need good data. Physiol Biochem Zool 77:502–21. McNab BK. 1997. On the utility of uniformity in the definition of basal rate of metabolism. Physiol Zool 70:718–20. Nakagawa S, Santos ESA. 2012. Methodological issues and advances in biological meta-analysis. Evol Ecol 26:1253–74. Promislow DEL, Harvey PH. 1990. Living fast and dying young: a comparative analysis of life-history variation among mammals. J Zool 220:417–37.


Although we have focused on two traits, this mixed-model framework can be expanded to include more traits if multivariate models are used (Dingemanse and Dochtermann 2013). However, as the number of traits included in the model increases, the samples size also must be increased.

11

12

Stearns SC. 1992. The evolution of life histories. New York, NY: Oxford University Press. Tieleman IB, Williams JB, Ricklefs RE, Klasing KC. 2005. Constitutive innate immunity is a component of the pace-of-life syndrome in tropical birds. P Roy Soc B Bio 272:1715–20. van de Pol MV, Wright J. 2009. A simple method for distinguishing within- versus between-subject effects using mixed models. Anim Behav 77:753–8. Van Noordwijk AJ, De Jong G. 1986. Acquisition and allocation of resources - their influence on variation in lifehistory tactics. Am Nat 128:137–42. Versteegh MA, Schwabl I, Jaquier S, Tieleman BI. 2012. Do immunological, endocrine and metabolic traits fall on a single Pace-of-Life axis? Covariation and constraints among physiological systems. J Evol Biol 25:1864–76. Wiersma P, Munõz-Garcia A, Walker A, Williams JB. 2007. Tropical birds have a slow pace of life. Proc Natl Acad Sci USA 104:9340–5. Williams GC. 1966. Natural selection, the costs of reproduction, and a refinement of Lack’s principle. Am Nat 100:687–90. Zera AJ, Harshman LG. 2001. The physiology of life history trade-offs in animals. Annu Rev Ecol Syst 32:95–126.


R Core Development Team. 2013. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing. Reznick D, Nunney L, Tessier A. 2000. Big houses, big cars, superfleas and the costs of reproduction. Trends Ecol Evol 15:421–5. Ricklefs RE. 2000. Density dependence, evolutionary optimization, and the diversification of avian life histories. Condor 102:9–22. Ricklefs RE, Wikelski M. 2002. The physiology/life-history nexus. Trends Ecol Evol 17:462–8. Roff DA, Fairbairn DJ. 2007. The evolution of trade-offs: where are we? J Evol Biol 20:433–47. Saether BE. 1988. Pattern of covariation between life-history traits of European birds. Nature 331:616–7. Schmid-Hempel P. 2003. Variation in immune defence as a question of evolutionary ecology. P Roy Soc B Bio 270:357–66. Schmid-Hempel P, Ebert D. 2003. On the evolutionary ecology of specific immune defence. Trends Ecol Evol 18:27–32. Sibly RM, Calow P. 1986. Physiological ecology of animals: An evolutionary approach. Oxford: Blackwell Scientific Publications. Speakman JR, Kro´l E, Johnson MS. 2004. The functional significance of individual variation in basal metabolic rate. Physiol Biochem Zool 77:900–15.


Multiple testing correction in linear mixed models.

PACIC Instrument: disentangling dimensions using published validation models.

Testing Predictive Developmental Hypotheses.

A note on testing families of hypotheses using graphical procedures.

Approximate testing in two-stage nonlinear mixed models.

Testing hypotheses on dioptric power.

Forelimb preferences in human beings and other species: multiple models for testing hypotheses on lateralization.

Predictors for physical activity in adolescent girls using statistical shrinkage techniques for hierarchical longitudinal mixed effects models.

Hypotheses, limits, models and life.

Effective genetic-risk prediction using mixed models.

Testing hypotheses of bat baculum function with 3D models derived from microCT.

Branch: an interactive, web-based tool for testing hypotheses and developing predictive models.

Inversion of hierarchical Bayesian models using Gaussian processes.

Hierarchical control using networks trained with higher-level forward models.

Assessing evidence and testing appropriate hypotheses.

Testing global memory models using ROC curves.

Using Bayes factors for testing hypotheses about intervention effectiveness in addictions research.

What drove reversions to quadrupedality in ornithischian dinosaurs? Testing hypotheses using centre of mass modelling.

Multiple hypotheses testing procedures in clinical trials and genomic studies.

[Use of hierarchical models in nephrology].

Testing evolutionary hypotheses for DNA barcoding failure in willows.

Modeling the oxygen uptake kinetics during exercise testing of patients with chronic obstructive pulmonary diseases using nonlinear mixed models.

Estimating Identification Disclosure Risk Using Mixed Membership Models.

Differential expression analysis for RNAseq using Poisson mixed models.