2

I am trying to determine which evolutionary model is best for my discrete data using the function fitDiscrete() from the geiger package.

These are the values that I get for the number of parameters (k), maximum log likelihood (lnL), AIC and AICc for each model:

     k   lnL        AIC       AICc
ER   1   -115.8006  233.6012  233.6637
ARD  90  -85.98459  351.9692  -303.2308
SYM  45  -97.23202  284.4640  491.464

The same dataset (n = 66), tree and single trait with 10 levels were used to create each model. The only difference is the evolutionary model fitted (equal rates (ER), all rates different (ARD) and symmetrical rates (SYM))

I am having trouble interpreting these results, however.

To start, for AIC, I'm fairly sure that I should select the model with the smallest AIC score, i.e the ER model.

For lnL, however, I have seen that the model with the "largest value" should be selected with this being interpreted as the value closest to 0 (https://www.r-phylo.org/wiki/HowTo/Ancestral_State_Reconstruction), i.e. the ARD model. I realise though that lnL values tend to be biased towards models with higher k values. To address this, I did do a likelihood test as suggested by the website above (chi-squared test), which came to p < 0.001. This would suggest that the ARD model should be preferred over the ER model, which contradicts what the AIC scores are telling me.

As for AICc, again, the "smallest value" should be selected but the negative sign mixed in with the positive ones has thrown me. Is this the smallest absolute value or the value value closest to negative infinity?

So, all in all, how can I tell which model should be preferred?

kjetil b halvorsen
  • 63,378
  • 26
  • 142
  • 467
Carolina Karoullas
  • 141
  • 3
  • How many degrees of freedom did you use in your chi-square test? With 89 degrees of freedom, the difference between the ER and ARD log likelihoods doesn't look significant to me... And is the ER model nested in the ARD model? – jbowman Nov 12 '20 at 19:18
  • I've had to double check the website above but the df should be the difference between the values for k so 89 here. The likelihood test is calculated through the following command (from the webiste above): 1-pchisq(2*abs(ER(lnL) - ARD(lnL)), 1). ER should be nested in ARD but the code isn't clear: ER – Carolina Karoullas Nov 12 '20 at 20:09
  • The degrees of freedom of the chi-square test is the difference in the dimensionality of the parameters, in this case, 89: https://en.wikipedia.org/wiki/Wilks%27_theorem , the example should make it clear. This should help reconcile the difference you are seeing! ... so `pchisq(2*abs(ER(lnL) - ARD(lnL)), 89)`. – jbowman Nov 12 '20 at 20:26
  • Ahh that would be it! I've just tried it and it all works now! I knew the df should be 89 but the example didn't specify where to insert that and I didn't realise that their model only had 1 df! Thank you :) – Carolina Karoullas Nov 12 '20 at 22:23
  • The AICc value for ARD cannot be correct since AICc is always larger than the corresponding AIC. The AICc adds an additional penalty for parameter number. See here: https://en.wikipedia.org/wiki/Akaike_information_criterion#Modification_for_small_sample_size – abstrusiosity Nov 12 '20 at 23:07
  • Hi that's what I thought too so I calculated it by hand following the formula on the same wikipedia link as you sent and I arrived at the same result as what R got: 351.9692+(2*90*90+2*90)/(66-90-1) = -303.2308. The number of parameters (k = 90) exceeds the number of observations (n = 66) hence the negative result. – Carolina Karoullas Nov 13 '20 at 14:31

0 Answers0