3 ## Phylogenetic Generalized Least Squares
5 ## Copyright 2004 Julien Dutheil, and 2006-2008 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 (!("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")
21 corMartins <- function(value, phy, form=~1, fixed=FALSE)
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")
35 corGrafen <- function(value, phy, form=~1, fixed=FALSE)
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")
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")
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.")
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]
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]
81 attr(object, "index") <- index
87 corMatrix.corBrownian <- 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 have not been initialized.")
93 tree <- attr(object, "tree")
94 mat <- vcv.phylo(tree, cor = corr)
97 matr <- matrix(nrow=n, ncol=n)
98 index <- attr(object, "index")
101 matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
105 corMatrix.corMartins <-
106 function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
118 matr <- matrix(nrow=n, ncol=n)
119 index <- attr(object, "index")
122 matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
126 corMatrix.corGrafen <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
137 matr <- matrix(nrow=n, ncol=n)
138 index <- attr(object, "index")
141 matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
145 coef.corBrownian <- function(object, unconstrained = TRUE, ...)
147 if (!("corBrownian" %in% class(object)))
148 stop("object is not of class \"corBrownian\".")
152 coef.corMartins <- function(object, unconstrained = TRUE, ...)
154 if (!("corMartins" %in% class(object)))
155 stop("object is not of class \"corMartins\".")
157 if(attr(object, "fixed")) {
160 return(as.vector(object))
163 aux <- as.vector(object)
164 names(aux) <- "alpha"
168 coef.corGrafen <- function(object, unconstrained = TRUE, ...)
170 if (!("corGrafen" %in% class(object)))
171 stop("object is not of class \"corGrafen\".")
173 if(attr(object, "fixed")) {
176 return(as.vector(object))
179 aux <- exp(as.vector(object))
184 ### removed node.sons() and node.leafnumber() (2006-10-12)
186 ### changed by EP (2006-10-12):
188 compute.brlen <- function(phy, method = "Grafen", power = 1, ...)
190 if (!inherits(phy, "phylo"))
191 stop('object "phy" is not of class "phylo"')
192 Ntip <- length(phy$tip.label)
194 Nedge <- dim(phy$edge)[1]
195 if (is.numeric(method)) {
196 phy$edge.length <- rep(method, length.out = Nedge)
199 if (is.function(method)) {
200 phy$edge.length <- method(Nedge, ...)
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
211 (xx[phy$edge[, 1]]/m)^power - (xx[phy$edge[, 2]]/m)^power
218 corPagel <- function(value, phy, form = ~1, fixed = FALSE)
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")
231 corMatrix.corPagel <-
232 function(object, covariate = getCovariate(object), corr = TRUE, ...)
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]
245 coef.corPagel <- function(object, unconstrained = TRUE, ...)
248 if (attr(object, "fixed")) return(numeric(0))
249 else return(object[1])
252 names(aux) <- "lambda"
256 corBlomberg <- function(value, phy, form = ~1, fixed = FALSE)
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")
269 corMatrix.corBlomberg <-
270 function(object, covariate = getCovariate(object), corr = TRUE, ...)
272 index <- attr(object, "index")
274 stop("object has not been initialized")
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)
285 coef.corBlomberg <- function(object, unconstrained = TRUE, ...)
288 if (attr(object, "fixed")) return(numeric(0))
289 else return(object[1])