[Picante-commits] r48 - pkg/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Apr 11 06:19:32 CEST 2008
Author: pdc
Date: 2008-04-11 06:19:31 +0200 (Fri, 11 Apr 2008)
New Revision: 48
Modified:
pkg/R/pic.circular.R
Log:
consider the direction of shift around the circular trait
Modified: pkg/R/pic.circular.R
===================================================================
--- pkg/R/pic.circular.R 2008-04-02 02:22:00 UTC (rev 47)
+++ pkg/R/pic.circular.R 2008-04-11 04:19:31 UTC (rev 48)
@@ -1,7 +1,7 @@
### pic.R (2008-02-02)
###
-### Phylogenetically Independent Contrasts on circular data
-### an extension of the ape pic function
+### Phylogenetically Independent Contrasts on circular data
+### an extension of the ape pic function
###
### Copyright 2008 Peter Cowan
###
@@ -10,7 +10,7 @@
### Orginal copyright
### pic.R (2006-10-29)
###
-### Phylogenetically Independent Contrasts
+### Phylogenetically Independent Contrasts
###
### Copyright 2002-2006 Emmanuel Paradis
###
@@ -20,87 +20,84 @@
`pic.circular` <-
function(x, phy, scaled = TRUE, var.contrasts = FALSE)
{
- if (!require(circular)) {stop("The 'circular' package is required for this function")}
- if (class(phy) != "phylo")
- stop("object 'phy' is not of class \"phylo\"")
- if (is.null(phy$edge.length))
- stop("your tree has no branch lengths: you may consider setting them equal
- to one, or using the function `compute.brlen'.")
- nb.tip <- length(phy$tip.label)
- nb.node <- phy$Nnode
- if (nb.node != nb.tip - 1)
- stop("'phy' is not rooted and fully dichotomous")
- if (length(x) != nb.tip)
- stop("length of phenotypic and of phylogenetic data do not match")
- if (any(is.na(x)))
- stop("the present method cannot (yet) be used directly with missing data:
- you may consider removing the species with missing data from your
- tree with the function `drop.tip'.")
- if(!is.circular(x))
- stop("This version for circular data only")
- phy <- reorder(phy, "pruningwise")
- phenotype <- as.circular(rep(NA, nb.tip + nb.node),
- units = 'radians', modulo = '2pi',
- rotation = 'counter', zero = 0, type = 'angles',
- template = 'none'
- )
+ if (!require(circular)) {stop("The 'circular' package is required for this function")}
+ if (class(phy) != "phylo")
+ stop("object 'phy' is not of class \"phylo\"")
+ if (is.null(phy$edge.length))
+ stop("your tree has no branch lengths: you may consider setting them equal
+ to one, or using the function `compute.brlen'.")
+ nb.tip <- length(phy$tip.label)
+ nb.node <- phy$Nnode
+ if (nb.node != nb.tip - 1)
+ stop("'phy' is not rooted and fully dichotomous")
+ if (length(x) != nb.tip)
+ stop("length of phenotypic and of phylogenetic data do not match")
+ if (any(is.na(x)))
+ stop("the present method cannot (yet) be used directly with missing data:
+ you may consider removing the species with missing data from your
+ tree with the function `drop.tip'.")
+ if(!is.circular(x))
+ stop("This version for circular data only")
+ phy <- reorder(phy, "pruningwise")
+ phenotype <- as.circular(rep(NA, nb.tip + nb.node),
+ units = 'radians', modulo = '2pi',
+ rotation = 'counter', zero = 0, type = 'angles',
+ template = 'none'
+ )
+ if (all(names(x) %in% phy$tip.label)) {
+ phenotype[1:nb.tip] <- x[phy$tip.label]
+ } else {
+ phenotype[1:nb.tip] <- x
+ warning('the names of argument "x" and the names of the tip labels
+ did not match: the former were ignored in the analysis.')
+ }
- if (is.null(names(x))) {
- phenotype[1:nb.tip] <- x
- } else {
- if (all(names(x) %in% phy$tip.label))
- phenotype[1:nb.tip] <- x[phy$tip.label]
- else {
- phenotype[1:nb.tip] <- x
- warning('the names of argument "x" and the names of the tip labels
- did not match: the former were ignored in the analysis.')
- }
- }
- ## No need to copy the branch lengths: they are rescaled
- ## in the C code, so it's important to leave the default
- ## `DUP = TRUE' of .C.
- contr <- var.con <- numeric(nb.node)
+ contr <- var.con <- numeric(nb.node)
- ## ans <- .C("pic", as.integer(nb.tip), as.integer(nb.node),
- ## as.integer(phy$edge[, 1]), as.integer(phy$edge[, 2]),
- ## as.double(phy$edge.length), as.double(phenotype),
- ## as.double(contr), as.double(var.con),
- ## as.integer(var.contrasts), as.integer(scaled),
- ## PACKAGE = "ape")
bl=phy$edge.length
- ## The "old" R code:
- for (i in seq(from = 1, by = 2, length.out = nb.node)) {
- j <- i + 1
- anc <- phy$edge[i, 1]
- des1 <- phy$edge[i, 2]
- des2 <- phy$edge[j, 2]
- sumbl <- bl[i] + bl[j]
- ic <- anc - nb.tip
- contr[ic] <- range.circular(c(phenotype[des1], phenotype[des2]))
- if (scaled) contr[ic] <- contr[ic]/sqrt(sumbl)
- if (var.contrasts) var.con[ic] <- sumbl
-
- sin_des1 <- sin(phenotype[des1]) * bl[j]
- cos_des1 <- cos(phenotype[des1]) * bl[j]
-
- sin_des2 <- sin(phenotype[des2]) * bl[i]
- cos_des2 <- cos(phenotype[des2]) * bl[i]
- phenotype[anc] <- as.circular(atan2(sin_des1 + sin_des2, cos_des1 + cos_des2),
- type = "angles", units = "radians", template = "none",
- rotation = "counter", zero = 0, modulo = "2pi"
- )
+ for (i in seq(from = 1, by = 2, length.out = nb.node)) {
+ j <- i + 1
+ anc <- phy$edge[i, 1]
+ des1 <- phy$edge[i, 2]
+ des2 <- phy$edge[j, 2]
+ sumbl <- bl[i] + bl[j]
+ ic <- anc - nb.tip
+ tempcontr <- phenotype[des1] - phenotype[des2]
+ abtemp <- abs(tempcontr)
+ if(abtemp > (2 * pi)) {
+ stop("ERROR. contrast between ", substitute(phenotype[des1] ),
+ ' and ', substitute(phenotype[des2]), "is", substitute(tempcontr))
+ }
+ if(abtemp > pi) {
+ if(phenotype[des1] > phenotype[des2]) {
+ tempcontr <- tempcontr - (2 * pi)
+ } else if(phenotype[des1] < phenotype[des2]) {
+ tempcontr <- tempcontr + (2 * pi)
+ } else {
+ stop("Fatal Error", substitute(phenotype[des1] ),
+ ' and ', substitute(phenotype[des2]), " is ", substitute(tempcontr))}
+ }
+ contr[ic] <- tempcontr
+ if (scaled) contr[ic] <- contr[ic]/sqrt(sumbl)
+ if (var.contrasts) var.con[ic] <- sumbl
+
+ sin_des1 <- sin(as.numeric(phenotype[des1])) * bl[j]
+ cos_des1 <- cos(as.numeric(phenotype[des1])) * bl[j]
- k <- which(phy$edge[, 2] == anc)
- bl[k] <- bl[k] + bl[i]*bl[j]/sumbl
-
- }
- # contr <- ans[[7]]
- # if (var.contrasts) {
- # contr <- cbind(contr, ans[[8]])
- # dimnames(contr) <- list(1:nb.node + nb.tip, c("contrasts", "variance"))
- # } else
- # TODO check that the var.contrasts = T results are as expected
- names(contr) <- 1:nb.node + nb.tip
- contr
+ sin_des2 <- sin(as.numeric(phenotype[des2])) * bl[i]
+ cos_des2 <- cos(as.numeric(phenotype[des2])) * bl[i]
+
+ phenotype[anc] <- as.circular(atan2(sin_des1 + sin_des2, cos_des1 + cos_des2),
+ type = "angles", units = "radians", template = "none",
+ rotation = "counter", zero = 0, modulo = "2pi"
+ )
+ k <- which(phy$edge[, 2] == anc)
+ bl[k] <- bl[k] + bl[i]*bl[j]/sumbl
+
+ }
+
+ # TODO check that the var.contrasts = T results are as expected
+ names(contr) <- 1:nb.node + nb.tip
+ contr
}
More information about the Picante-commits
mailing list