Evolutionary allometry of horn length in the mammalian family Bovidae reconciled by non-linear regression

Evolutionary allometry of horn length in the mammalian family Bovidae reconciled by non-linear... Abstract I used the equivalent of ANCOVA for non-linear observations to re-examine data for allometric variation in length of the horns in male and female bovid mammals. The data came from an earlier investigation that arrived at several important, but weakly supported, conclusions about the relationship between horn size and body size across species. The new analysis on untransformed observations revealed that patterns of variation in horn length are described better by three-parameter power equations than by two-parameter power functions that were the basis for the original study. The Y-intercepts and allometric exponents estimated in the new analysis are identical for both sexes, but the allometric coefficient is larger for males than for females, thereby accounting for the higher trajectory for the curve fitted to data for males. The fitted equations point to sexual dimorphism in length of the horns that increases in concert with increases in body size. However, the function describing horns in males may have been unduly influenced by statistical outliers, so the pattern of allometric variation in males is not settled. The analytical protocol described here is more versatile than conventional and phylogenetic ANCOVAs on log transformations. analysis of covariance, ANCOVA, sexual dimorphism INTRODUCTION A recent investigation on horns of bovid mammals arrived at several important conclusions concerning patterns of allometric variation in length of the horns and the evolution of sexual dimorphism across species (Tidière et al., 2017). The study was a major undertaking in which linear and quadratic polynomials were fitted by phylogenetic generalized least squares (PGLS) to logarithmic transformations of mean values for horn length and body mass for both males (91 species) and females (54 species). The rectilinear model apparently was the better fit to observations for females (Fig. 1A) but the quadratic model (with a negative coefficient for the quadratic term) was the better fit for males (Fig. 1B). The straight line fitted to transformations for females implied that the original data can be described by a two-parameter power equation Figure 1. View largeDownload slide Logarithmic transformations of mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Bivariate display of observations for females (N = 54 species) without the visual clutter introduced by data for males. The line follows the path of the linear polynomial reported by Tidière et al. (2017), but the function seems not to capture the main pattern in the observations (see also their fig. S6). (B) Bivariate display of observations for males (N = 91 species) without the visual clutter introduced by data for females. The curve follows the path of the quadratic polynomial reported by Tidière et al. (2017). (C) Bivariate display of observations for both males and females. The greatest difference between the sexes is for species weighing approximately 100 kg (Tidière et al., 2017). Figure 1. View largeDownload slide Logarithmic transformations of mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Bivariate display of observations for females (N = 54 species) without the visual clutter introduced by data for males. The line follows the path of the linear polynomial reported by Tidière et al. (2017), but the function seems not to capture the main pattern in the observations (see also their fig. S6). (B) Bivariate display of observations for males (N = 91 species) without the visual clutter introduced by data for females. The curve follows the path of the quadratic polynomial reported by Tidière et al. (2017). (C) Bivariate display of observations for both males and females. The greatest difference between the sexes is for species weighing approximately 100 kg (Tidière et al., 2017).  Y=a*Xb with a fixed exponent (see their equations 1 and 2). The curvilinear model for males, however, was interpreted to mean that the exponent in a two-parameter allometric equation describing untransformed data is a declining function of body size (i.e. ‘complex allometry’ sensu Strauss, 1993). When the functions fitted to logarithms were displayed graphically, sexual dimorphism appeared to be most pronounced for species of intermediate size (Fig. 1C). The latter finding suggested, in turn, that the size of horns in males of large species is limited by unknown constraints or, alternatively, that their horns are not targets of sexual selection (Tidière et al., 2017). Two procedural issues need to be revisited, however, and both potentially impact the most important conclusions from the investigation. First, the curvilinear path followed by the function fitted to transformations of data for males actually indicates that a two-parameter power equation is not a suitable statistical model for describing pattern in the original data and that a model with different functional form should have been pursued (e.g. Reeve, 1940; Sprent, 1972; Stumpf & Porter, 2012). Second, the statistical models fitted to logarithms were interpreted in the context of the original measurements but without the benefit of graphical validation of the models on the arithmetic scale. It always is best to validate and interpret analyses on the scale of measurement (Finney, 1989), because logarithmic transformation alters bivariate distributions in ways that may conceal flaws in the original data (Packard, 2011, 2014, 2017a). Here I use a relatively new statistical approach to allometry to re-examine the same data that were studied by Tidière et al. (2017). Accordingly, means for body mass and horn length for male and female bovids were downloaded from the online supplement for their report and submitted to a new analysis. The new treatment identifies patterns in distributions for both males and females that went undetected in the original investigation. Several of the interpretations from the earlier investigation are thereby affected. MATERIAL AND METHODS The data for bovids appear at first to be well-suited for study by either standard (e.g. Cochran, 1957) or phylogenetic ANCOVA (e.g. Fuentes-G et al., 2016), with the observations being grouped by sex in a one-way analysis. However, both of these protocols require that logarithmic transformations of the predictor and response variables follow a linear path in bivariate display (Cochran, 1957; Fuentes-G et al., 2016). As noted already, the requirement for linearity is not met by transformed observations for male bovids (Fig. 1B), and a question might even be raised about the linearity of transformed observations for females (Fig. 1A). Consequently, neither type of ANCOVA is appropriate for examining this particular data set. I therefore examined the data by using a procedure that combines features of the ANOVA with non-linear regression (Lolli et al., 2017; Packard, 2018). The overall approach is not new (e.g. Nevill & Holder, 1994), but was until recently constrained to fitting models that assume normal, homoscedastic error (the default in most routines for non-linear regression). However, the procedure now can be applied to untransformed data exhibiting random error that is normal and homoscedastic, normal and heteroscedastic, or lognormal and heteroscedastic on the original scale (Table 1). Because the procedure can be applied to observations for X and Y that are linearized by transformation and to those that are not (i.e. ‘complex allometry’), the protocol is extraordinarily versatile. Although the protocol cannot at present include an adjustment for phylogenetic affinity, this is a minor concern when ANCOVA on transformations would yield a poor fit and questionable conclusions (as in the case of the male bovids). [Phylogenetic methods were developed for use with the loglinear model that has been at the heart of allometric research since the time of Huxley (1932). Thus, corrections for phylogeny are based on the assumption that the original, untransformed observations for X and Y follow the path of a two-parameter power equation with lognormal, heteroscedastic error (Packard, 2017b). However, PGLS weights residuals (e.g. Grafen, 1989, 1992; Revell, 2010; Smaers & Rohlf, 2016), and phylogenetic corrections actually can be made to any model that is linear in its parameters (such as a quadratic polynomial) and that has residuals expressed as logarithms. If the observations for X and Y are not loglinear, however, the corrections are made to a model having the wrong functional form on the arithmetic scale (e.g. Reeve, 1940; Sprent, 1972; Stumpf & Porter, 2012).] Table 1. Different forms for random error in regression models Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  β is a vector of parameters in a two- or three-parameter power function, f, relating the response variable, Y, to the measure of size, X. Error for all the models is assumed to be normally distributed, to sum to zero and to have the specified variance. Variance for models with multiplicative, lognormal, heteroscedastic error is based on residuals expressed on the logarithmic scale (see http://support.sas.com/documentation/cdl/en/etsug/60372/HTML/default/viewer.htm#etsug_model_sect045.htm). View Large Table 1. Different forms for random error in regression models Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  β is a vector of parameters in a two- or three-parameter power function, f, relating the response variable, Y, to the measure of size, X. Error for all the models is assumed to be normally distributed, to sum to zero and to have the specified variance. Variance for models with multiplicative, lognormal, heteroscedastic error is based on residuals expressed on the logarithmic scale (see http://support.sas.com/documentation/cdl/en/etsug/60372/HTML/default/viewer.htm#etsug_model_sect045.htm). View Large Take the following equation as a starting point:  Y=(a*X(b+d*X2)+Y0)*exp(c*X2) where the first parenthetical expression represents a three-parameter power equation and the second parenthetical expression represents the overall influence of the group to which a subject is assigned. Y is the response variable (i.e. horn length), X is the predictor variable (i.e. body mass) and X2 is a categorical variable designating the group (e.g. 0 for females, 1 for males). The other terms in the equation are fitted constants (or parameters), with a being the allometric coefficient, b the allometric exponent when X2 is 0, c a measure of the overall influence of group, d the added influence imparted to the exponent when X2 is 1, and Y0 the intercept with the Y axis. When the equation is evaluated at both levels for X2 it yields a pair of three-parameter equations, one for females (X2 = 0) and another for males (X2 = 1). Note that X2 affects the exponent in the equation (via the parameter d) as well as the allometric coefficient and the intercept, both of which are located inside the first set of parentheses where they are conditioned by a multiplier which includes X2. The basic model and several modifications thereof (e.g. models without one or another of the possible parameters, or with parentheses located differently to change the sequence of steps in computation) were used to explore the data for bovids using the Model Procedure in SAS 9.3. Models with different forms for random error were included in the mix (Table 1). Code for the most relevant models is available in Supporting Information Table S1. Akaike’s Information Criterion (AIC) was used to distinguish among good models, plausible models, and weakly supported or implausible models (Burnham & Anderson, 2002; Anderson, 2008; Burnham, Anderson & Huyvaert, 2011). RESULTS Preliminary analyses indicated that models with normally distributed error (both homoscedastic and heteroscedastic) were poor choices for describing the data set, so I concentrated on fitting models with lognormal error (Table S1). I also determined that there was little support for two-parameter power equations or for two of the three-parameter equations in the pool of candidate models (∆AIC ≥ 9.6). Four models survived this screening (Table 2). Table 2. Equations for curvilinear Y vs. X together with fitted functions describing patterns of allometric variation in horn length of bovids General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  View Large Table 2. Equations for curvilinear Y vs. X together with fitted functions describing patterns of allometric variation in horn length of bovids General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  View Large Two of the models in the reduced pool of candidates were good (equivalent) fits, and two were plausible alternatives to the best fit (Table 2). However, the term d did not enter at a significant level into models 1 and 2, and the term Y0 did not enter at a significant level into model 3 (Table 3). When d was removed from models 1 and 2, the models approached model 4 (which had the lowest value for AIC). When Y0 was removed from model 3, it became an unsupported two-parameter equation. All the estimated parameters entered significantly into model 4 (Table 3). Moreover, the Breusch–Pagan test for homoscedasticity was acceptable for model 4 (P = 0.39), and so, too, were standardized residuals (Fig. 2A). Residuals also passed the Shapiro–Wilk test for normality (P = 0.11). Thus, the best model in the pool of candidate models is number 4 (Table 2). Table 3. Estimates for parameters in models 1–4 and their levels of significance   Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034    Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034  View Large Table 3. Estimates for parameters in models 1–4 and their levels of significance   Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034    Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034  View Large Figure 2. View largeDownload slide Mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Standardized residuals vs. predicted values from the fitting of model 4 (Table 2). Residuals are in logarithmic form because of the assumption regarding lognormal error (Table 1). The computational algorithm has reversed the sign for each residual. (B) Bivariate display comparing patterns of variation in horn length for males and females. Fitted functions are based on model 4 (Table 2). (C) Observations for females are displayed without the visual clutter introduced by data for males. The path for the curve is exactly what is expected for an equation that predicts geometric means for the response variable (Miller, 1984). Observations for 37 species that lack horns are included here for perspective, but they were not included in the statistical analysis. (D) Observations for males are displayed without the visual clutter introduced by data for females. The distribution for males is slightly unusual because of the numerous small species that have long horns. Nevertheless, the function shown here is based on the best statistical model in the pool of candidate models. Figure 2. View largeDownload slide Mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Standardized residuals vs. predicted values from the fitting of model 4 (Table 2). Residuals are in logarithmic form because of the assumption regarding lognormal error (Table 1). The computational algorithm has reversed the sign for each residual. (B) Bivariate display comparing patterns of variation in horn length for males and females. Fitted functions are based on model 4 (Table 2). (C) Observations for females are displayed without the visual clutter introduced by data for males. The path for the curve is exactly what is expected for an equation that predicts geometric means for the response variable (Miller, 1984). Observations for 37 species that lack horns are included here for perspective, but they were not included in the statistical analysis. (D) Observations for males are displayed without the visual clutter introduced by data for females. The distribution for males is slightly unusual because of the numerous small species that have long horns. Nevertheless, the function shown here is based on the best statistical model in the pool of candidate models. The best model is something of an abstraction at this point because it provides information for two different power equations with lognormal, heteroscedastic error. I therefore evaluated the model for both males and females and plotted the functions on a graph together with the original data (Fig. 2B). The function for females is a good visual fit that follows the expected path for a statistical model predicting geometric means (Fig. 2C). The curve for males has a higher trajectory than does the curve for females (Fig. 2B), but the distribution of data for males is slightly peculiar and raises a question about the suitability of the data for use in an allometric analysis (Fig. 2D). Nevertheless, these power functions are the best of the various functions in the pool of candidate models (Tables 2, S1), and they are based on the same data that were examined by Tidière et al. (2017). The equations for males and females differ only in the value for the allometric coefficient, a; the equations share the same values for the allometric exponent and for the intercept with the Y axis (Fig. 2C, D). DISCUSSION The findings presented here differ in many respects from those reported by Tidière et al. (2017). Here I enumerate the differences and offer a new perspective on the evolutionary allometry of horns in bovids. First, the patterns of allometric variation in horn length are described best by three-parameter power equations (Table 2) and not by two-parameter power equations that were the basis for the investigation by Tidière et al. (2017, see their equation 1, which established the starting point and frame of reference for the study). Simply stated, an explicit intercept is needed in both the allometric equations to capture the maximum amount of information that is available in the data (Burnham & Anderson, 2002; Anderson, 2008). The presence of an intercept in the function for males is what caused logarithmic transformations to be curvilinear in bivariate display (Fig. 1B) (Sartori & Ball, 2009; Packard, 2017b, 2017c). A requirement for an intercept may also explain the curvature that appears in transformed observations for females (Fig. 1A). Second, Y-intercepts in the power equations described here are negative, thereby introducing a complication for students of allometry who believe that allometric equations with negative intercepts defy biologically realistic interpretation (e.g. Klingenberg, 1998; Lemaître et al., 2015). However, the complication results from undue focus on the Y-intercept (‘animals with horns having negative lengths’) when focus should be on the X-intercept (Bales, 1996; Sartori & Ball, 2009; Packard, 2015, 2016, 2017b), which identifies body mass when length of horns is zero. A three-parameter equation indicates that the probability of possessing horns diminishes with body size and that horns are unlikely to be present in species weighing less than a threshold for body mass. The X-intercept (i.e. the threshold) is approximately 1.22 kg for males and 2.90 kg for females. The smallest males with horns have a mass of 2 kg (Neotragus pygmaeus), and the smallest females bearing horns have a mass of 5.1 kg (Cephalophus monticola). All male bovids possess horns, but all of them also exceed the putative threshold; consequently, 1.22 kg may actually be the lower limit for size of males in any species of bovid. In contrast, females in three species weigh on average less than the threshold of 2.90 kg, and none of them possesses horns. Indeed, horns are absent from females of 34 other species weighing more than 2.90 kg (Fig. 2C), so attaining the threshold size is no guarantee that natural selection will favour females that carry horns. The alternative to three-parameter power equations, of course, would be two-parameter power equations that of necessity pass through the origin, thereby indicating that all adult bovids, irrespective of their sex or body mass, should possess horns. Third, exponents in the allometric equations for males and females are identical, and they remain unchanged over the full range in body masses (Fig. 2C, D). The perception of Tidière et al. (2017) that the allometric exponent is a variable in the equation fitted to data for males was caused by performing analyses in the logarithmic domain instead of on the arithmetic scale (Packard, 2017c). Fourth, the pattern of variation in horn length differs between males and females (Fig. 2B), but this difference is the result of effects of sex on the allometric coefficient (a). Given the emphasis that is usually given in allometry to the influence of the allometric exponent, it probably will come as a surprise to many readers that the pronounced difference in trajectories for curves in Figure 2B is due to differences in the allometric coefficient and not the allometric exponent. Such a difference would probably go undetected in studies treating data for males and females separately because the allometric coefficient is commonly ignored. Finally, the equations based on model 4 point to sexual dimorphism in horns that becomes more pronounced as species become larger (Fig. 2B). This finding stands in contrast to that reported by Tidière et al. (2017), who contend that species of intermediate size are more dimorphic than species from the extremes of the size distribution. Of course, a study of allometry is no better than the data on which it is based, and the present treatment is no exception. Although observations for horn length in females appear to be quite good (Fig. 2C), those for horns in males are problematic (Fig. 2D). The observations for males seem to comprise two different groups: a group of 85 species weighing less than 300 kg and a second group of six species weighing in excess of ~600 kg. Data for horns in the former group form an off-vertical column at the left of the graph, whereas observations for the latter group lie well to the right and substantially below the fitted curve. Thus, observations for males may not have a unitary pattern amenable to allometric analysis, in which case it would be premature to discuss a pattern of allometric variation that may be spurious. On the other hand, the odd distribution may be nothing more than a trivial anomaly caused by a gap in the sample, in which case the fitted curve may provide a reasonably accurate description for underlying pattern in the data. The question of the bivariate distribution needs to be resolved. However, neither of the aforementioned alternatives lends support to the claim that sexual dimorphism in length of the horns is greatest for species of bovids in the mid-range for body size. SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher’s website: Table S1. SAS code for implementing ANCOVA when the relationship between Y and X is curvilinear. ACKNOWLEDGEMENTS I thank M. Tidière, A. Nevill and an anonymous referee for thoughtful comments and criticisms that helped me to sharpen the focus of the manuscript. This investigation was funded out of pocket and received no financial support from any governmental or institutional source. REFERENCES Anderson DR. 2008. Model based inference in the life sciences: a primer on evidence . New York: Springer. Google Scholar CrossRef Search ADS   Bales GS. 1996. Heterochrony in brontothere horn evolution: allometric interpretations and the effect of life history scaling. Paleobiology   22: 481– 495. Google Scholar CrossRef Search ADS   Burnham KP, Anderson DR. 2002. Model selection and multimodel inference, 2nd edn . New York: Springer. Burnham KP, Anderson DR, Huyvaert KP. 2011. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behavioral Ecology and Sociobiology   65: 23– 35. Google Scholar CrossRef Search ADS   Cochran WG. 1957. Analysis of covariance: its nature and uses. Biometrics   13: 261– 281. Google Scholar CrossRef Search ADS   Finney DJ. 1989. Was this in your statistics textbook? V. Transformation of data. Experimental Agriculture   25: 165– 175. Google Scholar CrossRef Search ADS   Fuentes-G JA, Housworth EA, Weber A, Martins EP. 2016. Phylogenetic ANCOVA: estimating changes in evolutionary rates as well as relationships between traits. The American Naturalist   188: 615– 627. Google Scholar CrossRef Search ADS PubMed  Grafen A. 1989. The phylogenetic regression. Philosophical Transactions of the Royal Society of London   326: 119– 157. Google Scholar CrossRef Search ADS PubMed  Grafen A. 1992. The uniqueness of the phylogenetic regression. Journal of Theoretical Biology   156: 405– 423. Google Scholar CrossRef Search ADS   Huxley JS. 1932. Problems of relative growth . London: Methuen & Co. Klingenberg CP. 1998. Heterochrony and allometry: the analysis of evolutionary change in ontogeny. Biological Reviews   73: 79– 123. Google Scholar CrossRef Search ADS PubMed  Lemaître JF, Vanpé C, Plard F, Pélabon C, Gaillard JM. 2015. Response to Packard: make sure we do not throw out the biological baby with the statistical bath water when performing allometric analyses. Biology Letters   11: 20150144. Google Scholar CrossRef Search ADS PubMed  Lolli L, Batterham AM, Kratochvíl L, Flegr J, Weston KL, Atkinson G. 2017. A comprehensive allometric analysis of 2nd digit length to 4th digit length in humans. Proceedings of the Royal Society B   284: 20170356. Google Scholar CrossRef Search ADS PubMed  Miller DR. 1984. Reducing transformation bias in curve fitting. American Statistician   38: 124– 126. Nevill AM, Holder RL. 1994. Modelling maximum oxygen uptake—a case study in non-linear regression model formulation and comparison. Journal of the Royal Statistical Society C   43: 653– 666. Packard GC. 2011. Unanticipated consequences of logarithmic transformation in bivariate allometry. Journal of Comparative Physiology B   181: 841– 849. Google Scholar CrossRef Search ADS   Packard GC. 2014. On the use of log-transformation versus nonlinear regression for analyzing biological power laws. Biological Journal of the Linnean Society   113: 1167– 1178. Google Scholar CrossRef Search ADS   Packard GC. 2015. Allometric variation in the antlers of cervids: a comment on Lemaître et al. Biology Letters   11: 20140923. Google Scholar CrossRef Search ADS PubMed  Packard GC. 2016. Relative growth by the elongated jaws of gars: a perspective on polyphasic loglinear allometry. Journal of Experimental Zoology B   326: 168– 175. Google Scholar CrossRef Search ADS   Packard GC. 2017a. The essential role for graphs in allometric analysis. Biological Journal of the Linnean Society   120: 468– 473. Packard GC. 2017b. Misconceptions about logarithmic transformation and the traditional allometric method. Zoology   123: 115– 120. Google Scholar CrossRef Search ADS   Packard GC. 2017c. Why allometric variation in mammalian metabolism is curvilinear on the logarithmic scale. Journal of Experimental Zoology A   327: 537– 541. Google Scholar CrossRef Search ADS   Packard GC. 2018. A new research paradigm for bivariate allometry: combining ANOVA and nonlinear regression. Journal of Experimental Biology   221: jeb177519. Reeve ECR. 1940. Relative growth in the snout of anteaters. A study in the application of quantitative methods to systematics. Proceedings of the Zoological Society of London A   110: 47– 80. Revell LJ. 2010. Phylogenetic signal and linear regression on species data. Methods in Ecology and Evolution   1: 319– 329. Google Scholar CrossRef Search ADS   Sartori AF, Ball AD. 2009. Morphology and postlarval development of the ligament of Thracia phaseolina (Bivalvia: Thraciidae), with a discussion of model choice in allometric studies. Journal of Molluscan Studies   75: 295– 304. Google Scholar CrossRef Search ADS   Smaers JB, Rohlf FJ. 2016. Testing species’ deviation from allometric predictions using the phylogenetic regression. Evolution   70: 1145– 1149. Google Scholar CrossRef Search ADS PubMed  Sprent P. 1972. The mathematics of size and shape. Biometrics   28: 23– 37. Google Scholar CrossRef Search ADS PubMed  Strauss RE. 1993. The study of allometry since Huxley. In: Huxley JS. Problems of relative growth  [new edition]. Baltimore: Johns Hopkins University Press, xlvii– lxxv. (http://www.faculty.biol.ttu.edu/Strauss/Pubs/Papers/1993StraussHuxley.pdf) Stumpf MPH, Porter MA. 2012. Critical truths about power laws. Science   335: 665– 666. Google Scholar CrossRef Search ADS PubMed  Tidière M, Lemaître JF, Pélabon C, Gimenez O, Gaillard JM. 2017. Evolutionary allometry reveals a shift in selection pressure on male horn size. Journal of Evolutionary Biology   30: 1826– 1835. Google Scholar CrossRef Search ADS PubMed  © 2018 The Linnean Society of London, Biological Journal of the Linnean Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Biological Journal of the Linnean Society Oxford University Press

Evolutionary allometry of horn length in the mammalian family Bovidae reconciled by non-linear regression

Loading next page...
 
/lp/ou_press/evolutionary-allometry-of-horn-length-in-the-mammalian-family-bovidae-Uzoxp0jKXf
Publisher
The Linnean Society of London
Copyright
© 2018 The Linnean Society of London, Biological Journal of the Linnean Society
ISSN
0024-4066
eISSN
1095-8312
D.O.I.
10.1093/biolinnean/bly052
Publisher site
See Article on Publisher Site

Abstract

Abstract I used the equivalent of ANCOVA for non-linear observations to re-examine data for allometric variation in length of the horns in male and female bovid mammals. The data came from an earlier investigation that arrived at several important, but weakly supported, conclusions about the relationship between horn size and body size across species. The new analysis on untransformed observations revealed that patterns of variation in horn length are described better by three-parameter power equations than by two-parameter power functions that were the basis for the original study. The Y-intercepts and allometric exponents estimated in the new analysis are identical for both sexes, but the allometric coefficient is larger for males than for females, thereby accounting for the higher trajectory for the curve fitted to data for males. The fitted equations point to sexual dimorphism in length of the horns that increases in concert with increases in body size. However, the function describing horns in males may have been unduly influenced by statistical outliers, so the pattern of allometric variation in males is not settled. The analytical protocol described here is more versatile than conventional and phylogenetic ANCOVAs on log transformations. analysis of covariance, ANCOVA, sexual dimorphism INTRODUCTION A recent investigation on horns of bovid mammals arrived at several important conclusions concerning patterns of allometric variation in length of the horns and the evolution of sexual dimorphism across species (Tidière et al., 2017). The study was a major undertaking in which linear and quadratic polynomials were fitted by phylogenetic generalized least squares (PGLS) to logarithmic transformations of mean values for horn length and body mass for both males (91 species) and females (54 species). The rectilinear model apparently was the better fit to observations for females (Fig. 1A) but the quadratic model (with a negative coefficient for the quadratic term) was the better fit for males (Fig. 1B). The straight line fitted to transformations for females implied that the original data can be described by a two-parameter power equation Figure 1. View largeDownload slide Logarithmic transformations of mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Bivariate display of observations for females (N = 54 species) without the visual clutter introduced by data for males. The line follows the path of the linear polynomial reported by Tidière et al. (2017), but the function seems not to capture the main pattern in the observations (see also their fig. S6). (B) Bivariate display of observations for males (N = 91 species) without the visual clutter introduced by data for females. The curve follows the path of the quadratic polynomial reported by Tidière et al. (2017). (C) Bivariate display of observations for both males and females. The greatest difference between the sexes is for species weighing approximately 100 kg (Tidière et al., 2017). Figure 1. View largeDownload slide Logarithmic transformations of mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Bivariate display of observations for females (N = 54 species) without the visual clutter introduced by data for males. The line follows the path of the linear polynomial reported by Tidière et al. (2017), but the function seems not to capture the main pattern in the observations (see also their fig. S6). (B) Bivariate display of observations for males (N = 91 species) without the visual clutter introduced by data for females. The curve follows the path of the quadratic polynomial reported by Tidière et al. (2017). (C) Bivariate display of observations for both males and females. The greatest difference between the sexes is for species weighing approximately 100 kg (Tidière et al., 2017).  Y=a*Xb with a fixed exponent (see their equations 1 and 2). The curvilinear model for males, however, was interpreted to mean that the exponent in a two-parameter allometric equation describing untransformed data is a declining function of body size (i.e. ‘complex allometry’ sensu Strauss, 1993). When the functions fitted to logarithms were displayed graphically, sexual dimorphism appeared to be most pronounced for species of intermediate size (Fig. 1C). The latter finding suggested, in turn, that the size of horns in males of large species is limited by unknown constraints or, alternatively, that their horns are not targets of sexual selection (Tidière et al., 2017). Two procedural issues need to be revisited, however, and both potentially impact the most important conclusions from the investigation. First, the curvilinear path followed by the function fitted to transformations of data for males actually indicates that a two-parameter power equation is not a suitable statistical model for describing pattern in the original data and that a model with different functional form should have been pursued (e.g. Reeve, 1940; Sprent, 1972; Stumpf & Porter, 2012). Second, the statistical models fitted to logarithms were interpreted in the context of the original measurements but without the benefit of graphical validation of the models on the arithmetic scale. It always is best to validate and interpret analyses on the scale of measurement (Finney, 1989), because logarithmic transformation alters bivariate distributions in ways that may conceal flaws in the original data (Packard, 2011, 2014, 2017a). Here I use a relatively new statistical approach to allometry to re-examine the same data that were studied by Tidière et al. (2017). Accordingly, means for body mass and horn length for male and female bovids were downloaded from the online supplement for their report and submitted to a new analysis. The new treatment identifies patterns in distributions for both males and females that went undetected in the original investigation. Several of the interpretations from the earlier investigation are thereby affected. MATERIAL AND METHODS The data for bovids appear at first to be well-suited for study by either standard (e.g. Cochran, 1957) or phylogenetic ANCOVA (e.g. Fuentes-G et al., 2016), with the observations being grouped by sex in a one-way analysis. However, both of these protocols require that logarithmic transformations of the predictor and response variables follow a linear path in bivariate display (Cochran, 1957; Fuentes-G et al., 2016). As noted already, the requirement for linearity is not met by transformed observations for male bovids (Fig. 1B), and a question might even be raised about the linearity of transformed observations for females (Fig. 1A). Consequently, neither type of ANCOVA is appropriate for examining this particular data set. I therefore examined the data by using a procedure that combines features of the ANOVA with non-linear regression (Lolli et al., 2017; Packard, 2018). The overall approach is not new (e.g. Nevill & Holder, 1994), but was until recently constrained to fitting models that assume normal, homoscedastic error (the default in most routines for non-linear regression). However, the procedure now can be applied to untransformed data exhibiting random error that is normal and homoscedastic, normal and heteroscedastic, or lognormal and heteroscedastic on the original scale (Table 1). Because the procedure can be applied to observations for X and Y that are linearized by transformation and to those that are not (i.e. ‘complex allometry’), the protocol is extraordinarily versatile. Although the protocol cannot at present include an adjustment for phylogenetic affinity, this is a minor concern when ANCOVA on transformations would yield a poor fit and questionable conclusions (as in the case of the male bovids). [Phylogenetic methods were developed for use with the loglinear model that has been at the heart of allometric research since the time of Huxley (1932). Thus, corrections for phylogeny are based on the assumption that the original, untransformed observations for X and Y follow the path of a two-parameter power equation with lognormal, heteroscedastic error (Packard, 2017b). However, PGLS weights residuals (e.g. Grafen, 1989, 1992; Revell, 2010; Smaers & Rohlf, 2016), and phylogenetic corrections actually can be made to any model that is linear in its parameters (such as a quadratic polynomial) and that has residuals expressed as logarithms. If the observations for X and Y are not loglinear, however, the corrections are made to a model having the wrong functional form on the arithmetic scale (e.g. Reeve, 1940; Sprent, 1972; Stumpf & Porter, 2012).] Table 1. Different forms for random error in regression models Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  β is a vector of parameters in a two- or three-parameter power function, f, relating the response variable, Y, to the measure of size, X. Error for all the models is assumed to be normally distributed, to sum to zero and to have the specified variance. Variance for models with multiplicative, lognormal, heteroscedastic error is based on residuals expressed on the logarithmic scale (see http://support.sas.com/documentation/cdl/en/etsug/60372/HTML/default/viewer.htm#etsug_model_sect045.htm). View Large Table 1. Different forms for random error in regression models Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  Type of error  Model form  Additive, normal, homoscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, σ2)  Additive, normal, heteroscedastic error  Yi = f (Xi, β) + εi εi ~ N(0, (σ2 * Xi2θ))  Multiplicative, lognormal, heteroscedastic error  Yi = f (Xi, β) * exp(εi) εi ~ N(0, σ2)  β is a vector of parameters in a two- or three-parameter power function, f, relating the response variable, Y, to the measure of size, X. Error for all the models is assumed to be normally distributed, to sum to zero and to have the specified variance. Variance for models with multiplicative, lognormal, heteroscedastic error is based on residuals expressed on the logarithmic scale (see http://support.sas.com/documentation/cdl/en/etsug/60372/HTML/default/viewer.htm#etsug_model_sect045.htm). View Large Take the following equation as a starting point:  Y=(a*X(b+d*X2)+Y0)*exp(c*X2) where the first parenthetical expression represents a three-parameter power equation and the second parenthetical expression represents the overall influence of the group to which a subject is assigned. Y is the response variable (i.e. horn length), X is the predictor variable (i.e. body mass) and X2 is a categorical variable designating the group (e.g. 0 for females, 1 for males). The other terms in the equation are fitted constants (or parameters), with a being the allometric coefficient, b the allometric exponent when X2 is 0, c a measure of the overall influence of group, d the added influence imparted to the exponent when X2 is 1, and Y0 the intercept with the Y axis. When the equation is evaluated at both levels for X2 it yields a pair of three-parameter equations, one for females (X2 = 0) and another for males (X2 = 1). Note that X2 affects the exponent in the equation (via the parameter d) as well as the allometric coefficient and the intercept, both of which are located inside the first set of parentheses where they are conditioned by a multiplier which includes X2. The basic model and several modifications thereof (e.g. models without one or another of the possible parameters, or with parentheses located differently to change the sequence of steps in computation) were used to explore the data for bovids using the Model Procedure in SAS 9.3. Models with different forms for random error were included in the mix (Table 1). Code for the most relevant models is available in Supporting Information Table S1. Akaike’s Information Criterion (AIC) was used to distinguish among good models, plausible models, and weakly supported or implausible models (Burnham & Anderson, 2002; Anderson, 2008; Burnham, Anderson & Huyvaert, 2011). RESULTS Preliminary analyses indicated that models with normally distributed error (both homoscedastic and heteroscedastic) were poor choices for describing the data set, so I concentrated on fitting models with lognormal error (Table S1). I also determined that there was little support for two-parameter power equations or for two of the three-parameter equations in the pool of candidate models (∆AIC ≥ 9.6). Four models survived this screening (Table 2). Table 2. Equations for curvilinear Y vs. X together with fitted functions describing patterns of allometric variation in horn length of bovids General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  View Large Table 2. Equations for curvilinear Y vs. X together with fitted functions describing patterns of allometric variation in horn length of bovids General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  General form for model  Fitted equation  Log likelihood  AIC  ∆ AIC  Model 1: Y = (a * X(b+d*X2) + Y0) * exp(c*X2)  Ŷ = (1.71 * X(0.65 − 0.19*X2) − 1.81) * exp(1.37*X2)  −597.6  1207.2  3.5  Model 2: Y = (a * X(b+d*X2) * exp(c*X2)) + Y0  Ŷ = (3.62 * X(0.52− 0.06*X2) * exp(0.59*X2)) – 6.76  −596.1  1204.1  0.4  Model 3: Y = (a * Xb + Y0) * exp(c*X2)  Ŷ = (2.83 * X0.53 – 2.87) * exp(0.51*X2)  −599.4  1208.8  5.1  Model 4: Y = (a * Xb) * exp(c*X2) + Y0  Ŷ = (4.36 * X0.47) * exp(0.41*X2) – 7.19  −596.9  1203.7  0  View Large Two of the models in the reduced pool of candidates were good (equivalent) fits, and two were plausible alternatives to the best fit (Table 2). However, the term d did not enter at a significant level into models 1 and 2, and the term Y0 did not enter at a significant level into model 3 (Table 3). When d was removed from models 1 and 2, the models approached model 4 (which had the lowest value for AIC). When Y0 was removed from model 3, it became an unsupported two-parameter equation. All the estimated parameters entered significantly into model 4 (Table 3). Moreover, the Breusch–Pagan test for homoscedasticity was acceptable for model 4 (P = 0.39), and so, too, were standardized residuals (Fig. 2A). Residuals also passed the Shapiro–Wilk test for normality (P = 0.11). Thus, the best model in the pool of candidate models is number 4 (Table 2). Table 3. Estimates for parameters in models 1–4 and their levels of significance   Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034    Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034  View Large Table 3. Estimates for parameters in models 1–4 and their levels of significance   Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034    Parameter  Estimate  SE  t-value  Probability  Model 1  a  1.71  0.58  2.94  0.004    b  0.65  0.08  7.82  < 0.001    c  1.37  0.51  2.72  0.007    d  −0.19  0.11  −1.72  0.087    Y0  −1.81  0.83  −2.17  0.032  Model 2  a  3.62  1.75  2.07  0.040    b  0.52  0.10  5.32  < 0.001    c  0.59  0.24  2.48  0.014    d  −0.06  0.06  −0.88  0.382    Y0  −6.76  3.12  −2.14  0.034  Model 3  a  2.83  0.94  3.01  0.003    b  0.53  0.07  7.96  < 0.001    c  0.51  0.10  4.90  < 0.001    Y0  −2.87  1.67  −1.72  0.088  Model 4  a  4.36  1.76  2.48  0.014    b  0.47  0.07  6.82  < 0.001    c  0.41  0.07  5.62  < 0.001    Y0  −7.19  3.36  −2.14  0.034  View Large Figure 2. View largeDownload slide Mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Standardized residuals vs. predicted values from the fitting of model 4 (Table 2). Residuals are in logarithmic form because of the assumption regarding lognormal error (Table 1). The computational algorithm has reversed the sign for each residual. (B) Bivariate display comparing patterns of variation in horn length for males and females. Fitted functions are based on model 4 (Table 2). (C) Observations for females are displayed without the visual clutter introduced by data for males. The path for the curve is exactly what is expected for an equation that predicts geometric means for the response variable (Miller, 1984). Observations for 37 species that lack horns are included here for perspective, but they were not included in the statistical analysis. (D) Observations for males are displayed without the visual clutter introduced by data for females. The distribution for males is slightly unusual because of the numerous small species that have long horns. Nevertheless, the function shown here is based on the best statistical model in the pool of candidate models. Figure 2. View largeDownload slide Mean horn length (cm) and mean body mass (kg) in bovid mammals. (A) Standardized residuals vs. predicted values from the fitting of model 4 (Table 2). Residuals are in logarithmic form because of the assumption regarding lognormal error (Table 1). The computational algorithm has reversed the sign for each residual. (B) Bivariate display comparing patterns of variation in horn length for males and females. Fitted functions are based on model 4 (Table 2). (C) Observations for females are displayed without the visual clutter introduced by data for males. The path for the curve is exactly what is expected for an equation that predicts geometric means for the response variable (Miller, 1984). Observations for 37 species that lack horns are included here for perspective, but they were not included in the statistical analysis. (D) Observations for males are displayed without the visual clutter introduced by data for females. The distribution for males is slightly unusual because of the numerous small species that have long horns. Nevertheless, the function shown here is based on the best statistical model in the pool of candidate models. The best model is something of an abstraction at this point because it provides information for two different power equations with lognormal, heteroscedastic error. I therefore evaluated the model for both males and females and plotted the functions on a graph together with the original data (Fig. 2B). The function for females is a good visual fit that follows the expected path for a statistical model predicting geometric means (Fig. 2C). The curve for males has a higher trajectory than does the curve for females (Fig. 2B), but the distribution of data for males is slightly peculiar and raises a question about the suitability of the data for use in an allometric analysis (Fig. 2D). Nevertheless, these power functions are the best of the various functions in the pool of candidate models (Tables 2, S1), and they are based on the same data that were examined by Tidière et al. (2017). The equations for males and females differ only in the value for the allometric coefficient, a; the equations share the same values for the allometric exponent and for the intercept with the Y axis (Fig. 2C, D). DISCUSSION The findings presented here differ in many respects from those reported by Tidière et al. (2017). Here I enumerate the differences and offer a new perspective on the evolutionary allometry of horns in bovids. First, the patterns of allometric variation in horn length are described best by three-parameter power equations (Table 2) and not by two-parameter power equations that were the basis for the investigation by Tidière et al. (2017, see their equation 1, which established the starting point and frame of reference for the study). Simply stated, an explicit intercept is needed in both the allometric equations to capture the maximum amount of information that is available in the data (Burnham & Anderson, 2002; Anderson, 2008). The presence of an intercept in the function for males is what caused logarithmic transformations to be curvilinear in bivariate display (Fig. 1B) (Sartori & Ball, 2009; Packard, 2017b, 2017c). A requirement for an intercept may also explain the curvature that appears in transformed observations for females (Fig. 1A). Second, Y-intercepts in the power equations described here are negative, thereby introducing a complication for students of allometry who believe that allometric equations with negative intercepts defy biologically realistic interpretation (e.g. Klingenberg, 1998; Lemaître et al., 2015). However, the complication results from undue focus on the Y-intercept (‘animals with horns having negative lengths’) when focus should be on the X-intercept (Bales, 1996; Sartori & Ball, 2009; Packard, 2015, 2016, 2017b), which identifies body mass when length of horns is zero. A three-parameter equation indicates that the probability of possessing horns diminishes with body size and that horns are unlikely to be present in species weighing less than a threshold for body mass. The X-intercept (i.e. the threshold) is approximately 1.22 kg for males and 2.90 kg for females. The smallest males with horns have a mass of 2 kg (Neotragus pygmaeus), and the smallest females bearing horns have a mass of 5.1 kg (Cephalophus monticola). All male bovids possess horns, but all of them also exceed the putative threshold; consequently, 1.22 kg may actually be the lower limit for size of males in any species of bovid. In contrast, females in three species weigh on average less than the threshold of 2.90 kg, and none of them possesses horns. Indeed, horns are absent from females of 34 other species weighing more than 2.90 kg (Fig. 2C), so attaining the threshold size is no guarantee that natural selection will favour females that carry horns. The alternative to three-parameter power equations, of course, would be two-parameter power equations that of necessity pass through the origin, thereby indicating that all adult bovids, irrespective of their sex or body mass, should possess horns. Third, exponents in the allometric equations for males and females are identical, and they remain unchanged over the full range in body masses (Fig. 2C, D). The perception of Tidière et al. (2017) that the allometric exponent is a variable in the equation fitted to data for males was caused by performing analyses in the logarithmic domain instead of on the arithmetic scale (Packard, 2017c). Fourth, the pattern of variation in horn length differs between males and females (Fig. 2B), but this difference is the result of effects of sex on the allometric coefficient (a). Given the emphasis that is usually given in allometry to the influence of the allometric exponent, it probably will come as a surprise to many readers that the pronounced difference in trajectories for curves in Figure 2B is due to differences in the allometric coefficient and not the allometric exponent. Such a difference would probably go undetected in studies treating data for males and females separately because the allometric coefficient is commonly ignored. Finally, the equations based on model 4 point to sexual dimorphism in horns that becomes more pronounced as species become larger (Fig. 2B). This finding stands in contrast to that reported by Tidière et al. (2017), who contend that species of intermediate size are more dimorphic than species from the extremes of the size distribution. Of course, a study of allometry is no better than the data on which it is based, and the present treatment is no exception. Although observations for horn length in females appear to be quite good (Fig. 2C), those for horns in males are problematic (Fig. 2D). The observations for males seem to comprise two different groups: a group of 85 species weighing less than 300 kg and a second group of six species weighing in excess of ~600 kg. Data for horns in the former group form an off-vertical column at the left of the graph, whereas observations for the latter group lie well to the right and substantially below the fitted curve. Thus, observations for males may not have a unitary pattern amenable to allometric analysis, in which case it would be premature to discuss a pattern of allometric variation that may be spurious. On the other hand, the odd distribution may be nothing more than a trivial anomaly caused by a gap in the sample, in which case the fitted curve may provide a reasonably accurate description for underlying pattern in the data. The question of the bivariate distribution needs to be resolved. However, neither of the aforementioned alternatives lends support to the claim that sexual dimorphism in length of the horns is greatest for species of bovids in the mid-range for body size. SUPPORTING INFORMATION Additional Supporting Information may be found in the online version of this article at the publisher’s website: Table S1. SAS code for implementing ANCOVA when the relationship between Y and X is curvilinear. ACKNOWLEDGEMENTS I thank M. Tidière, A. Nevill and an anonymous referee for thoughtful comments and criticisms that helped me to sharpen the focus of the manuscript. This investigation was funded out of pocket and received no financial support from any governmental or institutional source. REFERENCES Anderson DR. 2008. Model based inference in the life sciences: a primer on evidence . New York: Springer. Google Scholar CrossRef Search ADS   Bales GS. 1996. Heterochrony in brontothere horn evolution: allometric interpretations and the effect of life history scaling. Paleobiology   22: 481– 495. Google Scholar CrossRef Search ADS   Burnham KP, Anderson DR. 2002. Model selection and multimodel inference, 2nd edn . New York: Springer. Burnham KP, Anderson DR, Huyvaert KP. 2011. AIC model selection and multimodel inference in behavioral ecology: some background, observations, and comparisons. Behavioral Ecology and Sociobiology   65: 23– 35. Google Scholar CrossRef Search ADS   Cochran WG. 1957. Analysis of covariance: its nature and uses. Biometrics   13: 261– 281. Google Scholar CrossRef Search ADS   Finney DJ. 1989. Was this in your statistics textbook? V. Transformation of data. Experimental Agriculture   25: 165– 175. Google Scholar CrossRef Search ADS   Fuentes-G JA, Housworth EA, Weber A, Martins EP. 2016. Phylogenetic ANCOVA: estimating changes in evolutionary rates as well as relationships between traits. The American Naturalist   188: 615– 627. Google Scholar CrossRef Search ADS PubMed  Grafen A. 1989. The phylogenetic regression. Philosophical Transactions of the Royal Society of London   326: 119– 157. Google Scholar CrossRef Search ADS PubMed  Grafen A. 1992. The uniqueness of the phylogenetic regression. Journal of Theoretical Biology   156: 405– 423. Google Scholar CrossRef Search ADS   Huxley JS. 1932. Problems of relative growth . London: Methuen & Co. Klingenberg CP. 1998. Heterochrony and allometry: the analysis of evolutionary change in ontogeny. Biological Reviews   73: 79– 123. Google Scholar CrossRef Search ADS PubMed  Lemaître JF, Vanpé C, Plard F, Pélabon C, Gaillard JM. 2015. Response to Packard: make sure we do not throw out the biological baby with the statistical bath water when performing allometric analyses. Biology Letters   11: 20150144. Google Scholar CrossRef Search ADS PubMed  Lolli L, Batterham AM, Kratochvíl L, Flegr J, Weston KL, Atkinson G. 2017. A comprehensive allometric analysis of 2nd digit length to 4th digit length in humans. Proceedings of the Royal Society B   284: 20170356. Google Scholar CrossRef Search ADS PubMed  Miller DR. 1984. Reducing transformation bias in curve fitting. American Statistician   38: 124– 126. Nevill AM, Holder RL. 1994. Modelling maximum oxygen uptake—a case study in non-linear regression model formulation and comparison. Journal of the Royal Statistical Society C   43: 653– 666. Packard GC. 2011. Unanticipated consequences of logarithmic transformation in bivariate allometry. Journal of Comparative Physiology B   181: 841– 849. Google Scholar CrossRef Search ADS   Packard GC. 2014. On the use of log-transformation versus nonlinear regression for analyzing biological power laws. Biological Journal of the Linnean Society   113: 1167– 1178. Google Scholar CrossRef Search ADS   Packard GC. 2015. Allometric variation in the antlers of cervids: a comment on Lemaître et al. Biology Letters   11: 20140923. Google Scholar CrossRef Search ADS PubMed  Packard GC. 2016. Relative growth by the elongated jaws of gars: a perspective on polyphasic loglinear allometry. Journal of Experimental Zoology B   326: 168– 175. Google Scholar CrossRef Search ADS   Packard GC. 2017a. The essential role for graphs in allometric analysis. Biological Journal of the Linnean Society   120: 468– 473. Packard GC. 2017b. Misconceptions about logarithmic transformation and the traditional allometric method. Zoology   123: 115– 120. Google Scholar CrossRef Search ADS   Packard GC. 2017c. Why allometric variation in mammalian metabolism is curvilinear on the logarithmic scale. Journal of Experimental Zoology A   327: 537– 541. Google Scholar CrossRef Search ADS   Packard GC. 2018. A new research paradigm for bivariate allometry: combining ANOVA and nonlinear regression. Journal of Experimental Biology   221: jeb177519. Reeve ECR. 1940. Relative growth in the snout of anteaters. A study in the application of quantitative methods to systematics. Proceedings of the Zoological Society of London A   110: 47– 80. Revell LJ. 2010. Phylogenetic signal and linear regression on species data. Methods in Ecology and Evolution   1: 319– 329. Google Scholar CrossRef Search ADS   Sartori AF, Ball AD. 2009. Morphology and postlarval development of the ligament of Thracia phaseolina (Bivalvia: Thraciidae), with a discussion of model choice in allometric studies. Journal of Molluscan Studies   75: 295– 304. Google Scholar CrossRef Search ADS   Smaers JB, Rohlf FJ. 2016. Testing species’ deviation from allometric predictions using the phylogenetic regression. Evolution   70: 1145– 1149. Google Scholar CrossRef Search ADS PubMed  Sprent P. 1972. The mathematics of size and shape. Biometrics   28: 23– 37. Google Scholar CrossRef Search ADS PubMed  Strauss RE. 1993. The study of allometry since Huxley. In: Huxley JS. Problems of relative growth  [new edition]. Baltimore: Johns Hopkins University Press, xlvii– lxxv. (http://www.faculty.biol.ttu.edu/Strauss/Pubs/Papers/1993StraussHuxley.pdf) Stumpf MPH, Porter MA. 2012. Critical truths about power laws. Science   335: 665– 666. Google Scholar CrossRef Search ADS PubMed  Tidière M, Lemaître JF, Pélabon C, Gimenez O, Gaillard JM. 2017. Evolutionary allometry reveals a shift in selection pressure on male horn size. Journal of Evolutionary Biology   30: 1826– 1835. Google Scholar CrossRef Search ADS PubMed  © 2018 The Linnean Society of London, Biological Journal of the Linnean Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Biological Journal of the Linnean SocietyOxford University Press

Published: May 29, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off