# R/Cheverud.R In gjuggler/ape: Analyses of Phylogenetics and Evolution

## Cheverud.R (2004-10-29)

##    Cheverud's 1985 Autoregression Model

## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.

# This function is adapted from a MatLab code from
# Rholf, F. J. (2001) Comparative Methods for the Analysis of Continuous Variables: Geometric Interpretations.
# Evolution 55(11): 2143-2160
compar.cheverud <- function(y, W, tolerance=1e-6, gold.tol=1e-4)
{
## fix by Michael Phelan
diag(W) <- 0 # ensure diagonal is zero
## end of fix
y <- as.matrix(y)
if(dim(y)[2] != 1) stop("Error: y must be a single column vector.")
D <- solve(diag(apply(t(W),2,sum)))
Wnorm <- D %*% W #Row normalize W matrix
n <- dim(y)[1]
m <- dim(y)[2]
y <- y-matrix(rep(1, n)) %*% apply(y,2,mean) # Deviations from mean
Wy <- Wnorm %*% y

Wlam <- eigen(Wnorm)\$values # eigenvalues of W

# Find distinct eigenvalues
sorted <- sort(Wlam)
# Check real:
for (ii in 1:n) {
if(abs(Im(sorted[ii])) > 1e-12) {
warning(paste("Complex eigenvalue coerced to real:", Im(sorted[ii])))
}
sorted[ii] <- Re(sorted[ii]) # Remove imaginary part
}
sorted <- as.real(sorted)

Distinct <- numeric(0)
Distinct[1] <- -Inf
Distinct[2] <- sorted[1]
nDistinct <- 2
for(ii in 2:n) {
if(sorted[ii] - Distinct[nDistinct] > tolerance) {
nDistinct <- nDistinct + 1
Distinct[nDistinct] <- sorted[ii]
}
}

# Search for minimum of LL

likelihood <- function(rhohat) {
DetProd <- 1
for(j in 1:n) {
prod <- 1 - rhohat * Wlam[j]
DetProd <- DetProd * prod
}
absValDet <- abs(DetProd) #[abs to allow rho > 1]
logDet <- log(absValDet)
LL <- log(t(y) %*% y - 2 * rhohat * t(y) %*% Wy + rhohat * rhohat * t(Wy) %*% Wy) - logDet*2/n
return(LL)
}

GoldenSearch <- function(ax, cx) {
# Golden section search over the interval ax to cx
# Return rhohat and likelihood value.
r <- 0.61803399
x0 <- ax
x3 <- cx
bx <- (ax + cx)/2
if(abs(cx - bx) > abs(bx - ax)) {
x1 <- bx
x2 <- bx + (1-r)*(cx - bx)
} else {
x2 <- bx
x1 <- bx - (1-r)*(bx - ax)
}
f1 <- likelihood(x1)
f2 <- likelihood(x2)
while(abs(x3 - x0) > gold.tol*(abs(x1) + abs(x2))) {
if(f2 < f1) {
x0 <- x1
x1 <- x2
x2 <- r * x1 + (1 - r) * x3
f1 <- f2
f2 <- likelihood(x2)
} else {
x3 <- x2
x2 <- x1
x1 <- r * x2 + (1 - r) * x0
f2 <- f1
f1 <- likelihood(x1)
}
}
if(f1 < f2) {
likelihood <- f1
xmin <- x1
} else {
likelihood <- f2
xmin <- x2
}
return(list(rho=xmin, LL=likelihood))
}

LL <- Inf
for(ii in 2:(nDistinct -1)) {# Search between pairs of roots
# [ constrain do not use positive roots < 1]
ax <- 1/Distinct[ii]
cx <- 1/Distinct[ii+1]
GS <- GoldenSearch(ax, cx)
if(GS\$LL < LL) {
LL <- GS\$LL
rho <- GS\$rho
}
}
# Compute residuals:
res <- y - rho * Wy
return(list(rhohat=rho, Wnorm=Wnorm, residuals=res))
}

#For debugging:
#W<- matrix(c(
#  0,1,1,2,0,0,0,0,
#  1,0,1,2,0,0,0,0,
#  1,1,0,2,0,0,0,0,
#  2,2,2,0,0,0,0,0,
#  0,0,0,0,0,1,1,2,
#  0,0,0,0,1,0,1,2,
#  0,0,0,0,1,1,0,2,
#  0,0,0,0,2,2,2,0
#),8)
#W <- 1/W
#W[W == Inf] <- 0
#y<-c(-0.12,0.36,-0.1,0.04,-0.15,0.29,-0.11,-0.06)
#compar.cheverud(y,W)
#
#y<-c(10,8,3,4)
#W <- matrix(c(1,1/6,1/6,1/6,1/6,1,1/2,1/2,1/6,1/2,1,1,1/6,1/2,1,1), 4)
#compar.cheverud(y,W)
gjuggler/ape documentation built on May 17, 2019, 6:03 a.m.