One-way ANOVA with phylolm

[matthieu-haudiquet]

Hello,

Thanks again for the very useful phylolm() function :-).

I would like to perform a one-way ANOVA on the output of phylolm. Do you know if this is possible ?

Best,
Matthieu

[cecileane]

If you use formula y ~ x in phylolm(y ~ x, ...), where y is your response variable and x is one categorical predictor (with 2 or more levels), then you do a phylogenetic one-way ANOVA. The standard one-way ANOVA would be using lm(y ~ x). Or did I misunderstand your question?

[matthieu-haudiquet]

I'm sorry, I wasnt clear .
I am wondering if there's a straightforward way to get the same summary as the anova() function.
With lm, that would be :
mod <- lm(y ~ x)
anova(mod)

And not

mod <- lm(y ~ x)
summary(mod)

Best,
Matthieu

[cecileane]

I see. No, anova and drop1 don't work on phylolm objects. But you can use phylostep to compare models, and update also.

I don't recommend anova generally, because the p-values that it returns are sequential. For instance, the first test compares a model with an intercept only to a model with the first predictor only. Usually, that's not what we want to do (because ignoring the other predictors is typically bad), but it's easy to forget. I recommend drop1 instead. Its first p-value would compare a model with all predictors but the first one, with a model with all predictors (including the first one).

We can do these model comparisons with phylostep (which focuses on AIC) or manually with update and extracting log-likelihoods, like this for instance:

fit1 = phylolm(y ~ x1 + x2 + x3, phy=..., ...)
fit0 = update(fit1, . ~ . - x1) # model without x1 in the list of predictors
AIC(fit0) - AIC(fit1)
# for a likelihood ratio test:
x2 = 2*(logLik(fit1)$logLik - logLik(fit0)$logLik)
df = logLik(fit1)$df - logLik(fit0)$df
pchisq(x2, df, lower.tail=F) # approximate p-value