small changes in PGLS.R

git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@124 6e262413-ae40-0410-9e79-b911bd7a66b7

diff --git a/ChangeLog b/ChangeLog
index d1862f60b567b5ac81a855dbc0ddb53a86a9d0a0..8df6cc8b7aa06925ce5d54d8928a1742363128c8 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -7,8 +7,8 @@ NEW FEATURES
       text to be plotted in different fonts.
 
     o There are now replacement operators for [, [[, and $ for the class
-      "multiPhylo" (i.e., TREES[11:20] <- rmtree(10, 100)). They check
-      that the tip labels are the same in all trees.
+      "multiPhylo" (i.e., TREES[11:20] <- rmtree(10, 100)). They possibly
+      check that the tip labels are the same in all trees.
 
     o Objects of class "multiPhylo" can be built with c(): there are
       methods for the classes "phylo" and "multiPhylo".

diff --git a/R/PGLS.R b/R/PGLS.R
index df3e6b1d2601d984a3809ecacea3ac948423af42..510f4017ebb857a51eaba2275d45c3c329deac3f 100644 (file)
--- a/R/PGLS.R
+++ b/R/PGLS.R
@@ -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))