A new research paradigm for bivariate allometry: combining ANOVA and non-linear regression

The analysis of covariance (ANCOVA) is a powerful statistical tool that combines the analysis of variance (ANOVA) with linear regression (Fisher, 1932). The primary applications for ANCOVA are (1) to increase precision for the response variable (Y) in planned experiments on subjects that are randomly assigned to two or more treatments, and (2) to remove confounding effects of a second variable in descriptive studies comparing individuals in defined groups or classes (Cochran, 1957). In both these applications, primary interest is in variation of the response variable among groups; the relationship between the response variable and the potential covariate (X) is of no particular interest (Cochran, 1957). Thus, inclusion of the covariate in the analysis is intended to reduce ‘background noise’ in the response and thereby promote the detection of differences in Y among the groups (Packard and Boardman, 1987).

Biologists commonly recognize a third application for the ANCOVA – namely, its potential for use in studies of allometric variation, where interest is primarily in the relationship between the response variable and the covariate (White, 2003). In such investigations, ‘noise’ may actually be introduced by different structural relationships for data from different groups (Cochran, 1957). Thus, ANCOVA is used in this application to identify the scaling relationship between Y and X after removing the confounding influence of the different groups, or to identify different scaling relationships among the several groups.

The use of ANCOVA in studies of allometry has typically entailed logarithmic transformation of both the response variable and the scaling variable (i.e. the covariate), for two reasons. First, transformation has been standard practice in studies of bivariate allometry since early in the last century (Packard, 2017c). And second, statistical software for performing ANCOVA requires that a linear relationship exists between the response variable and the scaling variable (e.g. Cochran, 1957; White, 2003). Inasmuch as logarithmic transformation often linearizes a distribution that is curvilinear in arithmetic domain, transformation may thereby satisfy a requirement of the analytical procedure while simultaneously conforming to standard practice.

The study of allometric relationships by ANCOVA has been severely constrained, however, by the logarithmic transformation, because such analyses require that the distribution for Y and X in each group conforms to that of a two-parameter power equation with lognormal, heteroscedastic error on the original, arithmetic scale (Packard, 2017c). Many distributions follow the paths of three-parameter equations instead of two-parameter equations (Packard, 2015, 2016, 2017a), and random error may be normal and homoscedastic or normal and heteroscedastic instead of lognormal and heteroscedastic (Packard, 2016). Such distributions are not amenable to analysis by ANCOVA on transformations.

My objective here is to introduce readers to a method for performing the equivalent of ANCOVA on data that are curvilinear in arithmetic domain. The procedure can be applied to observations that are linearized by transformation and to those that are not (i.e. ‘complex allometry’, sensu Strauss, 1993). The procedure combines ANOVA with non-linear regression in a routine implemented in SAS (version 9.3 or later). The protocol was introduced by Lolli et al. (2017) to adjust data for curvilinear allometry in different groups of subjects. Here, I show how the protocol can be used also to evaluate and compare patterns of allometric variation in different groups. The procedure differs from the one used in earlier studies (see Nevill and Holder, 1994) in that the new protocol can assume error to be normal and homoscedastic, normal and heteroscedastic, or lognormal and heteroscedastic. For simplicity, I work with two groups, although more could be accommodated in more complicated designs.

where the first parenthetical expression represents the power function and the second parenthetical expression identifies the main group effect. Y is the response variable, X is the scaling variable and X2 is the group to which a subject is assigned (either 0 or 1 in the present example). 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 Y₀ the intercept with the Y-axis.

The equation is more complex than the one used in conventional ANCOVA (1) because curvature for a power function is affected by both the allometric coefficient and the allometric exponent and (2) because we also are concerned with whether or not an explicit, non-zero intercept is needed. Moreover, group membership directly affects the allometric coefficient (a) and the intercept (Y₀), both of which appear inside the first set of parentheses defining the three-parameter power equation; and group membership also affects the exponent in the equation via the additive effect of the parameter d. The basic design can be modified in various ways to explore a range of questions and statistical fits. For example, two-parameter power functions are readily accommodated by dropping Y₀ from the equation, and the fit of equations with a common exponent can be assessed by dropping the parameter d. By moving Y₀ outside the parentheses, the equation estimates a common intercept for both the fitted functions. In addition, different forms for random error can be readily incorporated into the full model (Lolli et al., 2017; Packard, 2016, 2017a).

The power and utility of this approach can be illustrated by applying it to a dataset taken from the current literature. Accordingly, values for field metabolic rate (FMR) and body mass for 80 species of placental mammal and 35 species of marsupial were recovered from the online appendix to the paper by Capellini et al. (2010). The observations were submitted to multiple analyses, starting with the full allometric equation described above. Subclass was treated as a categorical variable (0=marsupials, 1=placentals) in a variety of two- and three-parameter power equations having different forms for random error (i.e. normal and homoscedastic, normal and heteroscedastic, lognormal and heteroscedastic). However, a large majority of the models (including all models with normal distributions for error) were rejected when Akaike's information criterion (AIC) was found to be 15+ higher than AIC for the best model in the pool of candidates (Anderson, 2008; Burnham and Anderson, 2002; Burnham et al., 2011). SAS code for the 10 most relevant analyses (all with lognormal error) is available in Script 1 so that readers can gain an appreciation for the kinds of models that were explored.

Three statistical models merit consideration, and all three share important features (Table 1). First, all three include terms for Y₀, a, b and d (Table 1); and second, these parameters enter the models at significant levels (Table 2). Taken together, these findings point to important effects of group membership on the allometric exponent (via b and d), and to the necessity for including an explicit intercept in the fitted equations. Inclusion of a parameter for a discrete group effect (c) in models 1 and 2 enhances the goodness of fit, but the term may not be essential (see model 3).

