Hi all.
If you don't want to apply the fix, you can circumvent the problem by
first assuring that meserr is in the same order as the tree tip labels,
and then removing the names from the input vector of standard errors.
So, for instance:
se<-se[tree$tip.label]
names(se)<-NULL
fitContinuous(tree,x,model="lambda",meserr=se)
works. Note that the order of meserr should be, so far as I can tell,
the order of tree$tip.label and not the order of the data in x. This is
because the internally called function treedata(...,sort=T) sorts the
data into the order of the tip labels.
This also agrees with Matt's diagnosis of the problem because if
is.null(names(meserr)) then the code with the bug is circumvented.
In theory, one should be able to supply meserr as a matrix;
unfortunately there is a second bug whereby:
o<-match(rownames(td$data),rownames(meserr))
meserr<-meserr[o,]
should be changed to
o<-match(rownames(td$data),rownames(meserr))
meserr<-as.matrix(meserr[o,])
because the former inadvertently converts meserr to a vector, when a
matrix is subsequently needed.
- Liam
--
Liam J. Revell
University of Massachusetts Boston
web: http://faculty.umb.edu/liam.revell/
email: liam.rev...@umb.edu
blog: http://phytools.blogspot.com
On 2/2/2012 1:41 PM, Matt Johnson wrote:
I'm not sure if this is the best place to make a bug report, but I thought
others might be interested.
I was trying to replicate, with geiger, a demonstration of estimating lambda
with measurement error from Liam's blog:
http://phytools.blogspot.com/2011/11/including-measurement-error-in.html
Replicating this with fitContinuous, I received the following error:
fitContinuous(tree,xe,model="lambda",meserr=se)
Error in meserr[o, ] : incorrect number of dimensions
The problem is related to this bit of code in fitContinuous, used to determine
whether meserr is a single value, a vector, or a matrix:
else if (is.vector(meserr)) {
if (!is.null(names(meserr))) {
o<- match(rownames(td$data), names(meserr))
if (length(o) != ntax)
stop("meserr is missing some taxa from the tree")
meserr<- as.matrix(meserr[o,])
Removing the comma from the final line fixed the problem. The values I receive
from geiger after removing the comma are identical to the values I get from
phytools::phylosig for the same data, so I am assuming I haven't broken
anything important.
The full code I used is at the end of this message, for those interested.
Best,
Matt Johnson
library(geiger)
Loading required package: ape
Loading required package: MASS
Loading required package: mvtnorm
Loading required package: msm
Loading required package: ouch
Loading required package: subplex
library(phytools)
Loading required package: mnormt
Loading required package: phangorn
Loading required package: quadprog
Loading required package: igraph
Loading required package: rgl
source("mjContinuous.R") #fitContinuous with offending comma removed
set.seed(2)
gen.lambda=0.7
tree = drop.tip(birthdeath.tree(b=1,d=0,taxa.stop=501),"501")
se = sqrt(rchisq(n=500,df=2))
names(se) = tree$tip.label
x= fastBM(lambdaTree(tree,gen.lambda),sig2=1)
e = rnorm(n=500,sd=se)
xe=x+e
fitContinuous(tree,x,model="lambda")
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -1031.783
$Trait1$beta
[1] 1.221674
$Trait1$lambda
[1] 0.7613218
$Trait1$aic
[1] 2069.567
$Trait1$aicc
[1] 2069.615
$Trait1$k
[1] 3
fitContinuous(tree,xe,model="lambda")
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -1154.987
$Trait1$beta
[1] 1.412728
$Trait1$lambda
[1] 0.5495113
$Trait1$aic
[1] 2315.975
$Trait1$aicc
[1] 2316.023
$Trait1$k
[1] 3
fitContinuous(tree,xe,model="lambda",meserr=se)
Error in meserr[o, ] : incorrect number of dimensions
mjContinuous(tree,xe,model="lambda",meserr=se) #with problematic comma removed
Fitting lambda model:
$Trait1
$Trait1$lnl
[1] -1141.14
$Trait1$beta
[1] 1.170196
$Trait1$lambda
[1] 0.7484273
$Trait1$aic
[1] 2288.281
$Trait1$aicc
[1] 2288.329
$Trait1$k
[1] 3
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
