/auteur/R/likelihood.R

 1#primary function for computing the likelihood of data, given a root state, VCV matrix, and Brownian motion model 
 2#author: LJ HARMON 2009 and JM EASTMAN 2010
 3
 4bm.lik.fn <-
 5function(root, dat, vcv, SE) { 
 6# mod 12.02.2010 JM Eastman: using determinant()$modulus rather than det() to stay in log-space
 7	y=dat
 8	
 9	b <- vcv
10	if(any(SE!=0)) diag(b)=diag(b)+SE^2
11	w <- rep(root, nrow(b))
12	num <- -t(y-w)%*%solve(b)%*%(y-w)/2
13	den <- 0.5*(length(y)*log(2*pi) + as.numeric(determinant(b)$modulus))
14	list(
15		 root=root,
16		 lnL = (num-den)[1,1]
17		 )	
18}
19
20#primary function for computing the likelihood of data, using REML, under Brownian motion model 
21#author: LJ HARMON, LJ REVELL, and JM EASTMAN 2011
22
23bm.reml.fn <- 
24function(pruningwise.tre, pruningwise.rates, ic) {
25	new.tre=updatetree(pruningwise.tre, pruningwise.rates)
26	new.var=pic_variance(new.tre)
27	reml=dnorm(ic, mean=0, sd=sqrt(new.var), log=TRUE)
28	list(
29		 root=NA,
30		 lnL = sum(reml)
31		 )	
32}
33
34
35#compute expected PIC variance given tree: used for bm.reml.fn()  
36#author: JM EASTMAN 2011
37
38pic_variance <- function(phy)
39{
40	nb.tip <- Ntip(phy)
41    nb.node <- phy$Nnode
42	if (nb.node != nb.tip - 1) stop("'phy' is not rooted and fully dichotomous")
43	
44    phy <- reorder(phy, "pruningwise")
45	
46	ans <- .C("pic_variance", as.integer(nb.tip), as.integer(nb.node),
47              as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
48              as.double(phy$edge.length), double(nb.node),
49              PACKAGE = "auteur")
50	
51    var <- ans[[6]]
52    names(var) <- as.character(1:nb.node + nb.tip)
53    var
54}
55
56
57ou.lik.fn <- function(dat, phy, epochs, phyobject, ancestors, vcvSE, rate, alpha, regimes) {
58# mod 02.22.2010 JM Eastman
59	phyobject$regimes=regimes
60	ww <- ouweights(phy, epochs, phyobject, ancestors, alpha)
61	w <- ww$W
62	v <- oumat(vcvSE, alpha, rate)
63	y <- matrix(ww$optima, ncol=1)
64	e <- w%*%y - dat
65	dim(e) <- n <- length(dat)
66	q <- e%*%solve(v,e)
67	det.v <- determinant(v,logarithm=TRUE)
68	dev <- n*log(2*pi)+as.numeric(det.v$modulus)+q[1,1]
69	list(
70		 weight=w,
71		 vcov=v,
72		 resid=e,
73		 predict=y,
74		 lnL=-0.5*dev
75		 )
76}
77
78
79#general phylogenetic utility for determining whether to accept ('r'=1) or reject ('r'=0) a proposed model in Markov chain Monte Carlo sampling
80#author: JM EASTMAN 2010
81
82assess.lnR <-
83function(lnR) {
84	if(is.na(lnR) || abs(lnR)==Inf) {
85		error=TRUE
86		r=0
87	} else {
88		if(lnR < -20) {
89			r=0 
90		} else if(lnR > 0) {
91			r=1 
92		} else {
93			r=exp(lnR)
94		}
95		error=FALSE
96	}
97	return(list(r=r, error=error))
98}
99