The use of comparative methods (statistical methods that take into account the expected lack of independence of observations sampled across a phylogeny) has become quite common in recent years. The use of such methods is important because conventional statistical methods are based on the assumption that the observations have been randomly and independently sampled (Sokal and Rohlf 1995). The use of nonparametric and resampling methods does not get around this problem, as they are based on this same fundamental assumption. Martins and Garland (1991) distinguish various types of correlations that could be estimated. For the models investigated by them and used here, the appropriate correlation is what they call the “input correlation,” the correlation between the changes on different variables during evolution. Despite the common use of these methods, some aspects of their use and interpretation remain poorly understood and thus controversial. At least part of the controversy seems to be due to several related common misunderstandings of what the use of comparative methods can be expected to accomplish. Westoby et al. (1995b) said that comparative methods can be used to test the null hypothesis that the correlation between traits is due to phylogeny. Desdevises et al. (2003) viewed phylogeny as a confounding variable whose effects must be removed. Beauchamp and Fernandez-Juricic (2004) and many others said that comparative methods are able to partial out the effect of common ancestry or to remove trait variation that is correlated with phylogeny. The purpose of this paper is to show that these are not useful ways to interpret comparative methods.

The discussion that follows will use the properties of the phylogenetic generalized least squares method (PGLS) applied to continuous variables (e.g., Martins 1996; Martins and Hansen 1997). Some of its important properties are summarized below, but see Garland and Ives (2000) and Rohlf (2001) for more details of the properties of this method. The results presented below also apply to the phylogenetically independent contrasts method (PIC) proposed by Felsenstein (1985). This is because PIC represents an implementation of a special case of the PGLS method.

Statistical Analyses

There are two types of statistical analyses that are common in comparative studies of continuous traits. The first is regression, where the interest is in predicting the variation in a dependent variable as a function of one or more independent variables. It is based on the model y = Xb + ϵ, where y is the n × 1 vector of the dependent variable, X is the n × q matrix of the q independent variables (the first column is usually the unit vector in order to include the y-intercept in the model, the other columns vary depending of the independent variables of interest and the statistical design of the study), b is the q × 1 vector of partial regression coefficients, and ϵ is the error term. The number of observations, n, corresponds to the number of species or other taxonomic units. The acronym OTU (for operational taxonomic unit; Sokal and Sneath 1963) will be used here to refer to them.

The second is a multivariate analysis of covariances or correlations among two or more dependent variables. These analyses are based on the model Y = 1_nμ + E, where Y is an n × p matrix of n observations on p dependent variables, 1_n is a unit vector of length n, μ is a 1 × p vector of means, and E is an n × p matrix of multivariate normally distributed errors. Both rows and columns of E may have unequal variances and be correlated. The covariance matrix for E is most easily described if we express Y and E as vectors, y and ϵ, in which the columns of Y and E are stacked one above another to form vectors of length np. The simplest model for the np × np covariance matrix of the observations is Σ = Σ_p ⊗ Σ_n, where Σ_p is the covariance matrix for the p variables, Σ_n is the covariance matrix for the n observations, and ⊗ is the Kroneker or direct product. If correlations among the variables are of interest, then they are obtained from the estimated Σ_p matrix. For both types of analyses, the Σ_n matrix is assumed known, while the other parameters are estimated from the data.

In ordinary least-squares (OLS) the elements in ϵ are assumed to be independent and follow a normal distribution with a mean of zero and a variance of σ² (different dependent variables can, of course, have difference variances). More general methods are required when the covariances among the observations are not zero or the variances are unequal. The simplest approach is generalized least-squares (GLS) analysis. In a GLS model ϵ is assumed to follow a multivariate normal distribution with a mean of 0 (a null vector of length n) and an n × n covariance matrix σ²Σ_n. The diagonal elements of the σ²Σ_n matrix give the variances of each of the n OTUs (a function of the height of each OTU above the root), and the off-diagonals contain the expected covariances between pairs of OTUs (a function of the height above the root of their most recent common ancestor). The covariance matrix is a function of both the phylogeny and an assumed evolutionary model. Most applications are based on the simple continuous Brownian motion (BM) model, in which the diagonals of σ²Σ_n are proportional to the heights of the OTUs above the root and the off-diagonals are proportional to the heights of the most recent common ancestors of each pair of OTUs. The use of PGLS is not restricted to this evolutionary model. It can be used with other models if they allow one to construct a covariance matrix among the OTUs, for example, the Ornstein-Uhlenbeck and ACDC models. See Blomberg et al. (2003) for a description and a comparison of these methods. The choice of such alternative models does not effect the main conclusions of this paper, so only the BM model is discussed here.

