Emmanuel Paradis
email: emmanuel.paradis@mpl.ird.fr

ASYMMETRIES IN PHYLOGENETIC DIVERSIFICATION AND CHARACTER CHANGE CAN BE UNTANGLED

Keywords: .-diversification, extinction, maximum likelihood, phylogeny, speciation

 pointed out that uncareful use of some of these tools could lead to misleading conclusions on diversification or character evolution, and thus to difficulties in distinguishing both phenomena. I here present guidelines for the analyses of macroevolutionary data that may help to avoid these problems. The proper use of recently developed statistical methods may help to untangle diversification and character change, and so will allow us to address important evolutionary questions.

The last twenty years have witnessed a remarkable change in paradigm for evolutionists: the variation in speciation and extinction rates (i.e., the tempo of evolution), and the variation in the rates of character change (the mode of evolution) are now ideally studied with molecular phylogenies and data collected on recent species. These issues, traditionally called macroevolution, were the domain of paleontologists during several decades [START_REF] Simpson | The Major Features of Evolution[END_REF], whereas molecular data were used to address microevolutionary mechanisms (to have an idea of this changing paradigm, compare the two editions of the same book by [START_REF] Futuyma | Evolutionary Biology[END_REF][START_REF] Futuyma | Evolutionary Biology[END_REF].

In a recent paper, [START_REF] Maddison | Confounding asymmetries in evolutionary diversification and character change[END_REF] used simulated data following two scenarios to bring attention to some limits of this new paradigm. In the first scenario, some clades were simulated starting from the root, and a character evolved along the tree: two states were allowed (0 and 1) with a constant probability of change between them.

The speciation rate was related to the state of the character so that species in state 1 split at a higher rate than those in state 0. Maddison showed that the estimates of the ratio of character transition rates was correlated with the ratio of the speciation rates. In the second scenario, Maddison used a similar setting except that the speciation rate was constant, but the rates of transition of the character in one direction was four times higher than in the other. He then showed that the analysis of diversification led to infer, more frequently than by chance, that the speciation rate was different between the states of the character. [START_REF] Maddison | Confounding asymmetries in evolutionary diversification and character change[END_REF] concludes that, if simple data analyses are done, biased diversification with respect to a character may lead to wrong conclusions on the evolution of this character. On the other hand, biased character evolution may lead to falsely associate this character to biased diversification. In other words, this raises the question whether we may be able to untangle differential evolution of character from differential diversification. The question may be framed in more specific terms from Maddison's simulation study: if a character state is observed to be relatively rare in a clade, can we distinguish whether this state is associated with a low speciation rate, or evolution towards this state occurs at a low rate? (Or both?) This issue is of great importance because, if we can answer positively, we could point which effect most contributes to the generation of diversity. Heterogeneous diversification rates and heterogeneous character change rates are likely linked to different evolutionary mechanisms. Character states with a high transition rate may be the result of counter-selection, developmental instability, or low plasticity of the character. On the other hand, a high speciation rate related to a character state may be the result of its adaptative value, or its association with the rate of reproductive isolation.

My aim in the present paper is to show that appropriate data analysis methods can separate the effects of biased character change and biased diversification. This suggests that heterogeneous character evolution and heterogeneous diversification rates can be untangled in future analyses of macroevolutionary data. I point to some recommendations for data analyses in macroevolutionary studies with recent species, as well as further needs for future research in this area. (2006) showed that the maximum likelihood estimates of the ratio of character transition rates was correlated with the actual ratio in speciation rates. One may be tempted to interpret a high ratio of character transition rates as evidence for a strong bias in character change, thus concluding that rates are actually different. However, Maddison did not address the issue of testing whether these rates are significantly different. Indeed, looking at parameter estimates is not the ideal way to assess the validity of an hypothesis. The results from Maddison's analyses illustrate that these parameter estimates can be strongly biased when a wrong model is used, not that they are significantly different.

Methods

Analysis of Character Change

Maddison

With simulated data, we know the model used to generate the data, but with real data we have to test whether a given model is appropriate. In a likelihood framework, the likelihood ratio test (LRT) is the canonical way for hypothesis testing [START_REF] Lindsey | Parametric Statistical Inference[END_REF][START_REF] Ewens | Statistical Methods in Bioinformatics: An Introduction[END_REF]. To fix ideas, let us write explicitly the model used to simulate the data in Maddison's first scenario (this model is referred to as model 1 in the text); its rate matrix (often denoted Q) is:

0 1 0 1    -r r r -r    , (1)
where the row labels denote the initial states of the character (denoted x in this paper), the column labels the final ones, and r is the instantaneous rate of change, that is the probability that x changes of state during a very short interval of time. It is difficult to interpret the parameter r biologically, but its virtue is that it is independent of time. To obtain real probabilities, we have to fix a time interval, say t, and compute the matrix exponential e tQ (see below).

Maddison estimated the rates ratio using the following two-parameter model (referred to as model 2):

   -r 1 r 1 r 2 -r 2    . (2)
We know, in this particular situation, that this model has an extra parameter because the data were simulated with r 2 = r 1 . Models 2 and 1 can be compared with a LRT because they are nested: the latter is a particular case of the former [START_REF] Pagel | Detecting correlated evolution on phylogenies: a general method for the comparative analysis of discrete characters[END_REF][START_REF] Nosil | Testing hypotheses about ecological specialization using phylogenetic trees[END_REF]. The test follows a χ 2 distribution with df = 1 (the difference in number of parameters).

In real situations we cannot be sure that either model 1 or model 2 generated the observed data. If neither of them are the true model, the tests of hypotheses could be biased in the same way than Maddison showed that parameter estimates are biased.

It is thus necessary to assess the goodness-of-fit of a model in a general way. In the present context, the shape of the likelihood function is an indication of the poor or good fit of a model. If a model poorly fits the data, the likelihood function is relatively flat. It is possible to examine the shape of the likelihood function by plotting it against a range of parameter values, but an easier procedure is to look at the estimated standard-errors of the parameters which are derived, under the standard likelihood method, from the second derivatives of the likelihood function.

Thus, the smaller these standard-errors, the narrower the likelihood function. When no analytical expression of the second derivatives are possible, which is precisely the case for the Markovian models considered here, a numerical computation may be done using nonlinear optimization [START_REF] Schnabel | A modular system of algorithms for unconstrained minimization[END_REF]. The ratio of the parameter estimate, r, on its standard-error, se(r), can be used as an indication of the shape of the likelihood. This ratio is in fact a formulation of the Wald test which, under the null hypothesis that r is equal to zero, follows a standard normal distribution (see [START_REF] Rao | Linear Statistical Inference and its Applications[END_REF]:

r se(r) ∼ N (0, 1) if r = 0.
The LRT and the Wald test can thus be used to test the same hypothesis, but the latter is rarely used because it has generally poorer statistical performance than the LRT, particularly for small sample sizes [START_REF] Agresti | Categorical Data Analysis[END_REF][START_REF] Mcculloch | Generalized, Linear, and Mixed Models[END_REF].

However, both tests are expected to give the same results for large sample sizes [START_REF] Rao | Linear Statistical Inference and its Applications[END_REF][START_REF] Mcculloch | Generalized, Linear, and Mixed Models[END_REF].

I propose the following rule of thumb: when this ratio is less than two, it is likely that the model under consideration is not appropriate. A justification for this rule is that, under the assumption that the maximum likelihood estimates are normally distributed, (a standard assumption of the likelihood theory of estimation), then a 95% confidence interval may be calculated with r ± 1.96 × se(r). Consequently, if r/se(r) < 2 then zero is included in this interval suggesting the presence of a flat likelihood function. The method should be used as follows. First, perform the LRT comparing both models. If this test is significant, then examine the ratios of the rate estimates of model 2 on their standard-errors: if one of them is less than two, then it is likely that the LRT is biased and the null hypothesis is true.

Analysis of Diversification

Maddison (2006) inferred differences in speciation rates by comparing the number of species between sister-clades that are different with respect to the character, but all species within each clade have the same character state. This method does not use all available information as it considers only a subset of the nodes of the tree, and it ignores branch lengths. Several methods have been proposed to analyze diversification that essentially differ in the type of data under consideration (see reviews in [START_REF] Sanderson | Reconstructing shifts in diversification rates on phylogenetic trees[END_REF][START_REF] Mooers | Evolutionary process from phylogenetic tree shape[END_REF][START_REF] Pagel | Inferring the historical patterns of biological evolution[END_REF]. For instance, some methods consider tree topology and balance (e.g., [START_REF] Aldous | Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today[END_REF], whereas others consider the distribution of branch lengths (e.g., [START_REF] Pybus | Testing macro-evolutionary models using incomplete molecular phylogenies[END_REF]. When the topology and branch lengths of the analyzed phylogeny are available, it is possible to refine the analyses and use more elaborate methods. Recent developments have been done in the inference of diversification, particularly the Yule model with covariates which takes full account of all phylogenetic information (tree topology and branch lengths) to infer the effects of species traits on speciation rates [START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF]. In this model the speciation rate depends on a linear combination of some variables measured on each species. This approach is similar to a standard linear regression where the mean of the response is given by a linear combination of variables.

Consequently, a wide variety of models may be fitted to the same phylogenetic and species traits data. The latter could be continuous and/or discrete.

Because the Yule model with covariates is fitted by maximum likelihood, different models can be compared with a LRT if they are nested. In the present context of testing the effect of x on the speciation rate (λ), the general model is:

log λ 1 -λ = βx + α,
where α and β are parameters, so that the right-hand side of the above equation is equal to α if x = 0, or to α + β if x = 1 (see [START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF], for details). The LRT comparing this model to the standard Yule model tests the hypothesis of the effect of x on λ, i.e., whether β is significantly different from zero. The Wald test may also be performed by computing β/se( β). The same criterion proposed above for the analysis of character was applied as well.

Data Simulation

I simulated some data with known parameters in order to assess the statistical performance of the methods described above. The idea was to generate some data sets where the majority of the species are in a state due to, either different speciation rates, or different rates of character change, or both. The data were then analyzed to assess whether the two processes can be untangled. The simulations were started from the root of the tree, that is the initial bifurcation. At each time-step, each species present in the clade had a given probability (λ) to split into two, and a probability (p) to change the state of its character x. When a species split, the daughter-species inherited its value of x. Both λ and p depend on x, and so are hereafter denoted λ 0 , p 0 , λ 1 , and p 1 , where the subscript indicates the value of x.

Two of these parameters were fixed for all simulations: λ 0 = 10 -4 and p 0 = 5 × 10 -5 .

Three different combinations of parameters were used for λ 1 and p 1 : (1) λ 1 = 5λ 0 , p 1 = p 0 , (2) λ 1 = λ 0 , p 1 = p 0 /10, and (3) λ 1 = 2.5λ 0 , p 1 = p 0 /4.

These settings correspond to the three plausible biological scenarios leading to the abundance of a trait in a clade: (1) this trait is associated with a high speciation rate,

(2) species without this trait tend to evolve towards acquiring it, and (3) a mixture of both processes. The parameter values were chosen so that ca. 90% of species had

x = 1 at the end of the simulation. These were found analytically in setting (2), since speciation was homogeneous, using the matrix exponentiation explained above, and the fact that the expected number of species after t time-steps is given by 2e λt [START_REF] Kendall | On the generalized "birth-and-death" process[END_REF]. It was obviously unnecessary to consider here the null setting λ 1 = λ 0 and p 1 = p 0 because this would yield 50% of species in each state, and so little difficulty for data analysis. The first setting is similar to [START_REF] Maddison | Confounding asymmetries in evolutionary diversification and character change[END_REF] first scenario, whereas the second setting is close to his second one: the difference is that he used a ratio of 4 instead of 10.

The time-step of the simulations was transformed in time unit equal to 0.001.

Giving a probability of change of 5 × 10 -5 for t = 0.001, we can find by back-transformation with the matrix exponentiation and interpolation that the actual parameter of the first setting was r ≈ 0.05. All simulations were run until 100 species were present. This was replicated 100 times for each possible initial state at the root (0 or 1) and each combination of the parameters. The trees and values for x at the end of the simulation were saved (the files are available from the author). All simulations were programmed in R version 2.4.0 (R Development Core Team 2006).

The data were analyzed as if they were real data: the procedures described above were applied to all replicates using R's standard looping functions. Both LRTs and Wald tests were computed, as well as the proportion of cases where the tests agreed in rejecting the null hypothesis. The analyses of character change and of diversification were done with the package ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]).

Results

The results were slightly affected by the state of the root as it changed slightly the proportion of species in state 1 at the end of the simulation: 83.8% and 91.6%, for a root in state 0 or 1, respectively (data pooled over all simulations). This proportion had in fact a skewed distribution when calculated for each replicate; the corresponding medians were thus slightly larger: 85% (range: 28-99) and 92%

(63-99). For simplicity, the results below are presented for the two series of simulations pooled since they were overall consistent.

Analysis of Character Change

The proportions of replicates where the LRT comparing models 1 and 2 was significant at the 5% level were 0.18, 0.405, and 0.175 for the three scenarios of equal Table 1 rate of change, strong, and moderate biases, respectively (Table 1). On the other hand, the Wald test had relatively poor performance as the rejection rate of the null hypothesis was almost the same whatever the scenario (ca. 0.25). After examining the standard-errors of the parameter estimates, the proportions of cases where the LRT was significant and the ratios r1 /se(r 1 ) and r2 /se(r 2 ) were greater than 2 were lowered to 0.055, 0.26, and 0.1. Only the first one was not significantly different from 0.05 (two-sided exact binomial test: P = 0.744, P < 0.0001, and P = 0.003, respectively).

Figure 1 gives a summary of the distributions of the estimates under both models Fig. 1 for the three settings. As observed by [START_REF] Maddison | Confounding asymmetries in evolutionary diversification and character change[END_REF], a few cases resulted in extremely dispersed estimates, so I focused on the median and first and third quartiles. In the first setting, the median of the estimator r was very close to the true value (0.049, instead of 0.05). Remarkably, the dispersion of all parameter estimates and their standard-errors were smaller when model 2 was true.

The reconstruction of ancestral states is also informative. Under model 2 nearly all nodes (except the most terminal ones) had a relative likelihood of 0.5 for each state, meaning that these reconstructions were in fact indeterminate under this model. Under model 1 most nodes were estimated to be in state 1 with a high probability (relative likelihood greater than 0.9): this was particularly the case for the root even if the actual state was 0 (Fig. 2). Thus the models failed to estimate Fig. 2 correctly the ancestral states of the simulated character.

Analysis of Diversification

The simulated trees and phenotypic values were analyzed with the Yule model with covariates testing for the effect of x on λ. The method requires us to reconstruct the values of x at the nodes of the tree: this was done with a simple parsimony criterion.

The branch lengths were back-transformed to the original time scale of the simulations prior to analysis in order to ease the fitting procedure. The proportion of significant LRT were 0.415, 0.01, and 0.13. All these proportions are significantly different from 0.05 (two-sided exact binomial test: P ≤ 0.005 in all cases). The Wald test gave similar results, but with slightly larger rejection rates of the null hypothesis.

However, in all cases where the LRT was significant the Wald test was as well (Table 1). It is interesting to note that though the estimated ancestral states were inaccurate, the test was able to reject the null hypothesis in more than 5% of the cases, and so is robust, at least partially, to inexact ancestral character estimation (at least in the present situation).

Discussion

The combined use of the LRT and Wald test here revealed an interesting contrast between analyses. In the analysis of character change, the proposed criterion derived from the Wald test helped to keep the type I error rate at a reasonable value. In the analysis of diversification, the LRT performed well and there is, in the present scenarios, no need to use the proposed criterion. How to explain this difference?

There are three possible, nonexclusive explanations of the relatively poor performance of the models of character change. Firstly, the departures from the Markovian assumptions may be too strong so that the models considered here are too 'ill-defined' for the present data. Secondly, the maximization of the likelihood functions of these models is not straightforward as it requires iterative calculations of likelihoods along the nodes of the tree which may be difficult for current optimization methods. By contrast, the likelihood function of the Yule model with covariates is a linear function of the parameters [START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF]. Thirdly, the current formulation of the model of character change may not be adequate for likelihood maximization, and a reparametrization resulting in a linear combination of parameters might be more appropriate. [START_REF] Mooers | Reconstructing ancestor states with maximum likelihood: support for one-and two-rate models[END_REF] were cautious about the approximation of the LRT with a χ 2 distribution, and recommended using conservative values to assess the increase in likelihood when fitting model 2. They also called for more work to verify the validity of this approximation. In the present study, I observed in some cases a difficulty to fit the models which required to repeat the fitting procedure with different starting values. Some further study is needed to assess whether this is due to the inadequacy of these models, or a failure of our current algorithms of likelihood maximization.

All simulations were run with parameter values that gave equal number of species, and similar frequencies of the states of x. So the differences in the results of the tests were due to the internal structure of the data, not to the prevalence of species with x = 1. The Yule model with covariates will tend to be accepted when shorter branches are associated with a character state. Such an association may be significant even if a state is at low frequency. Another interesting property of the Yule model with covariates revealed here is that it is robust to inexact estimation of ancestral states. This is likely to be important giving difficulties in reconstructing ancestral states (see below).

The simulations strongly suggest that the methods considered here are statistically consistent: they had greater power to reject the null hypothesis when the difference between the parameters were larger. [START_REF] Mooers | Reconstructing ancestor states with maximum likelihood: support for one-and two-rate models[END_REF] suggested that most data sets are too small to yield enough statistical power to reject model 1.

Interestingly, most data sets they considered have many fewer than 100 species. With the increasing size of biological databases, and the advent of supertree techniques (e.g., [START_REF] Agnarsson | Sociality in theridiid spiders: repeated origins of an evolutionary dead end[END_REF][START_REF] Beck | A higher-level MRP supertree of placental mammals[END_REF][START_REF] Bininda-Emonds | The delayed rise of present-day mammals[END_REF], we are likely to have now the possibility to fit more complex models of character change.

A possible limitation of the current Markovian models of character change is the assumption of equilibrium [START_REF] Nosil | Testing hypotheses about ecological specialization using phylogenetic trees[END_REF]. When this is violated, and this is certainly true in real situations, then our methods will lose some power. However, the simulations presented here suggest that this is apparently a question of degree.

Though the method loses some power, it remains statistically consistent. In a more general context, the statistical consistency of data analysis methods is important in the face of confounding factors in realistic situations. The critical point is that the methods can still detect significant effects even when their assumptions are not met.

An important confounding factor not considered here is extinction and the variation of its rate. If extinction rate is related with a species character, then this will likely result in biases for both methods considered here because this character will be associated with long branches in the tree. With respect to the Yule model with covariates, a previous simulation study showed that random extinction resulted in a loss of statistical power but the method remained statistically consistent [START_REF] Paradis | Statistical analysis of diversification with species traits[END_REF].

The present paper considered separate analyses of character change and diversification with currently available methods. From a statistical point of view, this can be viewed as each analysis was conditioned on the other: the analysis of character change was done assuming constant and homogeneous diversification, whereas the analysis of diversification was done assuming given reconstruction of the ancestral states. An interesting perspective would be to develop a joint analysis of these processes. Since each analysis is done by maximum likelihood, it is possible to combine both likelihood functions in a joint likelihood function that would be maximized over all parameters (r and λ). It remains to be seen whether this would increase the statistical performance of our methods.

Estimation of Ancestral States

The results of the analyses of character change show that the inferred state at the root is likely to be misleading (see also Mooers and Schluter 1999, Fig. 1b). This failure to estimate correctly the state at the root is due to the assumption of equilibrium common to most Markovian models. A consequence of this assumption is that whatever the initial state, the probability to be in any state is given by its relative transition rate. For model 1, these probalilities are 0.5 since the transition rates are equal. This is concretely visualized by calculating e tQ whose values tend to 0.5 when t increases. Consequently, if one state is rare this implies that the initial state (i.e., the root) was in the other state, and so the system is in transition towards its equilibrium.

On the other hand, it was observed that the estimates of the rate of change were nearly unbiased. This suggests that, while we are able to correctly assess how frequently a character changed, our estimates of its ancestral states may not be meaningful. This is an important issue because there is more interest in ancestral states than on rates of change (e.g., [START_REF] Webster | Testing the accuracy of methods for reconstructing ancestral states of continuous characters[END_REF][START_REF] Oakley | Maximum likelihood models of trait evolution[END_REF]. This clearly needs further study because if the relative poor performance of inferring ancestral states is confirmed, this would require a revision of some previous results and ideas.

It must be pointed out that the ancestral states were estimated assuming uniform priors on the root (i.e., it could a priori be in any state). This can be relaxed, for instance, by assuming that the prior probabilities of the root are equal to the observed proportions of the states [START_REF] Maddison | Confounding asymmetries in evolutionary diversification and character change[END_REF]). It will be interesting to examine whether this has an effect on ancestral state estimates.

Conclusion

To conclude, I would echo [START_REF] Maddison | Confounding asymmetries in evolutionary diversification and character change[END_REF] view that uncareful analyses of evolutionary data can lead to wrong conclusions. Though the modern approach to macroevolution has certainly some limitations, there are reasons to be optimistic.

Future data analyses should include several tools based on model fitting, hypothesis testing, model checking, and parameter estimation. There is clearly a need to study the statistical behavior of our models in a wide range of realistic situations. There is also space for many methodological developments. One issue that still remains challenging is how extinction can be considered efficiently in these methods.

Fig. 1 .

 1 Fig. 1. Distribution summaries of the estimates of the rates of character change and their standard-errors. The dots indicate the median, and the error-bars the first and third quartiles. The panels correspond to the three simulated scenarios.

Fig. 2 .

 2 Fig.2. For each simulated data set, the relative likelihood that the root was in state x = 0 was calculated under both models of character change. This series of histograms shows the distribution of these values for each combination of parameters (as rows) and each actual initial state of the root (as pairs of columns). A value close to zero implies that the root state was inferred to be x = 1, a value close to one implies it was x = 0, and a value close to 0.5 implies that the inference of the root state was indeterminate by this model.

Table 1 .

 1 Proportions of cases where the null hypothesis was rejected (200 replications).A star indicates where the tested null hypothesis was actually true. The columns labeled 'Both' give the proportions of cases where the likelihood-ratio test (LRT) and the Wald test agreed on rejecting the null hypothesis.

			Character change	Diversification
	p 1	λ 1	LRT Wald	Both	LRT Wald Both
	p 0	5λ 0	0.18* 0.265* 0.055*	0.415 0.53	0.415
	0.1p 0	λ 0	0.405 0.31	0.26	0.01* 0.01* 0.01*
	0.25p 0 2.5λ 0 0.175 0.235	0.1	0.13	0.23	0.13

Acknowledgments

I am grateful to Wayne Maddison for reading a previous version of this paper, to

Gemain Chen for clarifying some issues on the Wald test, and to three anonymous referees and Arne Mooers for their constructive comments on the text. This research was supported by a project Spirales from the Institut de Recherche pour le Développement (IRD).


Emmanuel Paradis
email: paradis@isem.univ-montp2.fr

STATISTICAL ANALYSIS OF DIVERSIFICATION WITH SPECIES TRAITS

Keywords: .-diversification, extinction, maximum likelihood, phylogeny, generalized linear models, speciation

published or not. The documents may come

Statistical analysis of diversification with species traits

Bharathan 1993; [START_REF] Sanderson | Shifts in diversification rate with the origin of angiosperms[END_REF]Paradis 1998a;[START_REF] Bokma | Testing for equal rates of cladogenesis in diverse taxa[END_REF].

Other methods aim to test whether some traits are significantly correlated with species diversification [START_REF] Slowinski | Testing whether certain traits have caused amplified diversification: an improved method based on a model of random speciation and extinction[END_REF][START_REF] Barraclough | Sister-group analysis in identifying correlates of diversification[END_REF].

This category of methods is of particular relevance since they allow test of hypotheses about the effects of some biological traits on species diversification.

Examples of such traits include suspected evolutionary key-innovations [START_REF] Hunter | Key innovations and the ecology of macroevolution[END_REF], body size [START_REF] Gittleman | Body size and species-richness in carnivores and primates[END_REF][START_REF] Bokma | A statistical test of unbiased evolution of body size in birds[END_REF]Orme et al. 2002a), or other traits (e.g., [START_REF] Cardillo | Geographic range size, life history and rates of diversification in Australian mammals[END_REF][START_REF] Fisher | Extrinsic versus intrinsic factors in the decline and extinction of Australian marsupials[END_REF]. However, these methods suffer from the fact that they consider the product of evolution (namely species richness) rather than the process of evolution per se (speciation and/or extinction). Even analyzing species richness in a phylogenetic framework is not straightforward because of the possibility to use various indices of species richness differences, and the difficulty in combining several traits in the analysis [START_REF] Isaac | Phylogenetically nested comparisons for testing correlates of species richness: a simulation study of continuous variables[END_REF].

In this paper, I present a new method to analyze whether some traits affect diversification rate using observations of these traits on recent species and the phylogenetic tree of these species. This method is based on a new model of diversification which can be called the 'Yule model with covariates'. The parameters of the model are estimated by maximum likelihood. Several traits, either continuous or discrete, can be analyzed simultaneously. The statistical significance of the effects of these traits can be tested with likelihood ratio tests.

Since the present method is parametric, it is possible to predict how speciation rate varies with respect to the variables that have been found to be significant.

I develop below the model, the estimation procedure, and the statistical tests. I then report the results from a simulation study to assess the statistical properties of the method. To illustrate the use of the method, I analyzed some data on Primates.

Methods

The Model

For simplicity, I will assume that the lineage diversify following a birth-only process (i.e. no extinction). In this model, each species has an instantaneous probability of splitting in two daughter-species, called the speciation rate in this paper, and denoted λ. Let us assume that the value of λ for species i (λ i ) is determined by a set of species characters. For convenience, I will assume a linear relationship that can be written as:

λ i = β 1 x i1 + β 2 x i2 + . . . + β p x ip + α, (1)
where x 1 , . . . , x p are the species characters under consideration (x i1 . . . , x ip are the values of these characters for species i), and β 1 , . . . , β p , α are coefficients.

Since λ i is a probability, it can take values between 0 and 1; it is thus preferable to transform the left-hand side term above so that it varies between -∞ and +∞. Various functions can be used to this aim: I chose a logit function which is widely used in statistical analyses (e.g. in logistic regression):

ln λ i 1 -λ i = β 1 x i1 + β 2 x i2 + . . . + β p x ip + α. (2)
It is convenient to rewrite this equation in matrix form:

ln λ i 1 -λ i = x T i β, (3)
where β is a vector composed with the p + 1 coefficients, and x T i is the transposed vector of the p characters for species i binded with a 1 for the intercept (α). Thus, giving some values of x and β, the speciation rate is given by the inverse logit function:

λ i = 1 1 + e -x T i β .
(4)

Note that there is no direct time dependence in the speciation rate in this model: if the species characters vary through time, then λ will also vary, but if the former do not, then the latter will be constant through time.

The species characters (or variables) can be continuous or discrete. In the latter case, a suitable numeric coding may be used as common practice in generalized linear modeling [START_REF] Mccullagh | Generalized linear models[END_REF]. For instance, a discrete variable with two states would be coded as a variable taking the values 0 or 1 (see Appendix). Interactions between discrete variables, and between a continuous one and a discrete one can also be included in the model.

Parameter Estimation

The problem is to estimate the vector of coefficients β. I will assume that the available data are a phylogenetic tree that includes all species and with branch lengths proportional to time, and the values of the p species characters for the N species and for the N -1 nodes of the tree (assuming that the latter is fully dichotomous). A series of speciation events can be inferred from the tree; it is also inferred that no speciation occurred between two such consecutive events.

Let us denote t 1 , . . . , t q the times defined by the speciation events, and n 1 , . . . , n q the number of species living after each of these events (n q = N ). The first speciation event is given by the root of the tree, so that t 1 is the age of the root, and n 1 = 2. I will admit that there may be ties in these times (i.e.

speciation events occurring at the same time), so the number of species may increase by more than 1 through time, and q may be less than N -1. The probabilities of these events are proportional to the number of species living at a given time and the corresponding value of speciation rates; if two speciation events occurred at the same time (i.e. a tie), then this is the product of the number of species and the corresponding values of the two speciation rates (and so on with three or more speciation events).

The probability of a species not having a speciation event between times t and t is [START_REF] Cox | Analysis of survival data[END_REF]:

exp - t t λ(τ )dτ , (5)
where λ(τ ) is a function describing how the speciation rate changes between times t and t . The probability of n species not having any speciation event between times t and t is (see [START_REF] Darwin | The behaviour of an estimator for a simple birth and death process[END_REF]:

exp - n i=1 t t λ i (τ )dτ . (6)
We can now write the likelihood of the data under the model described above:

L = exp -t 2 t 1 λ 1 (τ )dτ + t 2 t 1 λ 2 (τ )dτ n 2 ν 2 × exp -t 3 t 2 λ 1 (τ )dτ + t 3 t 2 λ 2 (τ )dτ + t 3 t 2 λ 3 (τ )dτ n 3 ν 3 × . . . × exp - tq t q-1 λ 1 (τ )dτ + . . . + tq t q-1 λ N (τ )dτ ,
where ν 2 , . . . , ν q are the probabilities of the events observed at times t 2 , . . . , t q (namely one or more speciation events). A logarithmic transformation gives:

ln L = -t 2 t 1 λ 1 (τ )dτ + t 2 t 1 λ 2 (τ )dτ + ln n 2 + ln ν 2 + -t 3 t 2 λ 1 (τ )dτ + t 3 t 2 λ 2 (τ )dτ + t 3 t 2 λ 3 (τ )dτ + ln n 3 + ln ν 3 + . . . + - tq t q-1 λ 1 (τ )dτ + . . . + tq t q-1 λ N (τ )dτ .
Using the additive property of integrals, it is easy to see that all integrals add up over all branches. The log-likelihood thus becomes:

ln L = - 2N -2 j=1 t j t j λ(τ )dτ + q i=1 ln n i + q i=1 ln ν i , (7)
where t j and t j are the end-point times of the jth branch (2N -2 is the number of branch in a rooted dichotomous tree with N tips). A further simplification is

given by the fact that ν is a product of λ values, so we have:

ln L = - 2N -2 j=1 t j t j λ(τ )dτ + q i=1 ln n i + N -1 k=1 ln λ k . (8)
The way λ(τ ) is integrated over the branch lengths will depend on how the speciation rate varies through time. The simplest way is to assume a linear change:

t j t j λ(τ )dτ = λ t j + λ t j 2 l j , (9)
where l j is the length of the jth branch. Under this assumption, we end up with the log-likelihood:

ln L = - 2N -2 j=1 λ t j + λ t j 2 l j + q i=1 ln n i + N -1 k=1 ln λ k , (10)
which can be solved after substituting λ by the expression in equation 4.

Practically, maximizing ln L would require to find its first partial derivatives with respect to the different parameters; however these expression are too complex to be solved analytically. I used instead a numerical minimization method. The likelihood was transformed as the deviance (equal to -2 ln L);

finding the maximum likelihood was equivalent to minimizing the deviance.

This was done with the nonlinear minimization function of R [START_REF] Ihaka | R: a language for data analysis and graphics[END_REF] which uses a Newton-type algorithm for unconstrained minimization [START_REF] Schnabel | A modular system of algorithms for unconstrained minimization[END_REF]. This method allows one to find the minimum of a function even when its partial and second derivatives are unknown. However, these derivatives can be computed numerically by this algorithm, allowing to calculate the standard-errors of the MLEs with:

SE( β) =   - ∂ 2 ln L ∂β 2 β   -1 2 . ( 11
)
It is interesting to note that if we fix β = 0, then the estimated speciation rate λ, after logit back-transformation of the estimated parameter α using equation 4, is equal to the Kendall-Moran estimate (see Nee 2001, for details), as can be found with the function yule in APE [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]). The maximum likelihood found by both fitting procedures is the same.

Hypothesis Testing

To test whether the species variables x have a significant effect on speciation rate, it is possible to use standard statistical methodologies. Two models can be compared with a likelihood ratio test provided they are nested (i.e. one model is a particular case of the other). The hypothesis then tested is that the additional parameter(s) in the second model is (are) significantly different from zero. This likelihood ratio test follows a χ 2 distribution with a number of degrees of freedom equal to the difference in numbers of parameters between both models.

Diagnostics of Fit

Whereas statistical tests, such as likelihood ratio tests, indicate whether a model better describes the data than another one, they do not inform us whether this particular model adequately fits the data in some absolute ways. Diagnostics of fit are statistical tools to examine whether the observed data are well predicted by a particular model. For instance in a linear regression, residuals (the difference between the observed and the predicted values of the response) may indicate whether some potential predictors have not been included in the model, or whether the relationship may be nonlinear [START_REF] Venables | Modern applied statistics with S[END_REF].

With the method presented in this paper, it is possible to assess whether the branch lengths observed in a tree are well predicted by the Yule model with covariates using equation 5. This equation gives the probability of no speciation event, and so corresponds to the probability of observing a particular branch in the reconctructed tree. Under the above assumption on character change (i.e. linear change), equation 9 can be used practically.

Simulation Analysis

Some simulations were run to assess some statistical properties of the present method. Specifically, four questions were addressed. What is the type I error rate of the method? What is the statistical power of the method? What is the precision of the parameter estimators? What is the statistical robustness of the method with respect to extinction?

The evolution of a clade was numerically simulated using a non-homogeneous birth-death process where the speciation rate was determined by

λ = 1 1 + e βx-α ,
where α and β were parameters, and x was a continuous variable which evolved following a randow walk process, and the extinction rate µ was constant throughout time and all lineages.

The simulations were started with a single species. At each time step, each species living in the clade had a probability equal to µ to die. If it survived, each species had then a probability given by λ to generate two daughter-species, otherwise it simply survived to the next time step. This evolution process was simulated during 100 time-steps.

The evolution of x was so that at each time-step x t+1 = x t + with ∼ N (0, σ 2 ). All values of were independent among all lineages and throughout time. Such a model may not describe properly the evolution of many biological characters, but it was chosen because it can be simulated easily, and it conforms to the assumption of linear change along a branch of the tree.

The simulation parameters were given the following values: α = -3, β = {0, -1, -2, -3, -4}, σ = {0.01, 0.02}, and µ = {0, 0.0005, 0.001}. For each combination of parameters the simulation was replicated 100 times. Note that when β = 0, the simulated process was a homogeneous birth-death process with λ = 1/(1 + e 3 ) = 0.047. Fig. 1 illustrates how the speciation rate varied with Fig. 1 respect to x for the values of β used in these simulations.

Since the simulations were fully stochastic, so were the number of species living at the end of the simulation. The simulations with 0, 1, or 2 species at the end of the simulation were not considered. In some replications, the number of species were very large (up to 6,807,203 in the present study) which raised some problems because of the data manipulation subsequently needed. I thus considered only the trees with at most 2000 species. When a tree did not meet the selection criterion (between 3 and 2000 species), the simulation was repeated with the same parameter values until 100 trees were obtained for further analyses.

At the end of each simulation, the tree (topology with branch lengths) considering only the species living at the end of the simulation was output together with the values of x at its tips and nodes. Note that in real applications, the values at the nodes must be reconstructed separately. These data were analyzed with the method described in this paper. The estimates β, α, their standard-errors, the P -value of the likelihood ratio test of the hypothesis β = 0, the number of species living at the end of the simulation, and the depth of the tree (age of the most recent common ancestor of all living species) were stored for further analyses.

All analyses were done with R version 1.8.1 (R Development Core Team 2003). The simulations were done with a program written in C specially for the present study which was dynamically linked to R. The data analyses were performed with APE, a package written in R for phylogenetic analyses [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]).

Primates Data

A complete phylogeny of Primates has been recently reconstructed with the method of supertrees [START_REF] Purvis | A composite estimate of primate phylogeny[END_REF]. The phylogeny is ultrametric, and has branch lengths in million years (Fig. 2). The tree had several multichotomies Fig. 2 (191 nodes for 259 tips): it was necessary to accomodate this since the model used here considers only dichotomous trees. The simplest solution to this problem, and the one adopted here, is to treat these multichotomies as a series of dichotomies with zero-lengthed branches [START_REF] Purvis | Polytomies in comparative analyses of continuous characters[END_REF].

It was interesting to assess whether some biological traits could explain the variation in speciation rate in this order. An obvious candidate was body mass which relationship with diversification rate of Primates has been studied with another approach [START_REF] Gittleman | Body size and species-richness in carnivores and primates[END_REF]. Furthermore, one of the critical parameters that control speciation is the rate of fixation of new genetic mutations, and this rate is expected to be critically influenced by demography [START_REF] Gavrilets | Dynamics of speciation and diversification in a metapopulation[END_REF]. Small-bodied species have generally demographic features that make them likely to have a higher rate of fixation than large-bodied ones, such as largely fluctuating numbers [START_REF] Krebs | Population cycles in small mammals[END_REF], short dispersal distances [START_REF] Paradis | Patterns of natal and breeding dispersal in birds[END_REF], or short generation times [START_REF] Peters | The ecological implications of body size[END_REF]. It is therefore logical to expect a negative effect of body size on speciation rate. Data on body mass were taken from [START_REF] Smith | Body mass of late quaternary mammals[END_REF]. For the species with no data on body mass in this reference, the mean of the other species belonging to the same genus was calculated: this was done for 43 species. Body mass was still missing for three species which genus is monospecific (Cebuella pygmaea, Simias concolor, and Homo sapiens): data were found for these species in [START_REF] Nowak | Walker's mammals of the world[END_REF].

The values of body mass at the nodes of the phylogenies were estimated assuming a Brownian model of trait evolution [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF]. Under this model of trait evolution, the assumption in equation 9 is valid.

Results

Simulation Analysis

In eight cases out of the 3000 simulations the fitting algorithm did not converge which was clearly due to a small sample size since there were three species in these eight simulated trees.

The overall type I error rate at the nominal level of 5% calculated over all simulations with β = 0 was 0.022. For each combination of σ and µ, the estimated type I error rates ranged between 0.01 and 0.06. Given that these estimates are based on only 100 replications, there is no evidence that the probability of rejecting the null hypothesis when it is true overrates the nominal level of 5%, or that it is substantially affected by the extinction rate being different from zero.

The estimated power increased with greater values of β, larger values of σ, and smaller values of µ (Table 1). When β was large, the effect of the species Table 1 trait on speciation rate was stronger, and the likelihood ratio test detected this effect with a higher frequency. When σ was high the variation in the species trait was larger, and thus the relation between this trait and speciation rate was tighter and found statistically significant with a higher proportion. When µ > 0, the assumptions of the method were not fulfilled, and it may be expected that its statistical properties are affected. This was the case here, but this result depended on the values of the other parameters of the simulations. For instance, the likelihood ratio test was less powerful when β = -1 and µ = 0, than when β = -4 and µ = 0.001, even though the assumptions of the method were met in the former case.

An effect of sample size (number of species) was observed only when µ = 0 (Fig. 3). As could have been expected, the power of the test was higher for the Fig. 3 largest sample sizes. The relation was loose for β = -1, whereas the power was high (≈ 1) for β ≤ -2 when there was 150 or more species. It cannot be excluded that a relationship between power and sample size exists beyond the limit fixed in this study (2000) in the cases with µ > 0.

When considering the precision of the estimators, the results were quite sensitive to the depth of the tree: trees with a most recent common ancestor younger than 60 gave highly dispersed estimates of both parameters. Thus only trees with a depth of at least 60 were considered in the present analysis. The estimates of β were nearly unbiased when µ = 0, whereas a bias was observed when µ > 0, particularly when σ = 0.01 (Fig. 4). This bias was stronger for the Fig. 4 larger values of β. On the other hand, the estimates of α were nearly unbiased in all situations and their distribution showed no variation with respect to β or σ. The median of this distribution was -3.04 (first and third quartiles: -3.17, and -2.92; 98% of the values were distributed between -4 and -2).

Primates

Fitting the Yule model (i.e. the null model where λ is constant) gave a deviance of 95.20, and an estimate λ = 0.146 (SE = 0.009). Fitting a model with an effect of ln(body mass) on λ resulted in a deviance of 76.51. The likelihood ratio test of the effect of ln(body mass) was statistically significant: χ 2 1 = 18.69, P < 0.0001. Surprisingly, the effect of ln(body mass) was positive: β = 0.21 (SE = 0.03).

Previous analyses showed a strong contrast between different clades of Primates in terms of diversification rates [START_REF] Purvis | Macroevolutionary inferences from primate phylogeny[END_REF]Paradis 1998a).

Particularly, Old World monkeys (Catarrhini) appeared to have diversified at a much higher rate than the other groups of Primates. Given that Old World monkeys have on average larger body masses than New World ones (Platyrrhini), the positive effect of body mass on speciation rate could be due to an artefact. To assess this hypothesis, I fitted a model with a 'clade' effect where clade was a categorical variable with five categories: Ape, Catarrhini, Platyrrhini, Strepsirhini (Malagasy lemurs), and Tarsius. The states of this variable on the basal nodes of the phylogeny were reconstructed using a parsimony criterion. The deviance of this model was 11.46, and the likelihood ratio test comparing it with the null (Yule) model was χ 2 4 = 83.74, P = 0. Thus, the improvement in the model fit was dramatic.

I further tested whether ln(body mass) had a possible effect on this 'clade' effect model. The deviance of this 'clade + ln(body mass)' model was 7.21. The likelihood ratio test comparing this model with the 'clade' model was slightly significant: χ 2 1 = 4.25, P = 0.039. Interestingly, the effect of ln(body mass) was now negative (Table 2). I finally selected this last model for parameter estimation. Using equation 4 it is possible to compute the predicted values of Table 2 speciation rates according to the selected model (Fig. 5). Considering the Fig. 5 observed distributions of body mass in each clade of Primates, the 'clade' effect now explains the overall positive relationship between body mass and speciation rate. However, within each clade there is a negative effect of body mass on speciation rate.

A plot of the diagnostics of model fit shows that the terminal branches were particularly well predicted by the selected model (Fig. 2). Long or basal branches were less well predicted, though most branches had a probability greater than or equal to 0.2.

Discussion

Recent theoretical works on speciation and diversification have been taken in several directions such as characterizing the factors at the origin of speciation [START_REF] Kaneko | Sympatric speciation: compliance with phenotype diversification from a single genotype[END_REF][START_REF] Gavrilets | Models of speciation: what have we learned in 40 years?[END_REF][START_REF] Hochberg | Socially mediated speciation[END_REF][START_REF] Streelman | The stages of vertebrate evolutionary radiation[END_REF], the patterns of phylogenetic diversity [START_REF] Aldous | Darwin's log: a toy model of speciation and extinction[END_REF][START_REF] Aldous | Stochastic models and descriptive statistics for phylogenetic trees, from Yule to today[END_REF][START_REF] Losos | Stumped by trees? A generalized null model for patterns of organismal diversity[END_REF], or the evolution of genomic complexity [START_REF] Lynch | The origins of genome complexity[END_REF]. With the increasing quantity of various kinds of data (DNA sequences, ecological and biological data, . . . ), it is crucial to have some statistical methods to test the above theories. The goal of the method presented in this paper is to contribute to fill this need.

There are two main assumptions in the present method: extinction has been absent, and the fully resolved (species-level) phylogeny and the values of the traits under consideration at the tips and the nodes of the phylogeny are all known. These assumtions are discussed in the next four paragraphs.

The model used here (Yule model with covariates) assumes there is no extinction. This assumption is unlikely to be true in many situations as clearly shown by the fossil record of many taxonomic groups. The reason for this assumption is that it makes possible the development of the likelihood function presented above. It cannot be excluded that a likelihood function may be derived taking extinction into account, but this still needs to be found.

Interestingly, the simulation study showed that the type I error rate of the present method is not affected by random extinction: the probability of rejecting the null hypothesis when it was true was not increased when there were extinctions occurring at a constant rate. This suggests that the present method is unlikely to reveal false significant effect of a trait on diversification rate. However, we cannot generalize this result to the cases where extinction rate is heterogeneous. Particularly, it would be interesting to assess whether an extinction rate varying with respect to a species trait (e.g. in the same way modeled here for speciation rate) would be detected by the present method as a significant effect on speciation rate because of the 'signal' left on the tree by this heterogeneous extinction rate. This needs to be studied with more extensive simulations than done here.

The power of the test (probability of rejecting the null hypothesis when it is false) was affected by the value of extinction rate: the greater the latter, the smaller the former. However, it should be noted that this depended on the strength of the effect of the trait on the speciation rate: the test of the null hypothesis had a significant power when this effect was strong even though there were extinctions. It should be noted that the simulations were run with a relatively short timespan (100 time steps) which may explain the relative poor performance of the likelihood ratio test in this situation: it was observed that a significant effect was found mostly when β = -4. However, the analysis of the Primates data showed statistically significant estimated coefficients which were much smaller than four in absolute value (as the preliminary analyses of other data did). The fact that the depth of the tree was important for the results of the simulation suggests that the discreteness of the time scale may be critical here. Ideally, it would be needed to run some simulations with a more realistic time scale (> 1000) but these may be computationally hard to tract. Though the estimated power of the tests may not be generalized to real applications, it seems reasonable to assume that the patterns revealed by the simulations (increase in power with sample size, variation in the trait, and strength of the effect) may be generalized.

On a more general note, the estimation of extinction rates with phylogenies of recent species is a subject that should retain more attention. [START_REF] Kubo | Inferring the rates of branching and extinction from molecular phylogenies[END_REF] showed that the variance of the estimator of extinction rate obtained from molecular phylogenies is too large to be reliable. [START_REF] Paradis | Can extinction rates be estimated without fossils?[END_REF] showed that the estimates of extinction rates are negatively biased in a wide range of situations, and that correct estimation of this parameter is possible only when the reconstructed tree is close to the true historical tree of the clade. These results are intuitive since no extinction events are observed in a phylogeny of recent species. Most of the information in this kind of data relates to speciation rate. This justifies, at least partially, to focus the attention on speciation rate. I even suspect that it is not feasible to model extinction rates in the way done here for speciation rates. An additional consideration is that the extinct lineages may have parameters different from those of the extant ones such an a higher extinction rate. Thus, in addition to the problem that these lineages are not observed, it is possible that there is a problem of heterogeneity in the parameters as well.

The present method requires to know the fully resolved phylogeny of the species, and the values of the traits for these species at the tips and at the nodes of the phylogeny. It is generally not possible to have direct observations of the traits for the nodes of the trees, so they must be estimated using ancestral character reconstruction methods (Webster and Purvis 2002, and references therein). The assumption of a fully resolved tree can also be relaxed since a tree with multichotomies can be analyzed. Ideally, one needs to take such uncertainty into account in the estimation procedure so that the standard-errors of the parameters may be corrected. One possible approach is to repeat the analysis with different trees and different values of the ancestral traits, and then quantify the resulting variability on the parameter estimates as a component of their standard-errors due to this uncertainty [START_REF] Anderson | Suggestions for presenting the results of data analyses[END_REF]. Another, more sophisticated, approach is to use directly the confidence intervals on trees and ancestral trait values (see [START_REF] Schluter | Likelihood of ancestor states in adaptive radiation[END_REF] to compute a component of the variance of the estimators due to model selection uncertainty [START_REF] Buckland | Model selection: an integral part of inference[END_REF][START_REF] Burnham | Model selection and multimodel inference. A practical information-theoretic approach[END_REF].

The method in the form presented here considers only linear effects. Though this gives many possibilities such as including several variables and possible interactions between them, it is possible that nonlinear models may be more biologically realistic models of the effects of species traits on speciation rates. Such candidate models are threshold models where the effect of a variable differs depending on some threshold values. It is straightforward to extend the present method to nonlinear models: the only condition is that it must be possible to derive the value of speciation rate with respect to the variables included in the model (such as equation 4).

Throughout this paper, it has been assumed that continuous traits evolved in a simple way (Brownian motion), but another model of change may be used provided that the integration along each branch may be expressed in a way that it can be incorporated in the likelihood function. It should be noted that it is not required to know the exact temporal variation in species traits and speciation rate, but only the integral of the latter along each branch of the tree.

For instance, it does not matter whether the speciation rate changed linearly or with some random fluctuations along a branch, the integral will be the same in both cases.

The present method gives some emphasis on parameter estimation by contrast to hypothesis testing. With this respect, it is in agreement with some recent trends in statistical science and data analysis [START_REF] Johnson | Statistical sirens: the allure of nonparametrics[END_REF][START_REF] Johnson | The insignificance of statistical significance testing[END_REF][START_REF] Nelder | From statistics to statistical science (with discussion)[END_REF][START_REF] Anderson | Null hypothesis testing: problems, prevalence, and an alternative[END_REF][START_REF] Venables | Modern applied statistics with S[END_REF]. In evolutionary biology, null hypothesis testing is still extremely popular. I believe much could be gained from giving more emphasis on parameter estimation in this field.

This would allow the development of predictive models from the data (i.e. with an empirical ground). Null hypotheses usually focus on simple models whereas more complex models are often needed to describe biological processes. Another advantage of parameter estimation approaches is that the estimates (under some obvious conditions) are comparable among studies. Considering more carefully parameter estimation does not deny the importance of statistical tests as illustrated here where likelihood ratio tests were used to select models with the appropriate variables.

The results presented in this paper are an illustration of the potentialities of the present method, but they clearly need some comments. The data analyzed here were nearly the same than in [START_REF] Gittleman | Body size and species-richness in carnivores and primates[END_REF]. With respect to body size, the results from both studies were in the same direction: a positive effect on diversification was found for primates; however, the studies disagree with respect to the significance of this effect. This may come either from the different methods used, or from differences in the data. [START_REF] Gittleman | Body size and species-richness in carnivores and primates[END_REF] used the phylogenetically-based contrasts method [START_REF] Barraclough | Sister-group analysis in identifying correlates of diversification[END_REF] which is likely to explain the discrepancy with the present study.

An important issue raised by the Primates example is the treatment of multichotomies by the present method. The Yule model (with or without covariates) assumes that all speciation events are dichotomous, and thus that a phylogenetic tree under analysis should be fully dichotomous as well. However, many phylogenies at the species level have multichotomies. If the latter are the results of a series of speciation events so close in time that they cannot be resolved in the reconstructed phylogeny, then it is valid to treat them as a series of dichotomies with zero-lengthed branches. However, if the multichotomy is a consequence of a lack of knowledge (a so-called "soft polytomy"), then resolving the multichotomy in this way may be more problematic. The crucial point is that ancestral traits should be correctly estimated. Several studies suggest that unresolved multichotomies do not bias the estimation of ancestral traits [START_REF] Felsenstein | Phylogenies and the comparative method[END_REF][START_REF] Purvis | Polytomies in comparative analyses of continuous characters[END_REF][START_REF] Garland | Polytomies and phylogenetically independent contrasts: examination of the bounded degrees of freedom approach[END_REF].

On the other hand, errors in the topology of the tree is a more serious problem for the estimation of ancestral trait values [START_REF] Symonds | The effects of topological inaccuracy in evolutionary trees on the phylogenetic comparative method of independent contrasts[END_REF]. It seems that unresolved multichotomies will tend to hide actual relationships, and thus increase type II error rates, rather than increase type I error rates [START_REF] Garland | Polytomies and phylogenetically independent contrasts: examination of the bounded degrees of freedom approach[END_REF][START_REF] Symonds | The effects of topological inaccuracy in evolutionary trees on the phylogenetic comparative method of independent contrasts[END_REF]. This obviously needs further study.

It has been possible to show that the apparent overall positive effect of body mass on speciation rate was actually an artefact due to differences among the different clades of Primates. The most rapidly diversifying clades (Catarhini, Ape) have, on average, larger body sizes than the slowly diversifying clades (Tarsius, Strepsirhini): this created an apparent positive relationship between body size and speciation rate when the clade-effect was not taken into account. This among-clade heterogeneity in diversification has already been characterized in previous studies [START_REF] Purvis | Macroevolutionary inferences from primate phylogeny[END_REF]Paradis 1998a). Similar clade-specific diversification has been shown in Hawaiian birds [START_REF] Lovette | Clade-specific morphological diversification and adaptive radiation in Hawaiian songbirds[END_REF]). An interesting result from the present analysis is that, after taking into account this among-clade heterogeneity, the effect of body mass on speciation rate was negative, though only slightly statistically significant. It is possible that the effect of body mass, if any, on speciation rate is too weak in this mammal order to be characterized clearly, as was suggested by [START_REF] Gittleman | Body size and species-richness in carnivores and primates[END_REF].

Furthermore, recent assessments of the effect of body size on species diversity on a large scale have not given clear-cut results (Orme et al. 2002a,b).

To conclude, I believe the present method will bring the opportunity to assess many hypotheses about macroevolutionary processes. It will hopefully fill a gap between the approach of estimating speciation and extinction rates and In a model with a categorical variable (as color above) and a continuous variable (say x), it is possible to include, in addition to the main effects of both variables, an interaction term which is numerically coded as the product of the numeric variable coding for color and x:

red → 0 green → x
The model fit will give four parameter estimates: β1 (for x), β2 (for color), β3 (for the interaction), and α. The predicted values of speciation rate are then obtained with: red: logit(λ) = β1 x + α green: logit(λ) = ( β1 + β3 )x + β2 + α

This shows clearly that the interaction term is interpreted as a contrast in the effect of x with respect to the different categories of color (if β3 = 0, then only the intercept will be different). In the general case of a categorical variable with n categories, n -1 variables will be created with the product of the n -1 binary variables with x.

Similarly, a model with two categorical variables may have an interaction term numerically coded by all possible products of the individual numeric variables: a full model including the main and interaction effects of two categorical variables with n 1 and n 2 categories, respectively, will thus have n 1 n 2 -1 parameters (n 1 -1 for the main effect of the first variable, n 2 -1 for the second, and (n 1 -1)(n 2 -1) for the interaction). These values were used in the simulation study. The observed values of body mass for each species are indicated on the inner side of the x-axis.

Fig. 1 .

 1 Fig. 1. Examples of variation of speciation rate (λ) with respect to a trait (x).

Fig. 2 .

 2 Fig.2. Phylogeny of Primates from[START_REF] Purvis | A composite estimate of primate phylogeny[END_REF], and barplot of ln(body mass) for each species. The grey level indicates the probability that the given branch is observed according to the selected model.

Fig. 3 .

 3 Fig. 3. Relationship between the P-values of the likelihood ratio test of the hypothesis β = 0 and the number of species in the simulated trees for the different values of β and extinction rate (µ).

Fig. 4 .

 4 Fig. 4. Relationship between the true values of β and its estimates β for the different values of rate of trait evolution (σ) and extinction rate (µ). The dashed lines are x = y. Only the simulated trees with a depth of at least 60 were considered.

Fig. 5 .

 5 Fig. 5. Predicted values of speciation rate (λ) with respect to body mass for the different clades of primates according to the selected model. The bold parts of the curves show the range of observed values of body mass for each clade.

  Figure 2:

  Figure 4:

Table 1 :

 1 Estimated statistical power for each combination of strength of effect of species trait on speciation rate (β), trait evolution rate (σ), and extinction rate (µ).

	β	σ	µ
		0	0.0005 0.001
	-1 0.01 0.05 0.01	0.03
		0.02 0.15 0.05	0.01
	-2 0.01 0.17 0.03	0.02
		0.02 0.49 0.10	0.06
	-3 0.01 0.38 0.10	0.04
		0.02 0.68 0.26	0.12
	-4 0.01 0.58 0.05	0.09
		0.02 0.58 0.26	0.22

Table 2 :

 2 Parameter estimates for Primates.

	Effect		Estimate SE
	clade	Ape	-	-
		Catarrhini	1.22	0.24
		Platyrrhini	-0.36	0.27
		Strepsirhini	-0.86	0.30
		Tarsius	-1.91	0.46
	ln(body mass)		-0.17	0.06
	intercept		-0.32	0.60

Acknowledgments I am grateful to two anonymous referees and the Associate Editor for their constructive comments on a previous version of this paper. This research was financially supported by the Institut Français de la Biodiversité and the Centre National de la Recherche Scientifique. This is publication 2004-XXX of the Institut des Sciences de l' Évolution (Unité Mixte de Recherche 5554 du Centre National de la Recherche Scientifique).

Appendix

Some details are given here on how discrete (categorical) variables are handled in the present method, and how the corresponding parameters must be interpreted. For instance, consider color as such a variable which could take two values: red or green. To include this variable in a model, it will be substituted by a numeric in the following way: In the general case of a categorical variable with n categories, n -1 binary variables are created. For instance, with the primate data above, the four numeric variables were:

This led to the estimation of five parameters: β1 , β2 , β3 , β4 , and α. The predicted values of speciation rate were obtained with: