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 ## fix by Michael Phelan
16 diag(W) <- 0 # ensure diagonal is zero
19 if(dim(y)[2] != 1) stop("Error: y must be a single column vector.")
20 D <- solve(diag(apply(t(W),2,sum)))
21 Wnorm <- D %*% W #Row normalize W matrix
24 y <- y-matrix(rep(1, n)) %*% apply(y,2,mean) # Deviations from mean
27 Wlam <- eigen(Wnorm)$values # eigenvalues of W
29 # Find distinct eigenvalues
33 if(abs(Im(sorted[ii])) > 1e-12) {
34 warning(paste("Complex eigenvalue coerced to real:", Im(sorted[ii])))
36 sorted[ii] <- Re(sorted[ii]) # Remove imaginary part
38 sorted <- as.real(sorted)
40 Distinct <- numeric(0)
42 Distinct[2] <- sorted[1]
45 if(sorted[ii] - Distinct[nDistinct] > tolerance) {
46 nDistinct <- nDistinct + 1
47 Distinct[nDistinct] <- sorted[ii]
51 # Search for minimum of LL
53 likelihood <- function(rhohat) {
56 prod <- 1 - rhohat * Wlam[j]
57 DetProd <- DetProd * prod
59 absValDet <- abs(DetProd) #[abs to allow rho > 1]
60 logDet <- log(absValDet)
61 LL <- log(t(y) %*% y - 2 * rhohat * t(y) %*% Wy + rhohat * rhohat * t(Wy) %*% Wy) - logDet*2/n
65 GoldenSearch <- function(ax, cx) {
66 # Golden section search over the interval ax to cx
67 # Return rhohat and likelihood value.
72 if(abs(cx - bx) > abs(bx - ax)) {
74 x2 <- bx + (1-r)*(cx - bx)
77 x1 <- bx - (1-r)*(bx - ax)
81 while(abs(x3 - x0) > gold.tol*(abs(x1) + abs(x2))) {
85 x2 <- r * x1 + (1 - r) * x3
91 x1 <- r * x2 + (1 - r) * x0
103 return(list(rho=xmin, LL=likelihood))
107 for(ii in 2:(nDistinct -1)) {# Search between pairs of roots
108 # [ constrain do not use positive roots < 1]
110 cx <- 1/Distinct[ii+1]
111 GS <- GoldenSearch(ax, cx)
119 return(list(rhohat=rho, Wnorm=Wnorm, residuals=res))
135 #y<-c(-0.12,0.36,-0.1,0.04,-0.15,0.29,-0.11,-0.06)
136 #compar.cheverud(y,W)
139 #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)
140 #compar.cheverud(y,W)