See Rohlf (2001) for an explanation of how to construct the Σ_n matrix. For convenience, it will be assumed below that Σ_n is scaled so that the sum of its diagonal elements is equal to n, that is, tr(Σ) = n. Note that in these models, Σ_n specifies the pattern of covariance of the OTUs for all variables. The Σ_p matrix scales the variances of different variables and introduces any correlations between pairs of variables.

Effects of Incorrectly Specifying Σ_N

An important property of GLS estimates of regression coefficients is that they are unbiased even if one uses an incorrect Σ_n (Rao and Toutenburg 1999). This means that using an incorrect phylogeny or a phylogeny with incorrect branch lengths still results in unbiased estimates of regression coefficients. This does not imply that one will obtain the same values for a sample regression coefficient whether or not one correctly takes phylogeny into account. This point was also made by and Pagel (1993), Rohlf (2001) but perhaps its implications have not been given sufficient emphasis.

The effect of using an incorrect phylogeny (either topology or branch lengths) or an incorrect evolutionary model is to use an incorrect estimate, A, for the covariance matrix rather than Σ_n. The GLS estimate of β will become

where the matrices A and X^tAX are assumed to be of full rank (Rao and Toutenburg 1999). While the sample estimates, b and b_A, will usually be different, it is well-known that both are unbiased estimates of the same population parameters. Rao and Toutenburg (1999) gave the following equation for the covariance matrix of the regression coefficients when an incorrect covariance matrix is used This yields larger standard errors of the regression coefficients than if the correct covariance matrix were used. The expected error mean square becomes An important special case is when the lack of independence is ignored and OLS is used. This corresponds to using an identity matrix for A. The estimate of the covariance matrix of the slopes then simplifies to The variances are larger than those given by equation (2). For the simple case where X = 1_n, the ratio the expected variances of the mean using OLS to that using GLS is This means that OLS estimates will be less reliable than GLS estimates.

The expected value of the estimated error variance is

For the case where X = 1_n, the ratio the expected error mean square of OLS to that of GLS is This is consistent with equation (6) of Blomberg et al. (2003) when the scaling tr(Σ) = n is used. For positively correlated observations (as is the case phylogenetic applications), the estimated error variance will be biased downward (i.e., the model will seem to fit better than it actually does). What this means is that the standard errors computed using OLS will be smaller than they should be. Ricklefs and Starck (1996) reported that the PIC standard errors were about 1.2 to 2.0 times larger than those using OLS. This is not a failure of the PGLS method as suggested in that paper. This underestimation of the standard errors and the increased variability among the estimates is what leads to the increased rate of Type I errors when the effects of phylogeny are ignored.

Unfortunately, comparable analytical results for correlation coefficients are more difficult than for regression coefficients. Martins and Garland (1991) reported that several comparative methods were found to yield unbiased estimates of the input correlation even when incorrect phylogenies or branch lengths were used. However, the simulations performed here suggest that the estimates are slightly biased.

Simulations

Simulations were performed to illustrate the results given above for the estimation of regression coefficients. In addition, simulations were necessary to investigate the statistical properties of correlation coefficients because analytical results are much more difficult. Martins (1996) performed similar simulations for estimating correlation coefficients but used much smaller sample sizes that made it difficult to distinguish a bias in the estimates from the effects of sampling error. Díaz-Uriarte and Garland (1996) also reported on an extensive series of simulations but examined Type I error rates rather than bias.

