3 ## Phylogenetic Generalized Least Squares
5 ## Copyright 2004 Julien Dutheil, and 2006-2012 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 corBrownian <- function(value = 1, phy, form = ~1)
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")
21 corMartins <- function(value, phy, form = ~1, fixed = FALSE)
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")
35 corGrafen <- function(value, phy, form = ~1, fixed = FALSE)
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")
50 Initialize.corPhyl <- function(object, data, ...)
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"))
59 attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
61 ## Obtaining the covariate(s)
62 attr(object, "covariate") <- getCovariate(object, data = data)
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.")
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]
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]
80 attr(object, "index") <- index
86 corMatrix.corBrownian <-
87 function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
97 index <- attr(object, "index")
101 corMatrix.corMartins <-
102 function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
114 index <- attr(object, "index")
118 corMatrix.corGrafen <-
119 function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
130 index <- attr(object, "index")
134 coef.corBrownian <- function(object, unconstrained = TRUE, ...)
136 if (!("corBrownian" %in% class(object)))
137 stop('object is not of class "corBrownian"')
141 coef.corMartins <- function(object, unconstrained = TRUE, ...)
143 if (!("corMartins" %in% class(object)))
144 stop('object is not of class "corMartins"')
146 if (attr(object, "fixed")) {
149 return(as.vector(object))
152 aux <- as.vector(object)
153 names(aux) <- "alpha"
157 coef.corGrafen <- function(object, unconstrained = TRUE, ...)
159 if (!("corGrafen" %in% class(object)))
160 stop('object is not of class "corGrafen"')
162 if (attr(object, "fixed")) {
165 return(as.vector(object))
168 aux <- exp(as.vector(object))
173 ### removed node.sons() and node.leafnumber() (2006-10-12)
175 ### changed by EP (2006-10-12):
177 compute.brlen <- function(phy, method = "Grafen", power = 1, ...)
179 if (!inherits(phy, "phylo"))
180 stop('object "phy" is not of class "phylo"')
181 Ntip <- length(phy$tip.label)
183 Nedge <- dim(phy$edge)[1]
184 if (is.numeric(method)) {
185 phy$edge.length <- rep(method, length.out = Nedge)
188 if (is.function(method)) {
189 phy$edge.length <- method(Nedge, ...)
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
200 (xx[phy$edge[, 1]]/m)^power - (xx[phy$edge[, 2]]/m)^power
207 corPagel <- function(value, phy, form = ~1, fixed = FALSE)
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")
220 corMatrix.corPagel <-
221 function(object, covariate = getCovariate(object), corr = TRUE, ...)
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]
234 coef.corPagel <- function(object, unconstrained = TRUE, ...)
237 if (attr(object, "fixed")) return(numeric(0))
238 else return(object[1])
241 names(aux) <- "lambda"
245 corBlomberg <- function(value, phy, form = ~1, fixed = FALSE)
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")
258 corMatrix.corBlomberg <-
259 function(object, covariate = getCovariate(object), corr = TRUE, ...)
261 index <- attr(object, "index")
263 stop("object has not been initialized")
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)
274 coef.corBlomberg <- function(object, unconstrained = TRUE, ...)
277 if (attr(object, "fixed")) return(numeric(0))
278 else return(object[1])