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