A phylogenetic tree (see Fig. 1) with n = 10 OTUs was selected somewhat arbitrarily for this illustration. The tree was chosen so that the effects of using comparative methods could be detected. The tree was also chosen to allow some independent evolution after each OTU separated from its most recent common ancestor. In other words, it was chosen so that its corresponding phylogenetic covariance matrix, Σ_n, would be neither close to being diagonal nor close to being singular. The exact details of its structure are not very important here. A total of 10,000 simulated datasets were constructed for the regression and correlation analyses. To model a regression analysis, a vector of y values was generated using the model y = Xβ + ϵ, where the first column of the X matrix consisted of unit vector and the second column consisted of a fixed set of n equally spaced values scaled so that their variance was 1.0. Various values of β were used, and ϵ was a vector of length n of multivariate normally distributed deviates with a covariance matrix Σ_n. These deviates do not have to be created by a simulation along a phylogenetic tree (e.g., see Martins and Garland 1991) because a data matrix consistent with this stochastic model can be constructed directly from the model as shown in the Appendix.

Results

An expected result of the regression simulations was the demonstration that the averages would be very close to their parametric values for both OLS and PGLS methods. The means and variances are given in the figure captions. The main effect of using PGLS was to decrease the variability of the estimates (cf. Fig. 2A and 2B). The magnitude of the expected differences in their variances depends on how different the actual tree is from that assumed in the computations. For the present example, the observed ratio of variances of slopes, 2.9814, is very close to the ratio expected (2.9680) using equations (2) and (9). The estimates were also strongly correlated as illustrated in Figure 2C.

The correlation simulations yielded an analogous result. For ρ = 0 the average correlation was close to zero so the correction for bias, equation (5), had little effect. However, in the ρ = 0.7 simulation the average correlation estimated using PGLS was 0.6725. The average correlation increased to 0.6975 when the bias adjusted correlations were averaged. The average correlation using OLS was much lower, 0.6471 (0.6668 even after applying the bias correction). This difference, though small, is statistical significant (for 10,000 observation, the 5% critical values for ρ are 0.6899 and 0.7099). Consistent results (not shown) were obtained from a series of simulations over a range of correlations. The larger bias for OLS is due to the fact that the effective n that should be used in equation (5) is much smaller than the number of OTUs. As before, the most noticeable effect of using PGLS was the reduction in the variance of the estimated correlations (cf. Figs. 3A, 4A and 3B, 4B). The observed ratios of variances of correlations was 1.947 for ρ = 0 and 2.417 for ρ = 0.7. There does not seem to be a theoretical result with which these values can be compared. The estimates were also highly correlated, as illustrated in Figures 3D and 4D.

As found in previous studies, ignoring phylogeny results in a larger number of Type I errors (e.g., Martins and Garland 1991; Díaz-Uriarte and Garland 1996). For the β = 0 simulation, the Type I error rate using OLS was 37.5% at the nominal 5% level. The error rate dropped to 4.91% when PGLS was used. For the ρ = 0 simulation, the Type I error rate using OLS was 21.20% at the nominal 5% level. The error rate dropped to 5.13% when PGLS was used. The error rates using PGLS were not significantly different from 0.05.

Discussion

The fact that OLS and methods equivalent to PGLS yield unbiased estimates of regression coefficients has been known for some time (e.g., Pagel 1993). Correlations estimated using OLS and PGLS methods have often been found to be similar. For example, Martins and Garland (1991) reported that all of the estimation methods and evolutionary models they investigated yielded unbiased estimates of the correlation between a pair of traits. Though small, almost all (24 of the 27) of the deviations they reported in their table 4 were negative. The deviations for the PIC method are consistent with the expected bias in correlations due to small sample size (eq. 5). The deviations corresponding to the use of OLS (which they called TIPS) were greater, as found in this study.

It is interesting to note that the direction of the small bias when using OLS rather than PGLS is in the opposite direction that many seem to expect. For example, Westoby et al. (1995a,b) and Garland et al. (1999) said that comparative methods typically reduce the observed correlation. Ricklefs and Starck (1996) concluded that phylogenetic methods tended to result in weak correlations with broad confidence limits. Gregory (2003) found no consistent trend in his survey. These surveys were all small with unknown sampling biases, which makes it difficult to reach a conclusion about any empirical difference in average correlation. Papers discussing why comparative methods should be used often show examples where the correlation is reduced after comparative methods are applied (e.g., Felsenstein 1985; Gittleman and Kot 1990; Harvey 2000). Perhaps that is what has lead users to assume that these methods would remove the effects of phylogeny from data correlation.

