[R-sig-phylo] problems to initialize.corPhyl for ancestal characters estimation with ace(ape) under a gls model

Emmanuel Paradis Emmanuel.Paradis at mpl.ird.fr
Tue Jan 20 09:05:23 CET 2009


It's a bug in ace(). The fix is to change this line:

o <- gls(x ~ 1, correlation = Initialize(corStruct, data.frame(x)))

with this one:

o <- gls(x ~ 1, data.frame(x), correlation = corStruct)

This will be in the next release of ape (in the next few days).

Thanks for pointing this out.

EP

Le 17.01.2009 00:49, Sebastien Lavergne a écrit :
> 
> I am bumping into a puzzling warning message when trying to compute 
> ancestral characters under a gls model, using the function ace (from 
> library ape). It seems that when initializing the chosen correlation 
> structure, the function does not get to match tip-names and rownames of 
> the tips data set.
> However corStruct does not cause any problem of this type when used to 
> fit a gls regression model.
> See the example below.
> I hope I am not being too naive and that the answer is not 
> overwhelmingly simple.
> 
> Any help appreciated on this too.
> Regards
> Seb
> 
> 
> Here is the example to run through:
> 
>  > library(ape) ;  library(geiger) ; library(picante)
>  > tree <- rcoal(50)
>  > tip.data <- evolve.brownian(tree, value=2, var=0.5)
> 
> Tip names and data rownames match perfectly
>  > name.check(tree, tip.data)
> [1] "OK"
> 
> However, when trying to compute ancestral character values with a gls 
> model in ace, here is the outcome:
>  > glsnodes <- 
> ace(tip.data,tree,type="continuous",method="GLS",corStruct=corBrownian(phy=tree))$ace 
> 
> Warning message:
> In Initialize.corPhyl(X[[1L]], ...) :
>  Row names in dataframe do not match tree tip names. data taken to be in 
> the same order as in tree.
> 
> Interestingly, this does not happen when using the pic method
>  > picnodes <- ace(tip.data,tree,type="continuous",method="pic")$ace
> 
> So let's see if corStruct works well in another gls fit
>  > tip.data <- data.frame(cbind(tip.data), rnorm(50))
>  > colnames(tip.data) <- c("depVar", "varExpl")
>  > glsModel1 <- gls(depVar ~ varExpl , data=tip.data, method="ML", 
> correlation=corBrownian(phy=tree))
>  > glsModel1
> Generalized least squares fit by maximum likelihood
>  Model: depVar ~ varExpl
>  Data: tip.data
>  Log-likelihood: 31.67452
> Coefficients:
> (Intercept)     varExpl
> 1.955331756 0.001608361
> 
> If the data are randomized, the fitted model is strictly the same, 
> meaning that tips and rownames are match perfectly
>  > tip.data <- tip.data[sample(rownames(tip.data), 
> length(rownames(tip.data))),]
>  > glsModel2 <- gls(depVar ~ varExpl , data=tip.data, method="ML", 
> correlation=corBrownian(phy=tree))
>  > glsModel2
> Generalized least squares fit by maximum likelihood
>  Model: depVar ~ varExpl
>  Data: tip.data
>  Log-likelihood: 31.67452
> Coefficients:
> (Intercept)     varExpl
> 1.955331756 0.001608361
> 
> 

-- 
Emmanuel Paradis
IRD, Montpellier, France
   ph: +33 (0)4 67 16 64 47
  fax: +33 (0)4 67 16 64 40
http://ape.mpl.ird.fr/



More information about the R-sig-phylo mailing list