1 ## Cheverud.R (2004-10-29)
3 ## Cheverud's 1985 Autoregression Model
5 ## Copyright 2004 Julien Dutheil
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 # This function is adapted from a MatLab code from
11 # Rholf, F. J. (2001) Comparative Methods for the Analysis of Continuous Variables: Geometric Interpretations.
12 # Evolution 55(11): 2143-2160
13 compar.cheverud <- function(y, W, tolerance=1e-6, gold.tol=1e-4)
15 W <- W - diag(W) # ensure diagonal is zero
17 if(dim(y)[2] != 1) stop("Error: y must be a single column vector.")
18 D <- solve(diag(apply(t(W),2,sum)))
19 Wnorm <- D %*% W #Row normalize W matrix
22 y <- y-matrix(rep(1, n)) %*% apply(y,2,mean) # Deviations from mean
25 Wlam <- eigen(Wnorm)$values # eigenvalues of W
27 # Find distinct eigenvalues
31 if(abs(Im(sorted[ii])) > 1e-12) {
32 warning(paste("Complex eigenvalue coerced to real:", Im(sorted[ii])))
34 sorted[ii] <- Re(sorted[ii]) # Remove imaginary part
36 sorted <- as.real(sorted)
38 Distinct <- numeric(0)
40 Distinct[2] <- sorted[1]
43 if(sorted[ii] - Distinct[nDistinct] > tolerance) {
44 nDistinct <- nDistinct + 1
45 Distinct[nDistinct] <- sorted[ii]
49 # Search for minimum of LL
51 likelihood <- function(rhohat) {
54 prod <- 1 - rhohat * Wlam[j]
55 DetProd <- DetProd * prod
57 absValDet <- abs(DetProd) #[abs to allow rho > 1]
58 logDet <- log(absValDet)
59 LL <- log(t(y) %*% y - 2 * rhohat * t(y) %*% Wy + rhohat * rhohat * t(Wy) %*% Wy) - logDet*2/n
63 GoldenSearch <- function(ax, cx) {
64 # Golden section search over the interval ax to cx
65 # Return rhohat and likelihood value.
70 if(abs(cx - bx) > abs(bx - ax)) {
72 x2 <- bx + (1-r)*(cx - bx)
75 x1 <- bx - (1-r)*(bx - ax)
79 while(abs(x3 - x0) > gold.tol*(abs(x1) + abs(x2))) {
83 x2 <- r * x1 + (1 - r) * x3
89 x1 <- r * x2 + (1 - r) * x0
101 return(list(rho=xmin, LL=likelihood))
105 for(ii in 2:(nDistinct -1)) {# Search between pairs of roots
106 # [ constrain do not use positive roots < 1]
108 cx <- 1/Distinct[ii+1]
109 GS <- GoldenSearch(ax, cx)
117 return(list(rhohat=rho, Wnorm=Wnorm, residuals=res))
133 #y<-c(-0.12,0.36,-0.1,0.04,-0.15,0.29,-0.11,-0.06)
134 #compar.cheverud(y,W)
137 #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)
138 #compar.cheverud(y,W)