-## PGLS.R (2008-04-18)
+## PGLS.R (2008-05-08)
## Phylogenetic Generalized Least Squares
corBrownian <- function(value = 1, phy, form=~1)
{
- if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
+ if (!("phylo" %in% class(phy)))
+ stop("object \"phy\" is not of class \"phylo\"")
attr(value, "formula") <- form
attr(value, "fixed") <- TRUE
attr(value, "tree") <- phy
corMartins <- function(value, phy, form=~1, fixed=FALSE)
{
- if(length(value) > 1) stop("ERROR!!! Only one parameter is allowed in corPGLS structure.")
- if(value < 0) stop("ERROR!!! Parameter alpha must be positive.")
- if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
+ if(length(value) > 1)
+ stop("only one parameter is allowed in corPGLS structure.")
+ if(value < 0) stop("parameter alpha must be positive.")
+ if (!("phylo" %in% class(phy)))
+ stop("object \"phy\" is not of class \"phylo\"")
attr(value, "formula") <- form
attr(value, "fixed") <- fixed
attr(value, "tree") <- phy
corGrafen <- function(value, phy, form=~1, fixed=FALSE)
{
- if(length(value) > 1) stop("ERROR!!! Only one parameter is allowed in corGrafen structure.")
- if(value < 0) stop("ERROR!!! Parameter rho must be positive.")
+ if(length(value) > 1)
+ stop("only one parameter is allowed in corGrafen structure.")
+ if(value < 0) stop("parameter rho must be positive.")
value <- log(value) # Optimization under constraint, use exponential transform.
- if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
+ if (!("phylo" %in% class(phy)))
+ stop("object \"phy\" is not of class \"phylo\"")
attr(value, "formula") <- form
attr(value, "fixed") <- fixed
attr(value, "tree") <- phy
corMatrix.corBrownian <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
- if (!("corBrownian" %in% class(object))) stop("ERROR!!! Object is not of class \"corBrownian\".")
- if(!any(attr(object, "index"))) stop("ERROR!!! object have not been initialized.")
+ if (!("corBrownian" %in% class(object)))
+ stop("object is not of class \"corBrownian\".")
+ if(!any(attr(object, "index")))
+ stop("object have not been initialized.")
tree <- attr(object, "tree")
mat <- vcv.phylo(tree, cor = corr)
n <- dim(mat)[1]
return(matr)
}
-corMatrix.corMartins <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
+corMatrix.corMartins <-
+ function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
- if (!("corMartins" %in% class(object))) stop("ERROR!!! Object is not of class \"corMartins\".")
- if(!any(attr(object, "index"))) stop("ERROR!!! object have not been initialized.")
- tree <- attr(object, "tree")
- dist <- cophenetic.phylo(tree)
- mat <- exp(-object[1] * dist)
- if(corr) mat <- cov2cor(mat)
- n <- dim(mat)[1]
- # reorder matrix:
- matr <- matrix(nrow=n, ncol=n)
- index <- attr(object, "index")
- for(i in 1:n)
- for(j in i:n)
- matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
- return(matr)
+ if (!("corMartins" %in% class(object)))
+ stop("object is not of class \"corMartins\".")
+ if(!any(attr(object, "index")))
+ stop("object have not been initialized.")
+ tree <- attr(object, "tree")
+ dist <- cophenetic.phylo(tree)
+ mat <- exp(-object[1] * dist)
+ if (corr) mat <- cov2cor(mat)
+ n <- dim(mat)[1]
+ ## reorder matrix:
+ matr <- matrix(nrow=n, ncol=n)
+ index <- attr(object, "index")
+ for(i in 1:n)
+ for(j in i:n)
+ matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
+ matr
}
corMatrix.corGrafen <- function(object, covariate = getCovariate(object), corr = TRUE, ...)
{
- if (!("corGrafen" %in% class(object))) stop("ERROR!!! Object is not of class \"corGrafen\".")
- if(!any(attr(object, "index"))) stop("ERROR!!! object have not been initialized.")
- tree <- compute.brlen(attr(object, "tree"), method = "Grafen", power = exp(object[1]))
+ if (!("corGrafen" %in% class(object)))
+ stop("object is not of class \"corGrafen\".")
+ if (!any(attr(object, "index")))
+ stop("object have not been initialized.")
+ tree <- compute.brlen(attr(object, "tree"),
+ method = "Grafen", power = exp(object[1]))
mat <- vcv.phylo(tree, cor = corr)
n <- dim(mat)[1]
# reorder matrix:
matr <- matrix(nrow=n, ncol=n)
index <- attr(object, "index")
for(i in 1:n)
- for(j in i:n)
- matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
+ for(j in i:n)
+ matr[i,j] <- matr[j,i] <- mat[index[i], index[j]]
return(matr)
}
coef.corBrownian <- function(object, unconstrained = TRUE, ...)
{
- if (!("corBrownian" %in% class(object))) stop("ERROR!!! Object is not of class \"corBrownian\".")
+ if (!("corBrownian" %in% class(object)))
+ stop("object is not of class \"corBrownian\".")
return(numeric(0))
}
coef.corMartins <- function(object, unconstrained = TRUE, ...)
{
- if (!("corMartins" %in% class(object))) stop("ERROR!!! Object is not of class \"corMartins\".")
+ if (!("corMartins" %in% class(object)))
+ stop("object is not of class \"corMartins\".")
if(unconstrained) {
if(attr(object, "fixed")) {
return(numeric(0))
coef.corGrafen <- function(object, unconstrained = TRUE, ...)
{
- if (!("corGrafen" %in% class(object))) stop("ERROR!!! Object is not of class \"corGrafen\".")
+ if (!("corGrafen" %in% class(object)))
+ stop("object is not of class \"corGrafen\".")
if(unconstrained) {
if(attr(object, "fixed")) {
return(numeric(0))
compute.brlen <- function(phy, method = "Grafen", power = 1, ...)
{
- if (!"phylo" %in% class(phy))
- stop('object "phy" is not of class "phylo"')
+ if (!inherits(phy, "phylo"))
+ stop('object "phy" is not of class "phylo"')
Ntip <- length(phy$tip.label)
Nnode <- phy$Nnode
Nedge <- dim(phy$edge)[1]
}
}
-### by EP:
+## by EP:
corPagel <- function(value, phy, form = ~1, fixed = FALSE)
{
names(aux) <- "lambda"
aux
}
+
+corBlomberg <- function(value, phy, form = ~1, fixed = FALSE)
+{
+ if (value <= 0)
+ stop("the value of g must be greater than 0.")
+ if (!("phylo" %in% class(phy)))
+ stop('object "phy" is not of class "phylo"')
+ attr(value, "formula") <- form
+ attr(value, "fixed") <- fixed
+ attr(value, "tree") <- phy
+ class(value) <- c("corBlomberg", "corPhyl", "corStruct")
+ value
+}
+
+corMatrix.corBlomberg <-
+ function(object, covariate = getCovariate(object), corr = TRUE, ...)
+{
+ index <- attr(object, "index")
+ if (is.null(index))
+ stop("object has not been initialized")
+ if (object[1] <= 0)
+ stop("the optimization has reached a value <= 0 for parameter 'g':
+probably need to set 'fixed = TRUE' in corBlomberg().")
+ phy <- attr(object, "tree")
+ d <- (dist.nodes(phy)[length(phy$tip.label) + 1, ])^(1/object[1])
+ phy$edge.length <- d[phy$edge[, 2]] - d[phy$edge[, 1]]
+ mat <- vcv.phylo(phy, cor = corr)
+ mat[index, index]
+}
+
+coef.corBlomberg <- function(object, unconstrained = TRUE, ...)
+{
+ if (unconstrained) {
+ if (attr(object, "fixed")) return(numeric(0))
+ else return(object[1])
+ }
+ aux <- object[1]
+ names(aux) <- "g"
+ aux
+}