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