Hello,
I'm a newbie in phylogenetic analysis so excuse me in advance for
incorrect terms and so on (or if my question is too silly).
I've recently started to work on a project for which I need to use
Ornstein-Uhlenbeck models and I decided to use R and ouch since
the documentation was great.
For the first experiments I have been working on some toy trees but
I'm getting puzzling results: the same tree, simply
defined with different orders of the string which I use to represent
it in the Newick format, when analyzed with the Hansen
functions gives different results.
I have tried different combinations and in some cases a process
converges while the other one does not, in other
cases what change are the chosen optimal sigma (and thus the
loglikelihoods and so on). This happens
when I use a regime with different levels for every node in the tree
but not when I use a single regime.
I'm pasting here an example of this phenomenon. It could be a stupid
error on my side but I'm not able to spot it right now.
Thanks,
Elena Grassi
"
> library("ouch")
> library("ape")
> t1 <- read.tree(text="(E:10,(D:5,C:5):5);")
> t2 <- read.tree(text="((C:5,D:5):5,E:10);")
> out1 <- ape2ouch(t1)
> out2 <- ape2ouch(t2)
> regime1 <- out@nodes
> regime1 <- out1@nodes
> regime2 <- rep("global", length(outree@nodes))
> regime2 <- rep("global", length(out2@nodes))
> expr1 <- c(NA, NA, 3, 5, 10)
> expr2 <- c(NA, NA, 10, 5, 3)
> models_data1 <- data.frame(expr1, regime1, regime2, row.names=out1@nodes)
> models_data2 <- data.frame(expr2, regime1, regime2, row.names=out2@nodes)
> m1 <- hansen(models_data1["expr1"], out1, models_data1["regime1"], 1, 1)
> m2 <- hansen(models_data2["expr2"], out2, models_data2["regime1"], 1, 1)
unsuccessful convergence, code 10, see documentation for `optim'
Warning message:
In hansen(models_data2["expr2"], out2, models_data2["regime1"], :
unsuccessful convergence
> m1
call:
hansen(data = models_data1["expr1"], tree = out1, regimes =
models_data1["regime1"],
sqrt.alpha = 1, sigma = 1)
nodes ancestors times labels regime1 expr1
1 1 <NA> 0.0 1 NA
2 2 1 0.5 2 NA
3 3 2 1.0 C 3 3
4 4 2 1.0 D 4 5
5 5 1 1.0 E 5 10
alpha:
[,1]
[1,] 3.734194
sigma squared:
[,1]
[1,] 9.233622e-30
theta:
$expr1
1 2 3 4 5
0.5049171 1.3917440 3.3191037 5.6847686 10.2324134
loglik deviance aic aic.c sic dof
100.5415 -201.0831 -187.0831 -209.4831 -193.3928 7.0000
> m22 <- hansen(models_data2["expr2"], out2, models_data2["regime2"], 1, 1)
> m11 <- hansen(models_data1["expr1"], out1, models_data1["regime2"], 1, 1)
> m11@loglik
[1] -7.47607
> m22@loglik
[1] -7.47607
> sessionInfo()
R version 2.15.1 (2012-06-22)
Platform: i486-pc-linux-gnu (32-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
LC_TIME=en_US.UTF-8
[4] LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C LC_NAME=C
LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=en_US.UTF-8
LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] ape_3.0-5 ouch_2.8-2 subplex_1.1-3
loaded via a namespace (and not attached):
[1] gee_4.13-18 grid_2.15.1 lattice_0.20-10 nlme_3.1-104
tools_2.15.1
"
--
www.biocut.unito.it
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
