R/PGLS.R

   1 ## PGLS.R (2012-02-09)
   2
   3 ##   Phylogenetic Generalized Least Squares
   4
   5 ## Copyright 2004 Julien Dutheil, and 2006-2012 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 (!inherits(phy, "phylo"))
  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     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")
  25     if (value < 0) stop("the parameter alpha must be positive")
  26     if (!inherits(phy, "phylo"))
  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     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")
  39     if (value < 0) stop("parameter rho must be positive")
  40     value <- log(value) # Optimization under constraint, use exponential transform.
  41     if (!inherits(phy, "phylo"))
  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     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)) data <- parent.frame()
  67     ## Added by EP 29 May 2006:
  68     if (nrow(data) != length(phy$tip.label))
  69         stop("number of observations and number of tips in the tree are not equal.")
  70     ## END
  71     if (is.null(rownames(data))) {
  72         warning("No rownames supplied in data frame, data taken to be in the same order than in tree")
  73         attr(object, "index") <- 1:dim(data)[1]
  74     } else {
  75         index <- match(rownames(data), phy$tip.label)
  76         if (any(is.na(index))) {
  77             warning("Rownames in data frame do not match tree tip names; data taken to be in the same order as in tree")
  78             attr(object, "index") <- 1:dim(data)[1]
  79         } else {
  80             attr(object, "index") <- index
  81         }
  82     }
  83     object
  84 }
  85
  86 corMatrix.corBrownian <-
  87     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 has not been initialized.")
  93     tree <- attr(object, "tree")
  94     mat <- vcv.phylo(tree, corr = corr)
  95     n <- dim(mat)[1]
  96     ## reorder matrix:
  97     index <- attr(object, "index")
  98     mat[index, index]
  99 }
 100
 101 corMatrix.corMartins <-
 102     function(object, covariate = getCovariate(object), corr = TRUE, ...)
 103 {
 104     if (!("corMartins" %in% class(object)))
 105         stop('object is not of class "corMartins"')
 106     if (!any(attr(object, "index")))
 107         stop("object has not been initialized.")
 108     tree <- attr(object, "tree")
 109     dist <- cophenetic.phylo(tree)
 110     mat <- exp(-object[1] * dist)
 111     if (corr) mat <- cov2cor(mat)
 112     n <- dim(mat)[1]
 113     ## reorder matrix:
 114     index <- attr(object, "index")
 115     mat[index, index]
 116 }
 117
 118 corMatrix.corGrafen <-
 119     function(object, covariate = getCovariate(object), corr = TRUE, ...)
 120 {
 121     if (!("corGrafen" %in% class(object)))
 122         stop('object is not of class "corGrafen"')
 123     if (!any(attr(object, "index")))
 124         stop("object has not been initialized.")
 125     tree <- compute.brlen(attr(object, "tree"),
 126                           method = "Grafen", power = exp(object[1]))
 127     mat <- vcv.phylo(tree, corr = corr)
 128     n <- dim(mat)[1]
 129     ## reorder matrix:
 130     index <- attr(object, "index")
 131     mat[index, index]
 132 }
 133
 134 coef.corBrownian <- function(object, unconstrained = TRUE, ...)
 135 {
 136     if (!("corBrownian" %in% class(object)))
 137         stop('object is not of class "corBrownian"')
 138     numeric(0)
 139 }
 140
 141 coef.corMartins <- function(object, unconstrained = TRUE, ...)
 142 {
 143     if (!("corMartins" %in% class(object)))
 144         stop('object is not of class "corMartins"')
 145     if (unconstrained) {
 146         if (attr(object, "fixed")) {
 147             return(numeric(0))
 148         } else {
 149             return(as.vector(object))
 150         }
 151     }
 152     aux <- as.vector(object)
 153     names(aux) <- "alpha"
 154     aux
 155 }
 156
 157 coef.corGrafen <- function(object, unconstrained = TRUE, ...)
 158 {
 159     if (!("corGrafen" %in% class(object)))
 160         stop('object is not of class "corGrafen"')
 161     if (unconstrained) {
 162         if (attr(object, "fixed")) {
 163             return(numeric(0))
 164         } else {
 165             return(as.vector(object))
 166         }
 167     }
 168     aux <- exp(as.vector(object))
 169     names(aux) <- "rho"
 170     aux
 171 }
 172
 173 ### removed node.sons() and node.leafnumber()  (2006-10-12)
 174
 175 ### changed by EP (2006-10-12):
 176
 177 compute.brlen <- function(phy, method = "Grafen", power = 1, ...)
 178 {
 179     if (!inherits(phy, "phylo"))
 180         stop('object "phy" is not of class "phylo"')
 181     Ntip <- length(phy$tip.label)
 182     Nnode <- phy$Nnode
 183     Nedge <- dim(phy$edge)[1]
 184     if (is.numeric(method)) {
 185         phy$edge.length <- rep(method, length.out = Nedge)
 186         return(phy)
 187     }
 188     if (is.function(method)) {
 189         phy$edge.length <- method(Nedge, ...)
 190         return(phy)
 191     }
 192     if (is.character(method)) { # == "Grafen"
 193         tr <- reorder(phy, "pruningwise")
 194         xx <- .C("node_depth", as.integer(Ntip), as.integer(Nnode),
 195                  as.integer(tr$edge[, 1]), as.integer(tr$edge[, 2]),
 196                  as.integer(Nedge), double(Ntip + Nnode),
 197                  DUP = FALSE, PACKAGE = "ape")[[6]] - 1
 198         m <- Ntip - 1
 199         phy$edge.length <-
 200           (xx[phy$edge[, 1]]/m)^power - (xx[phy$edge[, 2]]/m)^power
 201         return(phy)
 202     }
 203 }
 204
 205 ## by EP:
 206
 207 corPagel <- function(value, phy, form = ~1, fixed = FALSE)
 208 {
 209     if (value < 0 || value > 1)
 210         stop("the value of lambda must be between 0 and 1.")
 211     if (!inherits(phy, "phylo"))
 212         stop('object "phy" is not of class "phylo"')
 213     attr(value, "formula") <- form
 214     attr(value, "fixed") <- fixed
 215     attr(value, "tree") <- phy
 216     class(value) <- c("corPagel", "corPhyl", "corStruct")
 217     value
 218 }
 219
 220 corMatrix.corPagel <-
 221     function(object, covariate = getCovariate(object), corr = TRUE, ...)
 222 {
 223     if (!any(attr(object, "index")))
 224         stop("object has not been initialized")
 225     mat <- vcv.phylo(attr(object, "tree"), corr = corr)
 226     index <- attr(object, "index")
 227     mat <- mat[index, index]
 228     tmp <- diag(mat)
 229     mat <- object[1]*mat
 230     diag(mat) <- tmp
 231     mat
 232 }
 233
 234 coef.corPagel <- function(object, unconstrained = TRUE, ...)
 235 {
 236     if (unconstrained) {
 237         if (attr(object, "fixed")) return(numeric(0))
 238         else return(object[1])
 239     }
 240     aux <- object[1]
 241     names(aux) <- "lambda"
 242     aux
 243 }
 244
 245 corBlomberg <- function(value, phy, form = ~1, fixed = FALSE)
 246 {
 247     if (value <= 0)
 248         stop("the value of g must be greater than 0.")
 249     if (!inherits(phy, "phylo"))
 250         stop('object "phy" is not of class "phylo"')
 251     attr(value, "formula") <- form
 252     attr(value, "fixed") <- fixed
 253     attr(value, "tree") <- phy
 254     class(value) <- c("corBlomberg", "corPhyl", "corStruct")
 255     value
 256 }
 257
 258 corMatrix.corBlomberg <-
 259     function(object, covariate = getCovariate(object), corr = TRUE, ...)
 260 {
 261     index <- attr(object, "index")
 262     if (is.null(index))
 263         stop("object has not been initialized")
 264     if (object[1] <= 0)
 265         stop("the optimization has reached a value <= 0 for parameter 'g':
 266 probably need to set 'fixed = TRUE' in corBlomberg().")
 267     phy <- attr(object, "tree")
 268     d <- (dist.nodes(phy)[length(phy$tip.label) + 1, ])^(1/object[1])
 269     phy$edge.length <- d[phy$edge[, 2]] - d[phy$edge[, 1]]
 270     mat <- vcv.phylo(phy, corr = corr)
 271     mat[index, index]
 272 }
 273
 274 coef.corBlomberg <- function(object, unconstrained = TRUE, ...)
 275 {
 276     if (unconstrained) {
 277         if (attr(object, "fixed")) return(numeric(0))
 278         else return(object[1])
 279     }
 280     aux <- object[1]
 281     names(aux) <- "g"
 282     aux
 283 }