# Model selection in PGLS?

I am using Phylogenetic General Least Squares in the R package 'caper'.

I have 4 predictor variables and I would like to know which are correlated with my response variable, while taking the phylogeny into account.

I am uncertain on how to assess the models:

1. one way is to include all my predictors in one model and get the F-stat using anova().
2. The other option is to assess multiple models with different combinations of the predictor variables and use AIC() to select the best model.

Which method should I use?

If your goal is to obtain p-values for the individual predictors, then using anova on your fitted pgls models does that. AIC will give you measures of the overall goodness of fit of each model, but doesn't tell you about the significance of individual parameters. In fact, you can have a model favored by AIC that includes parameters with non-significant p-values (see Why applying model selection using AIC gives me non-significant p-values for the variables.

On the other hand, it is arguably more useful to consider which variables are included in the model selected using AIC, because those variables with non-significant p-values may still be useful. In particular, when you have collinear predictor variables, the p-values can be large even for variables that add predictive power to the model (see AIC or p-value: which one to choose for model selection?).

In practice, I would suggest reporting results from both approaches. While I always appreciate seeing the different models along with AICs, many readers and reviewers care more about seeing p-values for individual predictors. It is easy enough to report AICs along with p-values for predictors, and you may get insights about your data from looking at the AICs and p-values together (e.g., if you have a model with significant variables, and then adding a new variable improves the AIC but causes the other model variables to become non-significant, it would suggest that the new variable is correlated with the old variables and you should not necessarily conclude that the old variables don't matter, even if they have non-significant p-values in the full model).

• Thanks for the clear answer. So should i use anova() for the model with all parameters or the best model selected by AIC? (I notice the output changes with spurious parameters). Cheers – AlexR May 16 '16 at 6:13
• Well what I was suggesting was to provide the p-values for all of the models/parameters, as well as the AIC. That way people can see how p-values change depending on the model specification. What do you mean "spurious parameters"? – Slow loris May 16 '16 at 15:26
• Oh ok, but reporting the p-values for all the parameters/models is a bit tedious if you are modelling many parameter combinations? "spurious" in reference to the parameters excluded in best model determined by AIC. – AlexR May 19 '16 at 7:51
• If your concern is reporting a ridiculous number of results, then you can choose to only report models above a certain Akaike weight. If your concern is doing everything by hand, then you might find help from an R package such as MuMln. You can also look into model averaging rather than model selection. – Slow loris May 19 '16 at 16:20