R/PGLS.R

   1 ## PGLS.R (2008-05-08)
   2
   3 ##   Phylogenetic Generalized Least Squares
   4
   5 ## Copyright 2004 Julien Dutheil, and 2006-2008 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 corBrownian <- function(value = 1, phy, form=~1)
  11 {
  12   if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
  13   attr(value, "formula") <- form
  14   attr(value, "fixed")   <- TRUE
  15   attr(value, "tree")    <- phy
  16   class(value) <- c("corBrownian", "corPhyl", "corStruct")
  17   return(value)
  18 }
  19
  20 corMartins <- function(value, phy, form=~1, fixed=FALSE)
  21 {
  22   if(length(value) > 1) stop("ERROR!!! Only one parameter is allowed in corPGLS structure.")
  23   if(value < 0) stop("ERROR!!! Parameter alpha must be positive.")
  24   if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
  25   attr(value, "formula") <- form
  26   attr(value, "fixed")   <- fixed
  27   attr(value, "tree")    <- phy
  28   class(value) <- c("corMartins", "corPhyl", "corStruct")
  29   return(value)
  30 }
  31
  32 corGrafen <- function(value, phy, form=~1, fixed=FALSE)
  33 {
  34   if(length(value) > 1) stop("ERROR!!! Only one parameter is allowed in corGrafen structure.")
  35   if(value < 0) stop("ERROR!!! Parameter rho must be positive.")
  36   value <- log(value) # Optimization under constraint, use exponential transform.
  37   if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
  38   attr(value, "formula") <- form
  39   attr(value, "fixed")   <- fixed
  40   attr(value, "tree")    <- phy
  41   class(value) <- c("corGrafen", "corPhyl", "corStruct")
  42   return(value)
  43 }
  44
  45 Initialize.corPhyl <- function(object, data, ...)
  46 {
  47   # The same as in Initialize corStruct:
  48   form <- formula(object)
  49   ## Obtaining the group information, if any
  50   if(!is.null(getGroupsFormula(form))) {
  51     attr(object, "groups") <- getGroups(object, form, data = data)
  52     attr(object, "Dim")    <- Dim(object, attr(object, "groups"))
  53   } else { # no groups
  54     attr(object, "Dim")    <- Dim(object, as.factor(rep(1, nrow(data))))
  55   }
  56   ## Obtaining the covariate(s)
  57   attr(object, "covariate") <- getCovariate(object, data = data)
  58
  59   # Specific to corPhyl:
  60   phy <- attr(object, "tree")
  61   if (is.null(data))
  62     data <- parent.frame()
  63   ## Added by EP 29 May 2006:
  64   if (nrow(data) != length(phy$tip.label))
  65     stop("number of observations and number of tips in the tree are not equal.")
  66   ## END
  67   if(is.null(rownames(data))) {
  68     warning("No row names supplied in dataframe, data taken to be in the same order as in tree.")
  69     attr(object, "index") <- 1:dim(data)[1]
  70   } else {
  71     index <- match(rownames(data), phy$tip.label)
  72     if(any(is.na(index))) {
  73       warning("Row names in dataframe do not match tree tip names. data taken to be in the same order as in tree.")
  74       attr(object, "index") <- 1:dim(data)[1]
  75     } else {
  76       attr(object, "index") <- index
  77     }
  78   }
  79   return(object)
  80 }
  81
  82 corMatrix.corBrownian <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
  83 {
  84   if (!("corBrownian" %in% class(object))) stop("ERROR!!! Object is not of class \"corBrownian\".")
  85   if(!any(attr(object, "index"))) stop("ERROR!!! object have not been initialized.")
  86   tree <- attr(object, "tree")
  87   mat <- vcv.phylo(tree, cor = corr)
  88   n <- dim(mat)[1]
  89   # reorder matrix:
  90   matr <- matrix(nrow=n, ncol=n)
  91   index <- attr(object, "index")
  92   for(i in 1:n)
  93     for(j in i:n)
  94       matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
  95   return(matr)
  96 }
  97
  98 corMatrix.corMartins <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
  99 {
 100   if (!("corMartins" %in% class(object))) stop("ERROR!!! Object is not of class \"corMartins\".")
 101   if(!any(attr(object, "index"))) stop("ERROR!!! object have not been initialized.")
 102   tree <- attr(object, "tree")
 103   dist <- cophenetic.phylo(tree)
 104   mat <- exp(-object[1] * dist)
 105   if(corr) mat <- cov2cor(mat)
 106   n <- dim(mat)[1]
 107   # reorder matrix:
 108   matr <- matrix(nrow=n, ncol=n)
 109   index <- attr(object, "index")
 110   for(i in 1:n)
 111     for(j in i:n)
 112       matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
 113   return(matr)
 114 }
 115
 116 corMatrix.corGrafen <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
 117 {
 118   if (!("corGrafen" %in% class(object))) stop("ERROR!!! Object is not of class \"corGrafen\".")
 119   if(!any(attr(object, "index"))) stop("ERROR!!! object have not been initialized.")
 120   tree <- compute.brlen(attr(object, "tree"), method = "Grafen", power = exp(object[1]))
 121   mat <- vcv.phylo(tree, cor = corr)
 122   n <- dim(mat)[1]
 123   # reorder matrix:
 124   matr <- matrix(nrow=n, ncol=n)
 125   index <- attr(object, "index")
 126   for(i in 1:n)
 127     for(j in i:n)
 128       matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
 129   return(matr)
 130 }
 131
 132 coef.corBrownian <- function(object, unconstrained = TRUE, ...)
 133 {
 134   if (!("corBrownian" %in% class(object))) stop("ERROR!!! Object is not of class \"corBrownian\".")
 135   return(numeric(0))
 136 }
 137
 138 coef.corMartins <- function(object, unconstrained = TRUE, ...)
 139 {
 140   if (!("corMartins" %in% class(object))) stop("ERROR!!! Object is not of class \"corMartins\".")
 141   if(unconstrained) {
 142     if(attr(object, "fixed")) {
 143       return(numeric(0))
 144     } else {
 145       return(as.vector(object))
 146     }
 147   }
 148   aux <- as.vector(object)
 149   names(aux) <- "alpha"
 150   return(aux)
 151 }
 152
 153 coef.corGrafen <- function(object, unconstrained = TRUE, ...)
 154 {
 155   if (!("corGrafen" %in% class(object))) stop("ERROR!!! Object is not of class \"corGrafen\".")
 156   if(unconstrained) {
 157     if(attr(object, "fixed")) {
 158       return(numeric(0))
 159     } else {
 160       return(as.vector(object))
 161     }
 162   }
 163   aux <- exp(as.vector(object))
 164   names(aux) <- "rho"
 165   return(aux)
 166 }
 167
 168 ### removed node.sons() and node.leafnumber()  (2006-10-12)
 169
 170 ### changed by EP (2006-10-12):
 171
 172 compute.brlen <- function(phy, method = "Grafen", power = 1, ...)
 173 {
 174     if (!"phylo" %in% class(phy))
 175       stop('object "phy" is not of class "phylo"')
 176     Ntip <- length(phy$tip.label)
 177     Nnode <- phy$Nnode
 178     Nedge <- dim(phy$edge)[1]
 179     if (is.numeric(method)) {
 180         phy$edge.length <- rep(method, length.out = Nedge)
 181         return(phy)
 182     }
 183     if (is.function(method)) {
 184         phy$edge.length <- method(Nedge, ...)
 185         return(phy)
 186     }
 187     if (is.character(method)) { # == "Grafen"
 188         tr <- reorder(phy, "pruningwise")
 189         xx <- .C("node_depth", as.integer(Ntip), as.integer(Nnode),
 190                  as.integer(tr$edge[, 1]), as.integer(tr$edge[, 2]),
 191                  as.integer(Nedge), double(Ntip + Nnode),
 192                  DUP = FALSE, PACKAGE = "ape")[[6]] - 1
 193         m <- Ntip - 1
 194         phy$edge.length <-
 195           (xx[phy$edge[, 1]]/m)^power - (xx[phy$edge[, 2]]/m)^power
 196         return(phy)
 197     }
 198 }
 199
 200 ## by EP:
 201
 202 corPagel <- function(value, phy, form = ~1, fixed = FALSE)
 203 {
 204     if (value < 0 || value > 1)
 205         stop("the value of lambda must be between 0 and 1.")
 206     if (!("phylo" %in% class(phy)))
 207         stop('object "phy" is not of class "phylo"')
 208     attr(value, "formula") <- form
 209     attr(value, "fixed") <- fixed
 210     attr(value, "tree") <- phy
 211     class(value) <- c("corPagel", "corPhyl", "corStruct")
 212     value
 213 }
 214
 215 corMatrix.corPagel <-
 216     function(object, covariate = getCovariate(object), corr = TRUE, ...)
 217 {
 218     if (!any(attr(object, "index")))
 219         stop("object has not been initialized")
 220     mat <- vcv.phylo(attr(object, "tree"), cor = corr)
 221     index <- attr(object, "index")
 222     mat <- mat[index, index]
 223     tmp <- diag(mat)
 224     mat <- object[1]*mat
 225     diag(mat) <- tmp
 226     mat
 227 }
 228
 229 coef.corPagel <- function(object, unconstrained = TRUE, ...)
 230 {
 231     if (unconstrained) {
 232         if (attr(object, "fixed")) return(numeric(0))
 233         else return(object[1])
 234     }
 235     aux <- object[1]
 236     names(aux) <- "lambda"
 237     aux
 238 }
 239
 240 corBlomberg <- function(value, phy, form = ~1, fixed = FALSE)
 241 {
 242     if (value <= 0)
 243         stop("the value of g must be greater than 0.")
 244     if (!("phylo" %in% class(phy)))
 245         stop('object "phy" is not of class "phylo"')
 246     attr(value, "formula") <- form
 247     attr(value, "fixed") <- fixed
 248     attr(value, "tree") <- phy
 249     class(value) <- c("corBlomberg", "corPhyl", "corStruct")
 250     value
 251 }
 252
 253 corMatrix.corBlomberg <-
 254     function(object, covariate = getCovariate(object), corr = TRUE, ...)
 255 {
 256     index <- attr(object, "index")
 257     if (is.null(index))
 258         stop("object has not been initialized")
 259     if (object[1] <= 0)
 260         stop("the optimization has reached a value <= 0 for parameter 'g':
 261 probably need to set 'fixed = TRUE' in corBlomberg().")
 262     phy <- attr(object, "tree")
 263     d <- (dist.nodes(phy)[length(phy$tip.label) + 1, ])^(1/object[1])
 264     phy$edge.length <- d[phy$edge[, 2]] - d[phy$edge[, 1]]
 265     mat <- vcv.phylo(phy, cor = corr)
 266     mat[index, index]
 267 }
 268
 269 coef.corBlomberg <- function(object, unconstrained = TRUE, ...)
 270 {
 271     if (unconstrained) {
 272         if (attr(object, "fixed")) return(numeric(0))
 273         else return(object[1])
 274     }
 275     aux <- object[1]
 276     names(aux) <- "g"
 277     aux
 278 }