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))) 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")
20 corMartins <- function(value, phy, form=~1, fixed=FALSE)
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")
32 corGrafen <- function(value, phy, form=~1, fixed=FALSE)
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")
45 Initialize.corPhyl <- function(object, data, ...)
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"))
54 attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
56 ## Obtaining the covariate(s)
57 attr(object, "covariate") <- getCovariate(object, data = data)
59 # Specific to corPhyl:
60 phy <- attr(object, "tree")
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.")
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]
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]
76 attr(object, "index") <- index
82 corMatrix.corBrownian <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
90 matr <- matrix(nrow=n, ncol=n)
91 index <- attr(object, "index")
94 matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
98 corMatrix.corMartins <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
108 matr <- matrix(nrow=n, ncol=n)
109 index <- attr(object, "index")
112 matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
116 corMatrix.corGrafen <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
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)
124 matr <- matrix(nrow=n, ncol=n)
125 index <- attr(object, "index")
128 matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
132 coef.corBrownian <- function(object, unconstrained = TRUE, ...)
134 if (!("corBrownian" %in% class(object))) stop("ERROR!!! Object is not of class \"corBrownian\".")
138 coef.corMartins <- function(object, unconstrained = TRUE, ...)
140 if (!("corMartins" %in% class(object))) stop("ERROR!!! Object is not of class \"corMartins\".")
142 if(attr(object, "fixed")) {
145 return(as.vector(object))
148 aux <- as.vector(object)
149 names(aux) <- "alpha"
153 coef.corGrafen <- function(object, unconstrained = TRUE, ...)
155 if (!("corGrafen" %in% class(object))) stop("ERROR!!! Object is not of class \"corGrafen\".")
157 if(attr(object, "fixed")) {
160 return(as.vector(object))
163 aux <- exp(as.vector(object))
168 ### removed node.sons() and node.leafnumber() (2006-10-12)
170 ### changed by EP (2006-10-12):
172 compute.brlen <- function(phy, method = "Grafen", power = 1, ...)
174 if (!"phylo" %in% class(phy))
175 stop('object "phy" is not of class "phylo"')
176 Ntip <- length(phy$tip.label)
178 Nedge <- dim(phy$edge)[1]
179 if (is.numeric(method)) {
180 phy$edge.length <- rep(method, length.out = Nedge)
183 if (is.function(method)) {
184 phy$edge.length <- method(Nedge, ...)
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
195 (xx[phy$edge[, 1]]/m)^power - (xx[phy$edge[, 2]]/m)^power
202 corPagel <- function(value, phy, form = ~1, fixed = FALSE)
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")
215 corMatrix.corPagel <-
216 function(object, covariate = getCovariate(object), corr = TRUE, ...)
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]
229 coef.corPagel <- function(object, unconstrained = TRUE, ...)
232 if (attr(object, "fixed")) return(numeric(0))
233 else return(object[1])
236 names(aux) <- "lambda"
240 corBlomberg <- function(value, phy, form = ~1, fixed = FALSE)
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")
253 corMatrix.corBlomberg <-
254 function(object, covariate = getCovariate(object), corr = TRUE, ...)
256 index <- attr(object, "index")
258 stop("object has not been initialized")
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)
269 coef.corBlomberg <- function(object, unconstrained = TRUE, ...)
272 if (attr(object, "fixed")) return(numeric(0))
273 else return(object[1])