Note that the advantage of using comparative methods such as PGLS is statistical, that is, if the assumptions of the model on which they are based are correct, then their routine application should result in a user experiencing the expected Type I error rates and power. Their use cannot, of course, ensure that a better estimate will be obtained in any particular dataset. In the simulations presented here, the estimates obtained using PGLS were closer to their true values than the OLS estimates in about 70% of the samples. In an application of these methods one does not know whether a dataset corresponds to one in which PGLS gives a closer estimate of the unknown parameter or is one of the 30% in which OLS yields a closer estimate. For that reason, one must be cautious about trying to interpret the fact that one type of correlation is larger than the other (even if they seem to be significantly different statistically).

Several authors have suggested applying comparative methods to some variables in a study but not to others. However, in the models on which comparative methods are based, the adjustments are to take into account covariances among the observations (whether univariate or multivariate) not the variables. Abouheif (1999) suggested first testing variables and not using comparative methods on them if their variation seemed to be independent of the phylogeny. This introduces the additional complication of using conditional tests. The second test is applied in a way that depends on whether the first test is significant (which could be just a Type I error). The statistical properties of such procedures can be difficult to analyze. Rheindt et al. (2004) also argue that comparative methods should not be applied to all variables but only for those variables for which there is evidence of a phylogenetic signal. They computed regressions where phylogeny had been taken into account for one variable but not the other. In both cases, the problem being addressed seems to be that different variables can have different amounts of nonphylogenetic variation (whether environmental or just measurement error) superimposed. A better solution would be to use methods explicitly based on more complex models. The model on which PGLS is based assumes that evolution can be described as a continuous multivariate random walk and that all variation is phylogenetic, that is, it is a function of the properties of the Σ_n matrix. The Σ_n matrix can include variation within a species (by extending the lengths of the terminal branches), but it applies to all variables equally. This obviously represents an oversimplification.

The phylogenetic mixed model (PMM) developed by Lynch (1991) and extended by Housworth et al. (2004), distinguishes between heritable (phylogenetic) and residual (environmental) components of variation. Felsenstein (2004) outlines another approach to incorporating the effects of sampling error. Of course, by introducing more parameters to be estimated from the data, it is expected that larger sample sizes (more extensive phylogenies) will be required. It is surprising that in this case the introduction of so few additional parameters seems to require much larger sample sizes. A problem with a requirement of larger phylogenies is that one has to assume the regressions or correlations being estimated are homogeneous over larger taxonomic groups (allowing for different relationships in different parts of the phylogeny may be counter-productive as it introduces additional parameters to be estimated). An example of this heterogeneity is given in Garland and Ives (2000). To avoid this, one might just add closely related species but that is not very efficient statistically since the observations are not independent. Housworth et al. (2004) report that the PMM also yields unbiased estimates so that much of the discussion given above about the statistical properties of the PGLS method should also apply to the PMM method. Paradis and Claude (2002) proposed a different generalization using generalized estimating equations. They retain the simple evolutionary model of the PGLS method but allow variables to follow distributions other than just the normal. As mentioned above, these methods are not limited to simple BM models of evolution as long as it is possible to derive an expected Σ_n matrix. Such model-based approaches seem more promising than the various ad hoc adjustments that have been proposed for Felsenstein's (1985) independent contrasts method.

While it does not seem useful to view comparative methods as methods to correct a correlation or regression coefficient for the effects of phylogeny, such methods should still be routinely used because allowing for known dependencies among the observations provides estimates that will have smaller standard errors. This means that when β or ρ is equal to zero they can be expected to lead to significance tests that have Type I error rates closer to the desired nominal rates and tests that are more powerful when these parameters are not equal to zero. This will be true even if the phylogeny is not known exactly. A tree with multifurcations corresponding to unresolved nodes (soft polytomies) and branch lengths that are only approximate should still provide better estimates of β and ρ than when phylogeny is completely ignored by using OLS methods (because a covariance matrix based on a preliminary phylogeny is likely to be closer to Σ_n than an identity matrix). This is an issue that can be investigated using equations (7) and (8) to compare the effects of using alternative phylogenies.