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/
