X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FPGLS.R;h=6b6f4afe05408b8a1885f2ac76aa904edf319234;hb=e9f75370481c37eb9d3d811ce7494f818b423136;hp=1b918e9a9574097322794c5ff5b17fd9b5277638;hpb=048cc125137aa5f9270ce02207ccaa8be02e53a8;p=ape.git

diff --git a/R/PGLS.R b/R/PGLS.R
index 1b918e9..6b6f4af 100644
--- a/R/PGLS.R
+++ b/R/PGLS.R
@@ -1,4 +1,4 @@
-## PGLS.R (2008-04-18)
+## PGLS.R (2008-05-08)
 
 ##   Phylogenetic Generalized Least Squares
 
@@ -9,7 +9,8 @@
 
 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
@@ -19,9 +20,11 @@ corBrownian <- function(value = 1, phy, form=~1)
 
 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
@@ -31,10 +34,12 @@ corMartins <- function(value, phy, form=~1, fixed=FALSE)
 
 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
@@ -81,8 +86,10 @@ Initialize.corPhyl <- function(object, data, ...)
 
 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]
@@ -95,49 +102,57 @@ corMatrix.corBrownian <- function(object, covariate = getCovariate(object), corr
   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))
@@ -152,7 +167,8 @@ coef.corMartins <- function(object, unconstrained = TRUE, ...)
 
 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))
@@ -171,8 +187,8 @@ coef.corGrafen <- function(object, unconstrained = TRUE, ...)
 
 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]
@@ -197,7 +213,7 @@ compute.brlen <- function(phy, method = "Grafen", power = 1, ...)
     }
 }
 
-### by EP:
+## by EP:
 
 corPagel <- function(value, phy, form = ~1, fixed = FALSE)
 {
@@ -236,3 +252,43 @@ coef.corPagel <- function(object, unconstrained = TRUE, ...)
     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
+}