[R-sig-phylo] compar.ou in ape package - warnings and error messages on bird.orders example dataset

Annat Haber annat22 at gmail.com
Mon Mar 8 04:12:36 CET 2010


Hi Anna and everyone,
I'd like to add to Anna's post, now that I had to figure it out for my own
stuff...
Regarding the colnames(Wend) error message:
The function documentation says that "node" needs to be specified as node
number, but in fact it needs to be specified as node name if there are node
names (that are not just as.character of their numbers).
Unless I'm missing something, the problem seems to be with these lines (the
first and the last)
    bt <- branching.times(phy)
    Tmax <- bt[1]
    Wend <- matrix(0, nb.tip, length(node) + 1)
    colnames(Wend) <- c(names(sort(bt[node])), as.character(root))

The thing is that branching.time doesn't remember the original node numbers,
but their names. In the bird.orders example the node numbers are also the
node names, which is why as.character(node) and specifying node=25 was doing
the trick for that particular example. But if the nodes have other names
that wouldn't work, unless you specify in the function call their names
rather than their numbers. If the node number happened to be smaller than
length(bt) then it wouldn't give an error message but it would give
erroneous results. It would fit the wrong model, obviously.

As for the warnings:
They seem to be related to data location and scale. My data happened to be
PC scores so obviously all colMeans were zero. Adding 1000 to everything
helped. Well, to the extent that 18 warnings are better than over 50... and
getting 1 standard error NA rather than all of them. I had to play with the
number a bit. 1000 was giving me " Error in solve.default(out$hessian) :
  system is computationally singular: reciprocal condition number =
4.63219e-23" and 10 was giving 2 standrad errors NA and a lot more warnings.
Scaling didn't help me here at all, but it did in other examples, like when
dealing with procrustes shape data, so I thought I'd mention it.

HTH,
and would greatly appreciate any input,
Annat


On 2/03/10 9:33 AM, "Anna Kostikova" <anna.kostikova at gmail.com> wrote:

> Just to follow up the warnings and errors on compar.ou function, I've
> tried to debug a bit a code, and found out that at least for the error
> message (where node = -2 is specified), the problem seems to be
> connected with the way function compar.ou process this bit of code:
> 
>> colnames(Wend) <- c(names(sort(bt[node])), as.character(root))
> 
> I am not sure why in the documentation example the node is set to -2,
> but anyway, even if set to a particular node number (for bird.orders
> corresponding node would be 25), the code still collapse at this
> point, due to the fact that "node" should be a character, not numeric
> vector.
> So, it might be worth adding one line of code before, as to get a
> correct result:
> 
>> node <- as.character(node)
>> colnames(Wend) <- c(names(sort(bt[node])), as.character(root))
> 
> Please, correct me if I've got it wrong and I would be very glad to
> hear opinions on that and also follow up on warnings (which are
> connected with the way minimization is done, but I am not clear what
> exactly goes wrong there).
> 
> Thanks a lot again,
> 
> Anna
> 
> 2010/3/2 Anna Kostikova <anna.kostikova at gmail.com>:
>> Dear list,
>> 
>> When running compar.ou example from APE package documentation (v.
>> 23-11-09), I've got couple of error messages.
>> What might be the reason for them:
>> 
>>> compar.ou(rnorm(23), bird.orders, alpha = 0.1)
>> $deviance
>> [1] 57.60627
>> 
>> $para
>>          estimate    stderr
>> sigma2  0.71937096 0.1500408
>> theta1 -0.04114436 0.1365702
>> 
>> $call
>> compar.ou(x = rnorm(23), phy = bird.orders, alpha = 0.1)
>> 
>> attr(,"class")
>> [1] "compar.ou"
>> Warning messages:
>> 1: In log(2 * pi * p[2]) : NaNs produced
>> 2: In nlm(function(p) dev(c(alpha, p)), p = c(1, rep(mean(x), ncol(Wstart))),
>>  :
>>  NA/Inf replaced by maximum positive value
>> 
>>> y <- c(rnorm(5, 0), rnorm(18, 5))
>>> compar.ou(y, bird.orders, node = -2, alpha = .1)
>> Error in dimnames(x) <- dn :
>>  length of 'dimnames' [2] not equal to array extent
>> 
>> Thanks a lot for any clue,
>> 
>> Best regards,
>> Anna
>> 
> 
> _______________________________________________
> R-sig-phylo mailing list
> R-sig-phylo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo



More information about the R-sig-phylo mailing list