X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2FPGLS.R;h=5378ed109e8112aec4938ae12c6eae148dfc46f6;hb=756ad52b92dc1ac2922cf62ce882469ad4cacc2c;hp=8c5ab74032e3ad75e5b44d008270834c7e5023d6;hpb=bab1a9f9a25774bb2c11430540b02f1b91f458bc;p=ape.git

diff --git a/R/PGLS.R b/R/PGLS.R
index 8c5ab74..5378ed1 100644
--- a/R/PGLS.R
+++ b/R/PGLS.R
@@ -1,168 +1,173 @@
-## PGLS.R (2008-04-18)
+## PGLS.R (2012-02-09)
 
 ##   Phylogenetic Generalized Least Squares
 
-## Copyright 2004 Julien Dutheil, and 2006-2008 Emmanuel Paradis
+## Copyright 2004 Julien Dutheil, and 2006-2012 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
-corBrownian <- function(value = 1, phy, form=~1)
+corBrownian <- function(value = 1, phy, form = ~1)
 {
-  if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
-  attr(value, "formula") <- form
-  attr(value, "fixed")   <- TRUE
-  attr(value, "tree")    <- phy
-  class(value) <- c("corBrownian", "corPhyl", "corStruct")
-  return(value)
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
+    attr(value, "formula") <- form
+    attr(value, "fixed") <- TRUE
+    attr(value, "tree") <- phy
+    class(value) <- c("corBrownian", "corPhyl", "corStruct")
+    value
 }
 
-corMartins <- function(value, phy, form=~1, fixed=FALSE)
+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\"")
-  attr(value, "formula") <- form
-  attr(value, "fixed")   <- fixed
-  attr(value, "tree")    <- phy
-  class(value) <- c("corMartins", "corPhyl", "corStruct")
-  return(value)
+    if (length(value) > 1)
+        stop("only one parameter is allowed")
+    if (value < 0) stop("the parameter alpha must be positive")
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
+    attr(value, "formula") <- form
+    attr(value, "fixed") <- fixed
+    attr(value, "tree") <- phy
+    class(value) <- c("corMartins", "corPhyl", "corStruct")
+    value
 }
 
