Discussion:
[R-sig-phylo] Question regarding log likelihood in OUCH package
Danielle Miller
2018-08-12 14:32:10 UTC
Permalink
Hi,

I’m interested in using the OUCH package to estimate BM and OU parameters for a specific trait among many different trees.
My goal is to determined which model is the most suitable for each tree, applying likelihood ratio test.

As I’m a new user in R when it comes to phylogenetic analysis, I started by running the documentation example (Hansen, documentation page 11) and was surprised to see that the loglikelihood was a positive number

BM:
$call
brown(data = otd[c("tarsusL", "beakD")], tree = ot)

$sigma.squared
[,1] [,2]
[1,] 0.02878091 0.08897504
[2,] 0.08897504 0.43711838

$theta
$theta$tarsusL
[1] 3.020419

$theta$beakD
[1] 1.826695


$loglik
[1] 9.90115

As this number is crucial for further analysis - Is this a transformation of the resulting log likelihood? (e.g. -2 * log(L) as described in the paper) or am I missing something here..?

In addition I have another issue, I have a tree constructed of ~400 viral genomes and their corresponding trait values. When I’m running the documentation script with my own data (in the same format)
I get the following error:
Error in solve.default(v, e) :
system is computationally singular: reciprocal condition number = 1.59061e-17

I guess it says that my variance covariance matrix is not inversable, hence I manually tried to adjust the retol parameter in the Hansen function in order to make it work (however I’ll need to second guess my results?), but I still get the same error.
Code example:
h1 <- hansen(
+ tree=ot,
+ data=otd[c("k5")],
+ regimes=otd["regimes"],
+ fit=TRUE,
+ sqrt.alpha=1,
+ sigma=1,
+ maxit=500000,
+ reltol=1e-20,
+ method="Nelder-Mead"
+ )

I’ll be thankful for any advice or answer,
Danielle





_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/
Liam J. Revell
2018-08-12 16:10:12 UTC
Permalink
Hi Danielle.

With regard to the first problem, that is not an error. Likelihoods
obtained from continuous probability densities are not probabilities
(these can only be obtained by integrating the density function on a
finite interval) and thus can take values >1 (and thus log-likelihoods >0).

In my experience, the most common reason for the second problem is a
tree with zero-length terminal edges (and thus covariances between
species that are equal to their respective variances). The solution to
this *should not* be to add a very small value to the offending terminal
edges, as this can give very high weight to the associated tips. It
might be by pruning one tip or the other from the analysis, or by using
some reasonable criterion to modify the terminal edge lengths.

All the best,

Liam J. Revell
Associate Professor, University of Massachusetts Boston
Profesor Asistente, Universidad Católica de la Ssma Concepción
web: http://faculty.umb.edu/liam.revell/, http://www.phytools.org
Post by Danielle Miller
Hi,
I’m interested in using the OUCH package to estimate BM and OU parameters for a specific trait among many different trees.
My goal is to determined which model is the most suitable for each tree, applying likelihood ratio test.
As I’m a new user in R when it comes to phylogenetic analysis, I started by running the documentation example (Hansen, documentation page 11) and was surprised to see that the loglikelihood was a positive number
$call
brown(data = otd[c("tarsusL", "beakD")], tree = ot)
$sigma.squared
[,1] [,2]
[1,] 0.02878091 0.08897504
[2,] 0.08897504 0.43711838
$theta
$theta$tarsusL
[1] 3.020419
$theta$beakD
[1] 1.826695
$loglik
[1] 9.90115
As this number is crucial for further analysis - Is this a transformation of the resulting log likelihood? (e.g. -2 * log(L) as described in the paper) or am I missing something here..?
In addition I have another issue, I have a tree constructed of ~400 viral genomes and their corresponding trait values. When I’m running the documentation script with my own data (in the same format)
system is computationally singular: reciprocal condition number = 1.59061e-17
I guess it says that my variance covariance matrix is not inversable, hence I manually tried to adjust the retol parameter in the Hansen function in order to make it work (however I’ll need to second guess my results?), but I still get the same error.
h1 <- hansen(
+ tree=ot,
+ data=otd[c("k5")],
+ regimes=otd["regimes"],
+ fit=TRUE,
+ sqrt.alpha=1,
+ sigma=1,
+ maxit=500000,
+ reltol=1e-20,
+ method="Nelder-Mead"
+ )
I’ll be thankful for any advice or answer,
Danielle
_______________________________________________
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
_______________________________________________
R-sig-phylo mailing list - R-sig-***@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-***@r-project.org/
Loading...