[R-sig-phylo] compar.ou in ape package - warnings and error messages on bird.orders example dataset
Emmanuel Paradis
Emmanuel.Paradis at mpl.ird.fr
Tue Mar 16 10:55:07 CET 2010
Annat (and Anna),
You're right about node numbers: I've fixed this problem and clarified
the help page. Now the node(s) can be specified with their labels as
well (if the tree has node labels). Besides, I fixed a bug in the
calculation of the deviance so the results should be more meaningful.
And the warnings should be avoided too (so transformations should not be
needed anymore). The fixed version is here:
https://svn.mpl.ird.fr/ape/dev/ape/R/compar.ou.R
and the help page:
https://svn.mpl.ird.fr/ape/dev/ape/man/compar.ou.Rd
Thanks for the helpful reports.
Emmanuel
Annat Haber wrote on 10/03/2010 13:55:
> Ok, I just realized that (obviously) the node number can't be smaller than
> length(bt) b/c n is always t-1 or smaller, if t is the number of tips and n
> is number of nodes, and nodes are numbered in the edge matrix as
> (t+1):(t+n). so there's no risk of getting erroneous results, it just
> wouldn't work.
> Of course, it would work fine regardless of node labels if you interpret
> "node number" to mean its position in the node vector, rather than its
> designated number in the edge matrix. The example in the help file uses its
> edge designated number (in what seems to be the old version of phylo),
> however, which is why this is confusing.
>
> ~ Annat
>
> On 7/03/10 9:12 PM, "Annat Haber" <annat22 at gmail.com> wrote:
>
>> 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
>
> _______________________________________________
> R-sig-phylo mailing list
> R-sig-phylo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>
--
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