Dear Emmanuel, Thank you for your answer. I did try your suggestion and the alpha estimate was 0.89. I created models with alpha= 0.1, 1, 50 and a Brownian Motion model, which looks like it is the best fitting model.
Model df AIC BIC logLik fitOU0.1 1 17 49.25612 81.06653 -7.628058 fitOU1 2 17 50.42329 82.23371 -8.211647 fitOU50 3 17 50.42330 82.23371 -8.211648 fitBrownian 4 17 42.82985 74.64027 -4.414926 When I do a stepwise regression, it's even more accentuated: > stepwise<-stepAIC(model,direction="both",trace=0) Model df AIC BIC logLik stepwise0.1 1 6 40.14718 51.37439 -14.07359 stepwise1 2 6 41.66723 52.89444 -14.83362 stepwise50 3 6 41.66724 52.89444 -14.83362 stepwiseBrownian 4 2 27.97681 31.71921 -11.98841 Also, the model fit does not increase with alpha >1. I tried the following: > fitPagel<-gls(fcl~mass+loc_dim+activity+feeding+loc_typ+agility,correlation=corPagel(value=1,phy=tree,fixed=FALSE),data=df,method="ML",weights=varFixed(~vf)) > fitPagel Generalized least squares fit by maximum likelihood Model: fcl ~ mass + loc_dim + activity + feeding + loc_typ + agility Data: df Log-likelihood: -3.274044 Coefficients: (Intercept) mass -0.337813659 0.058271471 loc_dim3D activityDiurnal/Nocturnal 0.519780792 -0.202720741 activityNocturnal feedingOccasional_predator 0.057773877 -0.181687109 feedingPredator loc_typFlyer -0.153543188 -0.007564622 loc_typFossorial loc_typScansorial 0.744247055 0.291622930 loc_typSemiaquatic loc_typTerrestrial 0.029219495 0.449648236 agilityMedium agilityMedium_fast -0.072365659 -0.246153152 agilityMedium_slow agilityVery_slow -0.042370670 -0.255200809 Correlation Structure: corPagel Formula: ~1 Parameter estimate(s): lambda 0.9388239 Variance function: Structure: fixed weights Formula: ~vf Degrees of freedom: 48 total; 32 residual Residual standard error: 0.02930323 #Comparison between models Model df AIC BIC logLik fitOU0.1 1 17 49.25612 81.06653 -7.628058 fitOU1 2 17 50.42329 82.23371 -8.211647 fitOU50 3 17 50.42330 82.23371 -8.211648 fitBrownian 4 17 42.82985 74.64027 -4.414926 fitPagel 5 18 42.54809 76.22971 -3.274044 Model df AIC BIC logLik stepwise0.1 1 6 40.14718 51.37439 -14.07359 stepwise1 2 6 41.66723 52.89444 -14.83362 stepwise50 3 6 41.66724 52.89444 -14.83362 stepwiseBrownian 4 2 27.97681 31.71921 -11.98841 stepwisePagel 5 4 28.65373 36.13854 -10.32687 Since the Pagel estimates lambda values very close to 1, I believe this means my traits are evolving accoring to a Brownian Motion model. I'll try to search in the literature for a discussion on the estimation of alpha. Thank you, once again. Best regards, Sérgio. ----- Mensagem original ----- > De: "Emmanuel Paradis" <emmanuel.para...@ird.fr> > Para: "Sergio Ferreira Cardoso" <sff.card...@campus.fct.unl.pt>, > "r-sig-phylo" <r-sig-phylo@r-project.org> > Enviadas: Terça-feira, 24 De Janeiro de 2017 19:21:46 > Assunto: Re: [R-sig-phylo] ape - corMartins > Hi Sérgio, > > It seems your results make good sense. alpha=0 is too far from the > "optimum" value (i.e., the value that maximises the log-lik), so the > optimisation fails to maximise the log-likelihood with respect to this > parameter. It has been written in the literature that it is difficult to > estimate reliably alpha in an OU model, and your results seem to confirm > this. > > Besides, alpha=0 is equivalent to a Brownian motion model. So if your > traits are not really evolving according to this model, the GLS fitting > procedure may be in difficulty. > > You may try: > > correlation = corMartins(value=0.9, phy=tree, fixed=FALSE) > > If this estimates alpha=0.99, then you may probably rely on this result > (after checking the estimates of the coefficients of course). > > Best, > > Emmanuel > > Le 20/01/2017 à 15:16, Sergio Ferreira Cardoso a écrit : >> Dear all, >> >> I tried to estimate a value for alpha parameter of corMartins >> (Ornstein-Uhlenbeck) but the output is a bit puzzling. I do as follows: >> >>> >> fitOU<-gls(fcl~mass+loc_dim+activity+feeding+loc_typ+agility,correlation=corMartins(value=0,phy=tree,fixed=FALSE),data=df,method="ML",weights=varFixed(~vf)) >>> Error in recalc.corStruct(object[[i]], conLin) : >> NA/NaN/Inf in foreign function call (arg 4) >> >> Enter a frame number, or 0 to exit >> >> 1: gls(fcl ~ mass + loc_dim + activity + feeding + loc_typ + agility, >> correlati >> 2: nlminb(c(coef(glsSt)), function(glsPars) -logLik(glsSt, glsPars), >> control = >> 3: objective(.par, ...) >> 4: logLik(glsSt, glsPars) >> 5: logLik.glsStruct(glsSt, glsPars) >> 6: recalc(object, conLin) >> 7: recalc.modelStruct(object, conLin) >> 8: recalc(object[[i]], conLin) >> 9: recalc.corStruct(object[[i]], conLin) >> >> It work when I turn the "value" to 1, but then the estimate for alpha is >> 0.99, which looks like it isn't actually estimating a value and is instead >> accepting my input as the alpha value. Do you think there is a reason for >> this to happen? >> >> Thank you in advance. >> >> Cheers, >> Sérgio. _______________________________________________ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/