significance in phylolm analysis

[paternogbc]

Dear @lamho86,

This might be a sily question, however I am a bit confused. How can I compute a ANOVA table of a model fitted with 'phylolm()' function?

I am perfoming some analysis and want to check the significance of a explanatory variable. Lets say:

mod <- phylolm( Y ~ X1 + X2, data, phy, model = "lambda"). summary(mod)

dont show a significant effect of the variable X2. However if I use gls approach (with corPagel) or the function pgls from the package caper, the variable is significant when looking in the ANOVA table.

And both methods (gls and caper) gives identical summary(mod) results provided by 'phylolm'.

Could you please give me some advice on that? Should I look for significance of a variable with summary or ANOVA?

Kind regards, Gustavo

[cecileane]

In R, the anova() methods typically use type 1 sums of squares, that is, they provide sequential tests: they first give the significance of the first predictor ignoring all others, then they give the significance of the second predictor after accounting for any variation explained by the first (but ignoring all others), etc. So the results of anova() can depend a lot on the order in which you specify the variables: the significance of X2 can be different if you use the formula Y~X1+X2 or if you use Y~X2+X1.

In contrast, the summary() methods (including that in phylolm) uses type III sums of squares, that is, they give the significance of each coefficient after accounting for the effect of all other coefficients. If predictors are correlated, you might get that none of them are significant (after accounting for the effect of the others). The p-values given by the summary methods are not affected by the order in which you give the predictors, in your formula.

In your example, did you check the results for X2 when you changed the order of predictors in the formula? If the results change, the difference in p-values might be coming from this difference of using type 1 versus type 3 sums of squares.

[paternogbc]

Dear @cecileane,

thanks a lot for you kind explanation. The F values do change if the order of variables change in formula (due to the Type I sums of squares as you explained). But still significant.

The problem is that my X2 variable is categorical (with two levels), so in summary the intercept is level A and the second level is presented below. I have a highly significant p-value for the intercept and a non-significant p-value for the second level.

I am interested in testing the significance of the the variable itself (X2) so which p-value should I consider? As I understand, these p-values in summary are testing if the levels have an effect size different then zero. Is that correct?

How can I be sure that the categorical variable X2 have a significant effect in Y?

I have also performed LRT between models with and without X2 and removing this variable causes a significant change in the model. But I am not sure which is the best way to report these p-values.

Sorry to load you with these questions, but since we are working on a package that uses phylolm, due to its much faster math, I am concern about which p-values should I use when categorical variables are used in the model.

Thanks once again.

[cecileane]

Do you have one or more interaction terms involving this predictor X2? Also, do you get identical results when you do summary() either with phylolm or with caper?

[paternogbc]

Dear @cecileane,

I have one interaction term X1:X2 where X1 is numeric and X2 is a factor with two levels. Yes the output from summary is identical with phylolm, caper or gls using model = "lambda".

[cecileane]

Hi Gustavo, The interaction is the cause of the discrepancy: the p-values from anova() and from summary() are not testing the same thing: one will assume no interaction, the other will account for an interaction.

• anova: the p-values here are for sequential tests, and (by default) we never include an interaction term without all the main terms in that interaction. So if your formula is Y ~ X1*X2, the anova table will first give you a p-value for X1 ignoring any effect of X2 and ignoring the interaction X1:X2. Then it will give your a p-value for X2 after accounting for any effect of X1, but assuming no interaction X1:X2. And finally, it will give you a p-value for the interaction X1:X2 after accounting for the main effects of X1 and of X2. Roughly, the test of X2 here assumes no interaction, so it roughly tests the average effect of X2, over all values of X1.
• summary: in contrast, with summary you will get a p-value to test the hypothesis that beta_i=0 for each coefficient beta_i, after accounting for all other terms and coefficients. So it will test a particular coefficient for X2 accounting for the interaction X1:X2 (anova would assume no interaction instead). That is probably the reason why your two p-values for X2 are not the same: one assumes no interaction, and the other is for the effect of X2 at a particular level of X1 (at X1=0 if X1 is continuous), assuming an interaction (that is: the effect of X2 depends on the value of X1).

Both p-values can be of interest! Which one should be used depends on the question.

[paternogbc]

Dear @cecileane,

thank you so much for this explanation. This makes much more clear why I am finding discrepancies between p-values.

This solved my problem!

Cheers

[chrisjlaw]

Hi Cécile,

I have a similar question to this thread. I ran a model with an interaction (Y ~ size*ecotype) and have been calculating the slopes and using their 95% bootstrap CI to see if they differ from 0. However, I'm still confused about the P values and was wondering if you could help me interpret what they mean in this example:

Thanks!

[cecileane]

The p-values have their usual meaning, the same as p-values from regular lm models in R. Many textbooks explain the interpretation of interaction terms and p-values.

Here, these p-values are conditional on the estimated phylogenetic covariance (e.g. conditional on the estimated λ value, if you used Pagel's λ model). In particular, p-values are calculated from the t-values, not from the bootstrap procedure.

For example, the p-value 0.8727 for lnhumerus_size corresponds to the two-tailed probability beyond the t-value -0.1608897, under a t-distribution: 2*pt(-0.1608897, df=60) assuming a residual degree of freedom df=60, say. This p-value tests the null hypothesis that the slope for lnhumerus_size is 0, assuming that all other coefficients are included in the model --including all coefficients for the lnhumerus_size:ecotype interaction term. That's not a good hypothesis to test, because the inclusion of this interaction means that lnhumerus_size is still part of your model anyway.