Re: [R-sig-phylo] ape - corMartins

2017-01-25 Thread Sergio Ferreira Cardoso
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

2017-01-24 Thread Emmanuel Paradis

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/