Package: ape
Version: 2.7-3
-Date: 2011-06-21
+Date: 2011-06-22
Title: Analyses of Phylogenetics and Evolution
Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong, Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel, Christoph Heibl, Daniel Lawson, Vincent Lefort, Pierre Legendre, Jim Lemon, Yvonnick Noel, Johan Nylander, Rainer Opgen-Rhein, Klaus Schliep, Korbinian Strimmer, Damien de Vienne
Maintainer: Emmanuel Paradis <Emmanuel.Paradis@ird.fr>
CHANGES IN APE VERSION 2.7-3
+NEW FEATURES
+
+ o mantel.test() has a new argument 'alternative' which is
+ "two-sided" by default. Previously, this test was one-tailed
+ with no possibility to change.
+
+
BUG FIXES
o Branch lengths were wrongly updated with bind.tree(, where = <tip>,
position = 0). (Thanks to Liam Revell for digging this bug out.)
+ o Simulation of OU process with rTraitCont() did not work correctly.
+ This now uses formula from Gillespie (1996) reduced to a BM
+ process when alpha = 0 avoiding division by zero.
+
CHANGES IN APE VERSION 2.7-2
-## mantel.test.R (2006-07-28)
+## mantel.test.R (2011-06-22)
## Mantel Test for Similarity of Two Matrices
m[col(m) <= row(m)]
}
-mantel.test <- function (m1, m2, nperm = 1000, graph = FALSE, ...)
+mantel.test <- function (m1, m2, nperm = 1000, graph = FALSE,
+ alternative = "two.sided", ...)
{
+ alternative <- match.arg(alternative, c("two.sided", "less", "greater"))
n <- nrow(m1)
realz <- mant.zstat(m1, m2)
nullstats <- replicate(nperm, mant.zstat(m1, perm.rowscols(m2, n)))
- pval <- sum(nullstats > realz)/nperm
+ pval <- switch(alternative,
+ "two.sided" = 2 * sum(abs(nullstats) > abs(realz)),
+ "less" = sum(nullstats < realz),
+ "greater" = sum(nullstats > realz))
+ pval <- pval / nperm
if (graph) {
plot(density(nullstats), type = "l", ...)
abline(v = realz)
}
- list(z.stat = realz, p = pval)
+ list(z.stat = realz, p = pval, alternative = alternative)
}
-## rTrait.R (2011-06-16)
+## rTrait.R (2011-06-15)
## Trait Evolution
rTraitCont <-
function(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0,
- ancestor = FALSE, root.value = 0, linear = TRUE, ...)
+ ancestor = FALSE, root.value = 0, ...)
{
if (is.null(phy$edge.length))
stop("tree has no branch length")
if (length(theta) == 1) theta <- rep(theta, N)
else if (length(theta) != N)
stop("'theta' must have one or Nedge(phy) elements")
- if (!linear) model <- model + 1L
}
.C("rTraitCont", as.integer(model), as.integer(N),
as.integer(anc - 1L), as.integer(des - 1L), el,
\alias{mantel.test}
\title{Mantel Test for Similarity of Two Matrices}
\usage{
-mantel.test(m1, m2, nperm = 1000, graph = FALSE, ...)
+mantel.test(m1, m2, nperm = 1000, graph = FALSE,
+ alternative = "two.sided", ...)
}
\arguments{
\item{m1}{a numeric matrix giving a measure of pairwise distances,
\item{nperm}{the number of times to permute the data.}
\item{graph}{a logical indicating whether to produce a summary graph
(by default the graph is not plotted).}
+ \item{alternative}{a character string defining the alternative
+ hypothesis: \code{"two.sided"} (default), \code{"less"},
+ \code{"greater"}, or any unambiguous abbreviation of these.}
\item{\dots}{further arguments to be passed to \code{plot()} (to add a
title, change the axis labels, and so on).}
}
\description{
This function computes Mantel's permutation test for similarity of two
matrices. It permutes the rows and columns of the two matrices
- randomly and calculates a Z-statistic.
+ randomly and calculates a \eqn{Z}-statistic.
}
\details{
- The function calculates a Z-statistic for the Mantel test, equal to
+ The function calculates a \eqn{Z}-statistic for the Mantel test, equal to
the sum of the pairwise product of the lower triangles of the
permuted matrices, for each permutation of rows and columns. It
- compares the permuted distribution with the Z-statistic observed for
- the actual data.
+ compares the permuted distribution with the \eqn{Z}-statistic observed
+ for the actual data.
If \code{graph = TRUE}, the functions plots the density estimate of
- the permutation distribution along with the observed Z-statistic as a
- vertical line.
+ the permutation distribution along with the observed \eqn{Z}-statistic
+ as a vertical line.
- The \code{...} argument allows the user to give further options to
+ The \code{\dots} argument allows the user to give further options to
the \code{plot} function: the title main be changed with \code{main=},
the axis labels with \code{xlab =}, and \code{ylab =}, and so on.
}
\value{
- \item{z.stat}{the Z-statistic (sum of rows*columns of lower triangle)
- of the data matrices.}
- \item{p}{P-value (quantile of the observed Z-statistic in the
- permutation distribution).}
+ \item{z.stat}{the \eqn{Z}-statistic (sum of rows*columns of lower
+ triangle) of the data matrices.}
+ \item{p}{\eqn{P}-value (quantile of the observed \eqn{Z}-statistic in
+ the permutation distribution).}
+ \item{alternative}{the alternative hypothesis.}
}
\references{
Mantel, N. (1967) The detection of disease clustering and a
Manly, B. F. J. (1986) \emph{Multivariate statistical methods: a primer.}
London: Chapman & Hall.
}
-\author{Original code in S by Ben Bolker \email{bolker@zoo.ufl.edu}, ported
+\author{
+ Original code in S by Ben Bolker \email{bolker@zoo.ufl.edu}, ported
to \R by Julien Claude \email{claude@isem.univ-montp2.fr}
}
\examples{
\title{Continuous Character Simulation}
\usage{
rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0,
- ancestor = FALSE, root.value = 0, linear = TRUE, ...)
+ ancestor = FALSE, root.value = 0, ...)
}
\arguments{
\item{phy}{an object of class \code{"phylo"}.}
values at the nodes as well (by default, only the values at the tips
are returned).}
\item{root.value}{a numeric giving the value at the root.}
- \item{linear}{a logical indicating which parameterisation of the OU
- model to use (see details).}
\item{\dots}{further arguments passed to \code{model} if it is a
function.}
}
\item{\code{"OU"}:}{an Ornstein-Uhlenbeck model is used. The above
indexing rule is used for the three parameters \code{sigma},
- \code{alpha}, and \code{theta}. This may be more interesting for the
- last one to model varying phenotypic optima.
-
- By default the following formula is used:
-
- \deqn{x_{t''} = x_{t'} - \alpha l (x_{t'} - \theta) + \sigma
- l \epsilon}{x(t'') = x(t') - alpha l (x(t') - theta) + sigma
- l epsilon}
-
- where \eqn{l (= t'' - t')} is the branch length, and \eqn{\epsilon
- \sim N(0, 1)}{\epsilon ~ N(0, 1)}. If \eqn{\alpha > 1}{alpha > 1},
- this may lead to chaotic oscillations. Thus an alternative
- parameterisation is used if \code{linear = FALSE}:
-
- \deqn{x_{t''} = x_{t'} - (1 - exp(-\alpha l)) * (x_{t'} - \theta) +
- \sigma l \epsilon}{x(t'') = x(t') - (1 - exp(-alpha l)) * (x(t') -
- theta) + sigma l epsilon}}
+ \code{alpha}, and \code{theta}. This may be interesting for the last
+ one to model varying phenotypic optima. The exact updating formula
+ from Gillespie (1996) are used which are reduced to BM formula if
+ \code{alpha = 0}.
\item{A function:}{it must be of the form \code{foo(x, l)} where
\code{x} is the trait of the ancestor and \code{l} is the branch
present, otherwise, ``Node1'', ``Node2'', etc.
}
\references{
+ Gillespie, D. T. (1996) Exact numerical simulation of the
+ Ornstein-Uhlenbeck process and its integral. \emph{Physical Review E},
+ \bold{54}, 2084--2091.
+
Paradis, E. (2006) \emph{Analyses of Phylogenetics and Evolution with
R.} New York: Springer.
}
/*sib is tree C being passed by B*/
/*D is tree below etest*/
double D_AB, D_CD, D_AC, D_BD;
- double thisWeight;
+ /* double thisWeight; */
sib = siblingEdge(etest);
left = etest->head->leftEdge;
right = etest->head->rightEdge;
//void bNNI(tree *T, double **avgDistArray, int *count)
void bNNI(tree *T, double **avgDistArray, int *count, double **D, int numSpecies)
{
- edge *e, *centerEdge;
+ edge *e; /* , *centerEdge; */
edge **edgeArray;
int *p, *location, *q;
int i,j;
-/* bipartition.c 2007-06-29 */
+/* bipartition.c 2011-06-23 */
-/* Copyright 2005-2007 Emmanuel Paradis, and 2007 R Development Core Team */
+/* Copyright 2005-2011 Emmanuel Paradis, and 2007 R Development Core Team */
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
SEXP prop_part(SEXP TREES, SEXP nbtree, SEXP keep_partitions)
{
- int i, j, k, l, KeepPartition, Ntree, Ntip, Nnode, Npart, NpartCurrent, *no;
+ int i, j, k, KeepPartition, Ntree, Ntip, Nnode, Npart, NpartCurrent, *no;
SEXP bp, ans, nbtip, nbnode, number;
PROTECT(nbtree = coerceVector(nbtree, INTSXP));
-/* delta_plot.c 2010-01-12 */
+/* delta_plot.c 2011-06-23 */
-/* Copyright 2010 Emmanuel Paradis
+/* Copyright 2010-2011 Emmanuel Paradis
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
int x, y, u, v; /* same notation than in Holland et al. 2002 */
int n = *size, nb = *nbins;
double dxy, dxu, dxv, dyu, dyv, duv, A, B, C, delta;
- int xy, xu, xv, yu, yv, uv, i, k;
+ int xy, xu, xv, yu, yv, uv, i;
for (x = 0; x < n - 3; x++) {
xy = x*n - x*(x + 1)/2; /* do NOT factorize */
-/* dist_dna.c 2011-02-18 */
+/* dist_dna.c 2011-06-23 */
/* Copyright 2005-2011 Emmanuel Paradis
double *BF, int *variance, double *var,
int *gamma, double *alpha)
{
- int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
- double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
+ int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
+ double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
L = *s;
double *BF, int *variance, double *var,
int *gamma, double *alpha)
{
- int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
- double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
+ int i1, i2, Nd, Ns1, Ns2, L, target, s1, s2;
+ double P1, P2, Q, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
PREPARE_BF_TN93
currently not available.
</FIXME> */
{
- int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4], ndim = 4, info, ipiv[16];
- double P12[16], P21[16], U[16];
+ int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4];
+ double P12[16], P21[16];
for (i1 = 1; i1 < *n; i1++) {
for (i2 = i1 + 1; i2 <= *n; i2++) {
-/* matexpo.c 2007-10-08 */
+/* matexpo.c 2011-06-23 */
-/* Copyright 2007 Emmanuel Paradis
+/* Copyright 2007-2011 Emmanuel Paradis
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
/* This function computes the exponential of a nr x nr matrix */
{
double *U, *vl, *WR, *Uinv, *WI, *work;
- int i, j, k, l, info, *ipiv, n = *nr, nc = n*n, lw = nc << 1, *ord;
+ int i, j, k, l, info, *ipiv, n = *nr, nc = n*n, lw = nc << 1;
char yes = 'V', no = 'N';
U = (double *)R_alloc(nc, sizeof(double));
work = (double *)R_alloc(lw, sizeof(double));
ipiv = (int *)R_alloc(nc, sizeof(int));
- ord = (int *)R_alloc(n, sizeof(int));
/* The matrix is not symmetric, so we use 'dgeev'.
We take the real part of the eigenvalues -> WR
-/* plot_phylo.c (2009-10-03) */
+/* plot_phylo.c (2011-06-23) */
-/* Copyright 2004-2009 Emmanuel Paradis
+/* Copyright 2004-2011 Emmanuel Paradis
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
void node_height(int *ntip, int *nnode, int *edge1, int *edge2,
int *nedge, double *yy)
{
- int i, k, n;
+ int i, n;
double S;
/* The coordinates of the tips have been already computed */
- k = 1;
S = 0;
n = 0;
for (i = 0; i < *nedge; i++) {
void node_height_clado(int *ntip, int *nnode, int *edge1, int *edge2,
int *nedge, double *xx, double *yy)
{
- int i, k, n;
+ int i, n;
double S;
node_depth(ntip, nnode, edge1, edge2, nedge, xx);
/* The coordinates of the tips have been already computed */
- k = 1;
S = 0;
n = 0;
for (i = 0; i < *nedge; i++) {
-/* rTrait.c 2010-07-26 */
+/* rTrait.c 2011-06-25 */
-/* Copyright 2010 Emmanuel Paradis */
+/* Copyright 2010-2011 Emmanuel Paradis */
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
}
break;
case 2 : for (i = *Nedge - 1; i >= 0; i--) {
+ if (alpha[i]) {
+ alphaT = alpha[i] * el[i];
+ M = exp(-alphaT);
+ S = sigma[i] * sqrt((1 - exp(-2 * alphaT))/(2 * alpha[i]));
+ } else {
+ M = 1;
+ S = sqrt(el[i]) * sigma[i];
+ }
GetRNGstate();
- x[edge2[i]] = x[edge1[i]] - alpha[i]*el[i]*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand();
- PutRNGstate();
- }
- break;
- case 3 : for (i = *Nedge - 1; i >= 0; i--) {
- GetRNGstate();
- x[edge2[i]] = x[edge1[i]] - (1 - exp(alpha[i]*el[i]))*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand();
+ x[edge2[i]] = x[edge1[i]] * M + theta[i] * (1 - M) + S * norm_rand();
PutRNGstate();
}
break;
-/* tree_build.c 2011-02-28 */
+/* tree_build.c 2011-06-23 */
/* Copyright 2008-2011 Emmanuel Paradis */
void decode_internal_edge(const char *x, int a, int b, char *lab, double *w)
{
- int i, k, co = a;
+ int co = a;
char *endstr, str[100];
while (x[co] != ':') co++;
projects / ape.git / commitdiff