-corGrafen <- 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.")
-  value <- log(value) # Optimization under constraint, use exponential transform.
-  if (!("phylo" %in% class(phy))) stop("ERROR!!! Object \"phy\" is not of class \"phylo\"")
-  attr(value, "formula") <- form
-  attr(value, "fixed")   <- fixed
-  attr(value, "tree")    <- phy
-  class(value) <- c("corGrafen", "corPhyl", "corStruct")
-  return(value)
+    if (length(value) > 1)
+        stop("only one parameter is allowed")
+    if (value < 0) stop("parameter rho must be positive")
+    value <- log(value) # Optimization under constraint, use exponential transform.
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo"')
+    attr(value, "formula") <- form
+    attr(value, "fixed") <- fixed
+    attr(value, "tree") <- phy
+    class(value) <- c("corGrafen", "corPhyl", "corStruct")
+    value
 }
 
 Initialize.corPhyl <- function(object, data, ...)
 {
-  # The same as in Initialize corStruct:
-  form <- formula(object)
-  ## Obtaining the group information, if any
-  if(!is.null(getGroupsFormula(form))) {
-    attr(object, "groups") <- getGroups(object, form, data = data)
-    attr(object, "Dim")    <- Dim(object, attr(object, "groups"))
-  } else { # no groups
-    attr(object, "Dim")    <- Dim(object, as.factor(rep(1, nrow(data))))
-  }
-  ## Obtaining the covariate(s)
-  attr(object, "covariate") <- getCovariate(object, data = data)
-
-  # Specific to corPhyl:
-  phy <- attr(object, "tree")
-  if (is.null(data))
-    data <- parent.frame()
-  ## Added by EP 29 May 2006:
-  if (nrow(data) != length(phy$tip.label))
-    stop("number of observations and number of tips in the tree are not equal.")
-  ## END
-  if(is.null(rownames(data))) {
-    warning("No row names supplied in dataframe, data taken to be in the same order as in tree.")
-    attr(object, "index") <- 1:dim(data)[1]
-  } else {
-    index <- match(rownames(data), phy$tip.label)
-    if(any(is.na(index))) {
-      warning("Row names in dataframe do not match tree tip names. data taken to be in the same order as in tree.")
-      attr(object, "index") <- 1:dim(data)[1]
+    ## The same as in Initialize corStruct:
+    form <- formula(object)
+    ## Obtaining the group information, if any
+    if (!is.null(getGroupsFormula(form))) {
+        attr(object, "groups") <- getGroups(object, form, data = data)
+        attr(object, "Dim") <- Dim(object, attr(object, "groups"))
+    } else { # no groups
+        attr(object, "Dim") <- Dim(object, as.factor(rep(1, nrow(data))))
+    }
+    ## Obtaining the covariate(s)
+    attr(object, "covariate") <- getCovariate(object, data = data)
+
+    ## Specific to corPhyl:
+    phy <- attr(object, "tree")
+    if (is.null(data)) data <- parent.frame()
+    ## Added by EP 29 May 2006:
+    if (nrow(data) != length(phy$tip.label))
+        stop("number of observations and number of tips in the tree are not equal.")
+    ## END
+    if (is.null(rownames(data))) {
+        warning("No rownames supplied in data frame, data taken to be in the same order than in tree")
+        attr(object, "index") <- 1:dim(data)[1]
     } else {
-      attr(object, "index") <- index
+        index <- match(rownames(data), phy$tip.label)
+        if (any(is.na(index))) {
+            warning("Rownames in data frame do not match tree tip names; data taken to be in the same order as in tree")
+            attr(object, "index") <- 1:dim(data)[1]
+        } else {
+            attr(object, "index") <- index
+        }
     }
-  }
-  return(object)
-}
-
-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.")
-  tree <- attr(object, "tree")
-  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]]
-  return(matr)
-}
-
-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)
-}
-
-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]))
-  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]]
-  return(matr)
+    object
+}
+
+corMatrix.corBrownian <-
+    function(object, covariate = getCovariate(object), corr = TRUE, ...)
+{
+    if (!("corBrownian" %in% class(object)))
+        stop('object is not of class "corBrownian"')
+    if (!any(attr(object, "index")))
+        stop("object has not been initialized.")
+    tree <- attr(object, "tree")
+    mat <- vcv.phylo(tree, corr = corr)
+    n <- dim(mat)[1]
+    ## reorder matrix:
+    index <- attr(object, "index")
+    mat[index, index]
+}
+
+corMatrix.corMartins <-
+    function(object, covariate = getCovariate(object), corr = TRUE, ...)
+{
+    if (!("corMartins" %in% class(object)))
+        stop('object is not of class "corMartins"')
+    if (!any(attr(object, "index")))
+        stop("object has 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:
+    index <- attr(object, "index")
+    mat[index, index]
+}
+
+corMatrix.corGrafen <-
+    function(object, covariate = getCovariate(object), corr = TRUE, ...)
+{
+    if (!("corGrafen" %in% class(object)))
+        stop('object is not of class "corGrafen"')
+    if (!any(attr(object, "index")))
+        stop("object has not been initialized.")
+    tree <- compute.brlen(attr(object, "tree"),
+                          method = "Grafen", power = exp(object[1]))
+    mat <- vcv.phylo(tree, corr = corr)
+    n <- dim(mat)[1]
+    ## reorder matrix:
+    index <- attr(object, "index")
+    mat[index, index]
 }
 
 coef.corBrownian <- function(object, unconstrained = TRUE, ...)
 {
-  if (!("corBrownian" %in% class(object))) stop("ERROR!!! Object is not of class \"corBrownian\".")
-  return(numeric(0))
+    if (!("corBrownian" %in% class(object)))
+        stop('object is not of class "corBrownian"')
+    numeric(0)
 }
 
 coef.corMartins <- function(object, unconstrained = TRUE, ...)
 {
-  if (!("corMartins" %in% class(object))) stop("ERROR!!! Object is not of class \"corMartins\".")
-  if(unconstrained) {
-    if(attr(object, "fixed")) {
-      return(numeric(0))
-    } else {
-      return(as.vector(object))
+    if (!("corMartins" %in% class(object)))
+        stop('object is not of class "corMartins"')
+    if (unconstrained) {
+        if (attr(object, "fixed")) {
+            return(numeric(0))
+        } else {
+            return(as.vector(object))
+        }
     }
-  }
-  aux <- as.vector(object)
-  names(aux) <- "alpha"
-  return(aux)
+    aux <- as.vector(object)
+    names(aux) <- "alpha"
+    aux
 }
 
 coef.corGrafen <- function(object, unconstrained = TRUE, ...)
 {
-  if (!("corGrafen" %in% class(object))) stop("ERROR!!! Object is not of class \"corGrafen\".")
-  if(unconstrained) {
-    if(attr(object, "fixed")) {
-      return(numeric(0))
-    } else {
-      return(as.vector(object))
+    if (!("corGrafen" %in% class(object)))
+        stop('object is not of class "corGrafen"')
+    if (unconstrained) {
+        if (attr(object, "fixed")) {
+            return(numeric(0))
+        } else {
+            return(as.vector(object))
+        }
     }
-  }
-  aux <- exp(as.vector(object))
-  names(aux) <- "rho"
-  return(aux)
+    aux <- exp(as.vector(object))
+    names(aux) <- "rho"
+    aux
 }
 
 ### removed node.sons() and node.leafnumber()  (2006-10-12)
@@ -171,8 +176,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]
@@ -203,7 +208,7 @@ corPagel <- function(value, phy, form = ~1, fixed = FALSE)
 {
     if (value < 0 || value > 1)
         stop("the value of lambda must be between 0 and 1.")
-    if (!("phylo" %in% class(phy)))
+    if (!inherits(phy, "phylo"))
         stop('object "phy" is not of class "phylo"')
     attr(value, "formula") <- form
     attr(value, "fixed") <- fixed
@@ -217,7 +222,7 @@ corMatrix.corPagel <-
 {
     if (!any(attr(object, "index")))
         stop("object has not been initialized")
-    mat <- vcv.phylo(attr(object, "tree"), cor = corr)
+    mat <- vcv.phylo(attr(object, "tree"), corr = corr)
     index <- attr(object, "index")
     mat <- mat[index, index]
     tmp <- diag(mat)
@@ -236,3 +241,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 (!inherits(phy, "phylo"))
+        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, corr = 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
+}