Re: [R-sig-phylo] ape - corMartins
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 BIClogLik fitOU0.11 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 BIClogLik stepwise0.1 1 6 40.14718 51.37439 -14.07359 stepwise12 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.3378136590.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.7442470550.291622930 loc_typSemiaquatic loc_typTerrestrial 0.0292194950.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 BIClogLik fitOU0.11 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 fitPagel5 18 42.54809 76.22971 -3.274044 Model df AIC BIClogLik stepwise0.1 1 6 40.14718 51.37439 -14.07359 stepwise12 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 stepwisePagel5 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
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/
[R-sig-phylo] ape - corMartins
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. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology Museu da Lourinhã, GEAL. LATR/IST/CTN - Campus Tecnológico e Nuclear. Lisboa, Portugal [[alternative HTML version deleted]] ___ 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/
