R/Cheverud.R

   1 ## Cheverud.R (2004-10-29)
   2
   3 ##    Cheverud's 1985 Autoregression Model
   4
   5 ## Copyright 2004 Julien Dutheil
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  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)
  14 {
  15   W <- W - diag(W) # ensure diagonal is zero
  16   y <- as.matrix(y)
  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
  20   n <- dim(y)[1]
  21   m <- dim(y)[2]
  22   y <- y-matrix(rep(1, n)) %*% apply(y,2,mean) # Deviations from mean
  23   Wy <- Wnorm %*% y
  24
  25   Wlam <- eigen(Wnorm)$values # eigenvalues of W
  26
  27   # Find distinct eigenvalues
  28   sorted <- sort(Wlam)
  29   # Check real:
  30   for (ii in 1:n) {
  31     if(abs(Im(sorted[ii])) > 1e-12) {
  32       warning(paste("Complex eigenvalue coerced to real:", Im(sorted[ii])))
  33           }
  34     sorted[ii] <- Re(sorted[ii]) # Remove imaginary part
  35   }
  36   sorted <- as.real(sorted)
  37
  38         Distinct <- numeric(0)
  39   Distinct[1] <- -Inf
  40   Distinct[2] <- sorted[1]
  41   nDistinct <- 2
  42   for(ii in 2:n) {
  43     if(sorted[ii] - Distinct[nDistinct] > tolerance) {
  44       nDistinct <- nDistinct + 1
  45       Distinct[nDistinct] <- sorted[ii]
  46     }
  47   }
  48
  49   # Search for minimum of LL
  50
  51   likelihood <- function(rhohat) {
  52     DetProd <- 1
  53     for(j in 1:n) {
  54       prod <- 1 - rhohat * Wlam[j]
  55       DetProd <- DetProd * prod
  56     }
  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
  60     return(LL)
  61   }
  62
  63   GoldenSearch <- function(ax, cx) {
  64     # Golden section search over the interval ax to cx
  65     # Return rhohat and likelihood value.
  66     r <- 0.61803399
  67     x0 <- ax
  68     x3 <- cx
  69     bx <- (ax + cx)/2
  70     if(abs(cx - bx) > abs(bx - ax)) {
  71       x1 <- bx
  72       x2 <- bx + (1-r)*(cx - bx)
  73     } else {
  74       x2 <- bx
  75       x1 <- bx - (1-r)*(bx - ax)
  76     }
  77     f1 <- likelihood(x1)
  78     f2 <- likelihood(x2)
  79     while(abs(x3 - x0) > gold.tol*(abs(x1) + abs(x2))) {
  80       if(f2 < f1) {
  81         x0 <- x1
  82         x1 <- x2
  83         x2 <- r * x1 + (1 - r) * x3
  84         f1 <- f2
  85         f2 <- likelihood(x2)
  86       } else {
  87         x3 <- x2
  88         x2 <- x1
  89         x1 <- r * x2 + (1 - r) * x0
  90         f2 <- f1
  91         f1 <- likelihood(x1)
  92       }
  93     }
  94     if(f1 < f2) {
  95       likelihood <- f1
  96       xmin <- x1
  97     } else {
  98       likelihood <- f2
  99       xmin <- x2
 100     }
 101     return(list(rho=xmin, LL=likelihood))
 102   }
 103
 104   LL <- Inf
 105   for(ii in 2:(nDistinct -1)) {# Search between pairs of roots
 106     # [ constrain do not use positive roots < 1]
 107     ax <- 1/Distinct[ii]
 108     cx <- 1/Distinct[ii+1]
 109     GS <- GoldenSearch(ax, cx)
 110     if(GS$LL < LL) {
 111       LL <- GS$LL
 112       rho <- GS$rho
 113     }
 114   }
 115   # Compute residuals:
 116   res <- y - rho * Wy
 117   return(list(rhohat=rho, Wnorm=Wnorm, residuals=res))
 118 }
 119
 120 #For debugging:
 121 #W<- matrix(c(
 122 #  0,1,1,2,0,0,0,0,
 123 #  1,0,1,2,0,0,0,0,
 124 #  1,1,0,2,0,0,0,0,
 125 #  2,2,2,0,0,0,0,0,
 126 #  0,0,0,0,0,1,1,2,
 127 #  0,0,0,0,1,0,1,2,
 128 #  0,0,0,0,1,1,0,2,
 129 #  0,0,0,0,2,2,2,0
 130 #),8)
 131 #W <- 1/W
 132 #W[W == Inf] <- 0
 133 #y<-c(-0.12,0.36,-0.1,0.04,-0.15,0.29,-0.11,-0.06)
 134 #compar.cheverud(y,W)
 135 #
 136 #y<-c(10,8,3,4)
 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)
 139