Although model 3 is a plausible alternative, the best fits were for model 1 and model 2, which are statistically equivalent (Table 1). However, low AICs for these models do not guarantee that the fits to the data are good, because all the models in the candidate pool may have been bad (Anderson, 2008; Burnham and Anderson, 2002). Models 1 and 2 also are abstractions in that each of them (with five estimated parameters) actually describes the pattern in two different subclasses of mammal (each requiring a three-parameter power function). The only way to assess goodness of fit for the models, therefore, is to evaluate the predictive equations for each of the mammalian groups and present them graphically (Anscombe, 1973; Packard, 2017b).

The predictive equations derived from models 1 and 2 track observations for marsupials and placentals quite well (Fig. 1A,C), and standardized residuals for the models indicate that overall fits are good (Fig. 1B,D). Thus, it is reasonable to conclude from analyses of the available data that patterns of allometric variation in both subclasses of mammal require a positive intercept with the Y-axis. A single intercept for both functions is indicated by model 1, but the two proximate intercepts predicted by model 2 are equally likely. The presence of an intercept in both the functions indicates that the distributions for both marsupials and placentals are curvilinear on the logarithmic scale (see Packard, 2017d; Sartori and Ball, 2009). Such curvilinearity in the log domain is commonly taken to mean that a two-parameter power function describing untransformed data has an exponent which varies with body size (e.g. Strauss, 1993), yet this is clearly not the case with FMRs of mammals (Fig. 1A,C).

It also is safe to say that the allometric exponent is higher for placentals than for marsupials, thereby contributing to differences in trajectories of the curves (Fig. 1A,C). However, models 1 and 2 make rather different predictions for the allometric exponent for marsupials (Fig. 1A,C). Model 1 predicts an exponent of 0.66, which is close to the 0.67 that is expected with geometric scaling, whereas model 2 predicts an exponent of 0.71, which is tantalizingly close to the 0.75 that is expected with quarter-power scaling. Nevertheless, the models are equivalent by AIC. It is noteworthy, however, that the location and trajectory for each of the mean functions is affected by the allometric coefficient (a) as well as by the allometric exponent, and that the two parameters interact in reciprocal ways in the alternative models to bring about virtually identical outcomes (Fig. 1A,C). The very fact that these dissimilar estimates could come from equivalent statistical models fitted to a single dataset is not altogether surprising because curvilinearity is expressed differently in the alternative formulations: the two functions estimated by model 1 have the same intercept whereas different intercepts are required for functions estimated by model 2 (see Arah, 2008, for discussion on covariates).

With all the attention that is given in bivariate allometry to the value for the exponent in the allometric equation (or to the slope for a straight line fitted to log transformations), it is easy to forget that the trajectory for the fitted curve (and therefore predictions for the response variable) depends also on the value assumed by the allometric coefficient. Indeed, most investigations simply ignore the coefficient altogether. The protocol used here simultaneously recovers information on both the allometric coefficient and the allometric exponent, so the method represents an important advance over methods based on separate statistical models fitted to data for each of the groups. For example, the significant term for d in model 1 indicates that the exponent differs between allometric equations for marsupials and placentals (Table 2), while the significant term for c in the model points to allometric coefficients that also differ between the mammalian groups (Table 2). Simply stated, the group effect (indicated by c) modifies the expression of the allometric coefficient (a) in equations for both marsupials and placentals.

An earlier treatment of allometric variation in FMR of marsupials was inconclusive, with three equally good alternative models being fitted to the observations (Packard, 2017a). When the same data were included in a combined analysis together with data for placentals, greater discriminative power was achieved, and a three-parameter equation was favored (Table 1). The point, of course, is that the pattern of allometric variation for marsupials could be resolved when the data were combined with those for placentals in what amounts to being a one-way ANCOVA.

To my knowledge, no method is currently available in SAS to incorporate an adjustment for phylogenetic relationships into non-linear regressions with different forms for random error. This circumstance will doubtless be a matter of concern to workers who are engaged in studies of interspecific variation, but it should not be a deterrent to adopting the methodology described here (at least for preliminary screening of potential models). The ability to explore a variety of different power equations with different forms for random error (both additive and multiplicative) should have priority over adjusting for phylogeny in conventional ANCOVA on logarithmic transformations (e.g. Fuentes-G et al., 2016). After all, phylogenetic methods merely tweak residuals in the proposed statistical model so that it fits the data better than does a model without a correction for phylogeny (Revell, 2010). Incorporating phylogenetic signal into such analyses is of little consequence if the underlying descriptive equation is wrong.

In summary, the method described here for performing the equivalent of non-linear ANCOVA provides investigators with a powerful tool for empirical studies of allometric variation in two or more groups of organisms. Two- and (or) three-parameter power functions can be fitted to subjects in each group, and the different parameters in the functions can be evaluated and compared. Random error can be modeled as normal and homoscedastic, normal and heteroscedastic, or lognormal and heteroscedastic, thereby avoiding some of the controversy that has plagued allometric research in recent years (Packard, 2017c). In addition to the obvious application of the protocol to purely descriptive research, the approach has the potential for testing explicit hypotheses about variation in the response variable among groups (Lolli et al., 2017).

I am greatly indebted to Alan Batterham (Teesside University, UK), who shared with me his insights concerning the handling of categorical variables in non-linear regression. I am also grateful to multiple referees who offered constructive criticisms that helped me to improve the manuscript.