-## dist.topo.R (2010-01-22)
+## dist.topo.R (2010-01-25)
## Topological Distances, Tree Bipartitions,
## Consensus Trees, and Bootstrapping Phylogenies
if (method == "score" && (is.null(x$edge.length) || is.null(y$edge.length)))
stop("trees must have branch lengths for Billera et al.'s distance.")
nx <- length(x$tip.label)
+ x <- unroot(x)
+ y <- unroot(y)
bp1 <- .Call("bipartition", x$edge, nx, x$Nnode, PACKAGE = "ape")
bp1 <- lapply(bp1, function(xx) sort(x$tip.label[xx]))
ny <- length(y$tip.label) # fix by Otto Cordero
--- /dev/null
+'lmorigin' <-
+function(formula, data=NULL, origin=TRUE, nperm=999, method=NULL, silent=FALSE)
+#
+# This program computes a multiple linear regression and performs tests
+# of significance of the equation parameters using permutations.
+#
+# origin=TRUE: the regression line can be forced through the origin. Testing
+# the significance in that case requires a special permutation procedure.
+#
+# Permutation methods: raw data or residuals of full model
+# Default method in regression through the origin: raw data
+# Default method in ordinary multiple regression: residuals of full model
+# - In ordinary multiple regression when m = 1: raw data
+#
+# Pierre Legendre, March 2009
+{
+if(!is.null(method)) method <- match.arg(method, c("raw", "residuals"))
+if(is.null(method) & origin==TRUE) method <- "raw"
+if(is.null(method) & origin==FALSE) method <- "residuals"
+if(nperm < 0) stop("Incorrect value for 'nperm'")
+
+## From the formula, find the variables and the number of observations 'n'
+toto <- lm(formula, data)
+mf <- match.call(expand.dots = FALSE)
+m <- match(c("formula", "data", "offset"), names(mf), 0)
+mf <- mf[c(1, m)]
+mf$drop.unused.levels <- TRUE
+mf[[1]] <- as.name("model.frame")
+mf <- eval(mf, parent.frame())
+var.names = colnames(mf) # Noms des variables
+y <- as.matrix(mf[,1])
+colnames(y) <- var.names[1]
+X <- as.matrix(mf[,-1])
+n <- nrow(mf)
+m <- ncol(X)
+
+a <- system.time({
+mm<- m # No. regression coefficients, possibly including the intercept
+if(m == 1) method <- "raw"
+if(nrow(X) != n) stop("Unequal number of rows in y and X")
+
+if(origin) {
+ if(!silent) cat("Regression through the origin",'\n')
+ reg <- lm(y ~ 0 + X)
+ } else {
+ if(!silent) cat("Multiple regression with estimation of intercept",'\n')
+ reg <- lm(y ~ X)
+ mm <- mm+1
+ }
+
+if(!silent) {
+ if(nperm > 0) {
+ if(method == "raw") {
+ cat("Permutation method =",method,"data",'\n')
+ } else {
+ cat("Permutation method =",method,"of full model",'\n')
+ }
+ }
+ }
+t.vec <- summary(reg)$coefficients[,3]
+p.param.t <- summary(reg)$coefficients[,4]
+df1 <- summary(reg)$fstatistic[[2]]
+df2 <- summary(reg)$fstatistic[[3]]
+F <- summary(reg)$fstatistic[[1]]
+y.res <- summary(reg)$residuals
+# b.vec <- summary(reg)$coefficients[,1]
+# r.sq <- summary(reg)$r.squared
+# adj.r.sq <- summary(reg)$adj.r.squared
+# p.param.F <- pf(F, df1, df2, lower.tail=FALSE)
+
+if(df1 < m) stop("\nCollinearity among the X variables. Check using 'lm'")
+
+# Permutation tests
+if(nperm > 0) {
+
+ nGT.F <- 1
+ nGT1.t <- rep(1,mm)
+ nGT2.t <- rep(1,mm)
+ sign.t <- sign(t.vec)
+
+ for(i in 1:nperm) # Permute raw data. Always use this method for F-test
+ {
+ if(origin) { # Regression through the origin
+ dia.bin <- diag((rbinom(n,1,0.5)*2)-1)
+ y.perm <- dia.bin %*% sample(y)
+ reg.perm <- lm(y.perm ~ 0 + X)
+ } else { # Multiple linear regression
+ y.perm <- sample(y,n)
+ reg.perm <- lm(y.perm ~ X)
+ }
+
+ # Permutation test of the F-statistic
+ F.perm <- summary(reg.perm)$fstatistic[1]
+ if(F.perm >= F) nGT.F <- nGT.F+1
+
+ # Permutation tests of the t-statistics: permute raw data
+ if(method == "raw") {
+ t.perm <- summary(reg.perm)$coefficients[,3]
+ if(nperm <= 5) cat(t.perm,'\n')
+ for(j in 1:mm) {
+ # One-tailed test in direction of sign
+ if(t.perm[j]*sign.t[j] >= t.vec[j]*sign.t[j]) nGT1.t[j] <- nGT1.t[j]+1
+ # Two-tailed test
+ if( abs(t.perm[j]) >= abs(t.vec[j]) ) nGT2.t[j] <- nGT2.t[j]+1
+ }
+ }
+ }
+
+ if(method == "residuals") {
+ # Permute residuals of full model
+ for(i in 1:nperm) {
+ if(origin) { # Regression through the origin
+ dia.bin <- diag((rbinom(n,1,0.5)*2)-1)
+ y.perm <- dia.bin %*% sample(y.res)
+ reg.perm <- lm(y.perm ~ 0 + X)
+ } else { # Multiple linear regression
+ y.perm <- sample(y.res,n)
+ reg.perm <- lm(y.perm ~ X)
+ }
+
+ # Permutation tests of the t-statistics: permute residuals
+ t.perm <- summary(reg.perm)$coefficients[,3]
+ if(nperm <= 5) cat(t.perm,'\n')
+ for(j in 1:mm) {
+ # One-tailed test in direction of sign
+ if(t.perm[j]*sign.t[j] >= t.vec[j]*sign.t[j]) nGT1.t[j] <- nGT1.t[j]+1
+ # Two-tailed test
+ if( abs(t.perm[j]) >= abs(t.vec[j]) ) nGT2.t[j] <- nGT2.t[j]+1
+ }
+ }
+ }
+ # Compute the permutational probabilities
+ p.perm.F <- nGT.F/(nperm+1)
+ p.perm.t1 <- nGT1.t/(nperm+1)
+ p.perm.t2 <- nGT2.t/(nperm+1)
+
+ ### Do not test intercept by permutation of residuals in multiple regression
+ if(!origin & method=="residuals") {
+ if(silent) { # Note: silent==TRUE in simulation programs
+ p.perm.t1[1] <- p.perm.t2[1] <- 1
+ } else {
+ p.perm.t1[1] <- p.perm.t2[1] <- NA
+ }
+ }
+ }
+
+})
+a[3] <- sprintf("%2f",a[3])
+if(!silent) cat("Computation time =",a[3]," sec",'\n')
+#
+if(nperm == 0) {
+
+ out <- list(reg=reg, p.param.t.2tail=p.param.t, p.param.t.1tail=p.param.t/2, origin=origin, nperm=nperm, var.names=var.names, call=match.call())
+
+ } else {
+
+ out <- list(reg=reg, p.param.t.2tail=p.param.t, p.param.t.1tail=p.param.t/2, p.perm.t.2tail=p.perm.t2, p.perm.t.1tail=p.perm.t1, p.perm.F=p.perm.F, origin=origin, nperm=nperm, method=method, var.names=var.names, call=match.call())
+
+ }
+#
+class(out) <- "lmorigin"
+out
+}
--- /dev/null
+'print.lmorigin' <-
+ function(x, ...)
+{
+ if(x$origin) {
+ cat("\nRegression through the origin",'\n')
+ } else {
+ cat("\nMultiple regression with estimation of intercept",'\n')
+ }
+ cat("\nCall:\n")
+ cat(deparse(x$call),'\n')
+ if(x$origin) { names <- x$var.names[-1] }
+ else { names <- c("(Intercept)",x$var.names[-1]) }
+
+ cat("\nCoefficients and parametric test results \n",'\n')
+ res <- as.data.frame(cbind(summary(x$reg)$coefficients[,1], summary(x$reg)$coefficients[,2], summary(x$reg)$coefficients[,3], summary(x$reg)$coefficients[,4]))
+ rownames(res) <- names
+ colnames(res) <- c("Coefficient","Std_error","t-value","Pr(>|t|)")
+ printCoefmat(res, P.values=TRUE, signif.stars=TRUE)
+
+ if(x$nperm > 0) {
+ cat("\nTwo-tailed tests of regression coefficients\n",'\n')
+ res2 <- as.data.frame(cbind(summary(x$reg)$coefficients[,1], x$p.param.t.2tail, x$p.perm.t.2tail))
+ rownames(res2) <- names
+ colnames(res2) <- c("Coefficient","p-param","p-perm")
+ nc <- 3
+ printCoefmat(res2, P.values=TRUE, signif.stars=TRUE, has.Pvalue = 3 && substr(colnames(res2)[3],1,6) == "p-perm")
+
+ cat("\nOne-tailed tests of regression coefficients:",'\n')
+ cat("test in the direction of the sign of the coefficient\n",'\n')
+ res1 <- as.data.frame(cbind(summary(x$reg)$coefficients[,1], x$p.param.t.1tail, x$p.perm.t.1tail))
+ rownames(res1) <- names
+ colnames(res1) <- c("Coefficient","p-param","p-perm")
+ nc <- 3
+ printCoefmat(res1, P.values=TRUE, signif.stars=TRUE, has.Pvalue = 3 && substr(colnames(res1)[3],1,6) == "p-perm")
+
+ }
+ cat("\nResidual standard error:", summary(x$reg)$sigma, "on", summary(x$reg)$df[2],"degrees of freedom",'\n')
+ cat("Multiple R-square:", summary(x$reg)$r.squared," Adjusted R-square:", summary(x$reg)$adj.r.squared,'\n')
+
+ F <- summary(x$reg)$fstatistic[[1]]
+ df1 <- summary(x$reg)$fstatistic[[2]]
+ df2 <- summary(x$reg)$fstatistic[[3]]
+ p.param.F <- pf(F, df1, df2, lower.tail=FALSE)
+ cat("\nF-statistic:", F, "on", df1, "and", df2, "DF:\n")
+ cat(" parametric p-value :", p.param.F,'\n')
+ if(x$nperm > 0) {
+ cat(" permutational p-value:", x$p.perm.F,'\n')
+ if(x$method == "raw") {
+ cat("after",x$nperm,"permutations of",x$method,"data",'\n','\n')
+ } else {
+ cat("after",x$nperm,"permutations of",x$method,"of full model",'\n','\n')
+ }
+ }
+ invisible(x)
+}
}
\details{
Two methods are available: the one by Penny and Hendy (1985), and the
- branch length score by Kuhner and Felsenstein (1994).
+ branch length score by Kuhner and Felsenstein (1994). The trees are
+ always considered as unrooted.
The topological distance is defined as twice the number of internal
branches defining different bipartitions of the tips (Penny and Hendy
--- /dev/null
+\name{lmorigin}
+\alias{lmorigin}
+\alias{print.lmorigin}
+\alias{lmorigin.ex1}
+\alias{lmorigin.ex2}
+\title{ Multiple regression through the origin }
+\description{
+Function \code{\link{lmorigin}} computes a multiple linear regression and performs tests of significance of the equation parameters (F-test of R-square and t-tests of regression coefficients) using permutations.
+
+The regression line can be forced through the origin. Testing the significance in that case requires a special permutation procedure. This option was developed for the analysis of independent contrasts, which requires regression through the origin. A permutation test, described by Legendre & Desdevises (2009), is needed to analyze contrasts that are not normally distributed.
+}
+\usage{
+lmorigin(formula, data, origin=TRUE, nperm=999, method=NULL, silent=FALSE)
+}
+
+\arguments{
+ \item{formula }{ A formula specifying the bivariate model, as in
+ \code{\link{lm}} and \code{\link{aov}}. }
+ \item{data}{ A data frame containing the two variables specified in the formula. }
+ \item{origin}{ \code{origin = TRUE} (default) to compute regression through the origin; \code{origin = FALSE} to compute multiple regression with estimation of the intercept. }
+ \item{nperm}{ Number of permutations for the tests. If \code{nperm =
+ 0}, permutation tests will not be computed. The default value is \code{nperm = 999}. For large data files, the permutation test is rather slow since the permutation procedure is not compiled. }
+ \item{method}{ \code{method = "raw"} computes t-tests of the regression coefficients by permutation of the raw data. \code{method = "residuals"} computes t-tests of the regression coefficients by permutation of the residuals of the full model. If \code{method = NULL}, permutation of the raw data is used to test the regression coefficients in regression through the origin; permutation of the residuals of the full model is used to test the regression coefficients in ordinary multiple regression. }
+ \item{silent}{ Informative messages and the time to compute the tests will not be written to the R console if silent=TRUE. Useful when the function is called by a numerical simulation function. }
+}
+
+\details{
+The permutation F-test of R-square is always done by permutation of the raw data. When there is a single explanatory variable, permutation of the raw data is used for the t-test of the single regression coefficient, whatever the method chosen by the user. The rationale is found in Anderson & Legendre (1999).
+
+The \code{print.lmorigin} function prints out the results of the parametric tests (in all cases) and the results of the permutational tests (when nperm > 0).
+}
+
+\value{
+
+ \item{reg }{The regression output object produced by function \code{lm}. }
+ \item{p.param.t.2tail }{Parametric probabilities for 2-tailed tests of the regression coefficients. }
+ \item{p.param.t.1tail }{Parametric probabilities for 1-tailed tests of the regression coefficients. Each test is carried out in the direction of the sign of the coefficient. }
+ \item{p.perm.t.2tail }{Permutational probabilities for 2-tailed tests of the regression coefficients. }
+ \item{p.perm.t.1tail }{Permutational probabilities for 1-tailed tests of the regression coefficients. Each test is carried out in the direction of the sign of the coefficient. }
+ \item{p.perm.F }{Permutational probability for the F-test of R-square. }
+ \item{origin }{TRUE is regression through the origin has been computed, FALSE if multiple regression with estimation of the intercept has been used. }
+ \item{nperm }{Number of permutations used in the permutation tests. }
+ \item{method }{Permutation method for the t-tests of the regression coefficients: \code{method = "raw"} or \code{method = "residuals"}. }
+ \item{var.names }{Vector containing the names of the variables used in the regression. }
+ \item{call }{The function call.}
+}
+
+\author{ Pierre Legendre, Universite de Montreal }
+
+\references{
+Anderson, M. J. and Legendre, P. (1999) An empirical comparison of permutation methods for tests of partial regression coefficients in a linear model. \emph{Journal of Statistical Computation and Simulation}, \bold{62}, 271--303.
+
+Legendre, P. and Desdevises, Y. (2009) Independent contrasts and regression through the origin. \emph{Journal of Theoretical Biology}, \bold{259}, 727--743.
+
+Sokal, R. R. and Rohlf, F. J. (1995) \emph{Biometry - The principles and
+ practice of statistics in biological research. Third edition.} New
+ York: W. H. Freeman.
+}
+
+\examples{
+## Example 1 from Sokal & Rohlf (1995) Table 16.1
+## SO2 air pollution in 41 cities of the USA
+data(lmorigin.ex1)
+out <- lmorigin(SO2 ~ ., data=lmorigin.ex1, origin=FALSE, nperm=99)
+out
+
+## Example 2: Contrasts computed on the phylogenetic tree of Lamellodiscus
+## parasites. Response variable: non-specificity index (NSI); explanatory
+## variable: maximum host size. Data from Table 1 of Legendre & Desdevises
+## (2009).
+data(lmorigin.ex2)
+out <- lmorigin(NSI ~ MaxHostSize, data=lmorigin.ex2, origin=TRUE, nperm=99)
+out
+
+## Example 3: random numbers
+y <- rnorm(50)
+X <- as.data.frame(matrix(rnorm(250),50,5))
+out <- lmorigin(y ~ ., data=X, origin=FALSE, nperm=99)
+out
+
+}
+
+\keyword{ multivariate }
\details{
This function implements two methods for correcting for negative values in principal coordinate analysis (PCoA). Negative eigenvalues can be produced in PCoA when decomposing distance matrices produced by coefficients that are not Euclidean (Gower and Legendre 1986, Legendre and Legendre 1998).
-Another function in this library, \code{pcoa.all}, was developed to provide principal coordinate decomposition for the function \code{PCNM} of this library. It computes the eigenvectors corresponding to all eigenvalues, positive and negative, using matrix algebra.
-
In \code{pcoa}, when negative eigenvalues are present in the decomposition results, the distance matrix D can be modified using either the Lingoes or the Cailliez procedure to produce results without negative eigenvalues.
In the Lingoes (1971) procedure, a constant c1, equal to twice absolute value of the largest negative value of the original principal coordinate analysis, is added to each original squared distance in the distance matrix, except the diagonal values. A newe principal coordinate analysis, performed on the modified distances, has at most (n-2) positive eigenvalues, at least 2 null eigenvalues, and no negative eigenvalue.
In all cases, only the eigenvectors corresponding to positive eigenvalues are shown in the output list. The eigenvectors are scaled to the square root of the corresponding eigenvalues. Gower (1966) has shown that eigenvectors scaled in that way preserve the original distance (in the D matrix) among the objects. These eigenvectors can be used to plot ordination graphs of the objects.
-We recommend not to use PCoA to produce ordinations from the chord, chi-square, abundance profile, or Hellinger distances. It is easier to first transform the community composition data using the following transformations,available in the \code{decostand} function of the \code{vegan} package, and then to carry out a principal component analysis (PCA) on the transformed data:
+We recommend not to use PCoA to produce ordinations from the chord, chi-square, abundance profile, or Hellinger distances. It is easier to first transform the community composition data using the following transformations, available in the \code{decostand} function of the \code{vegan} package, and then carry out a principal component analysis (PCA) on the transformed data:
\describe{
\item{ }{Chord transformation: decostand(spiders,"normalize") }
}
\references{
-Cailliez, F. 1983. The analytical solution of the additive constant problem. \emph{Psychometrika}, \bold{48}, 305--308.
+Cailliez, F. (1983) The analytical solution of the additive constant problem. \emph{Psychometrika}, \bold{48}, 305--308.
-Gower, J. C. 1966. Some distance properties of latent root and vector methods used in multivariate analysis. \emph{Biometrika}, \bold{53}, 325--338.
+Gower, J. C. (1966) Some distance properties of latent root and vector methods used in multivariate analysis. \emph{Biometrika}, \bold{53}, 325--338.
-Gower, J. C. and P. Legendre. 1986. Metric and Euclidean properties of
+Gower, J. C. and Legendre, P. (1986) Metric and Euclidean properties of
dissimilarity coefficients. \emph{Journal of Classification}, \bold{3}, 5--48.
-Legendre, P. and E. D. Gallagher. 2001. Ecologically meaningful transformations for ordination of species data. \emph{Oecologia}, \bold{129}, 271--280.
+Legendre, P. and Gallagher, E. D. (2001) Ecologically meaningful transformations for ordination of species data. \emph{Oecologia}, \bold{129}, 271--280.
-Legendre, P. and L. Legendre. 1998. \emph{Numerical ecology, 2nd English
+Legendre, P. and Legendre, L. (1998) \emph{Numerical ecology, 2nd English
edition.} Amsterdam: Elsevier Science BV.
-Lingoes, J. C. 1971. Some boundary conditions for a monotone analysis of symmetric matrices. \emph{Psychometrika}, \bold{36}, 195--203.
+Lingoes, J. C. (1971) Some boundary conditions for a monotone analysis of symmetric matrices. \emph{Psychometrika}, \bold{36}, 195--203.
}
\author{ Pierre Legendre, Universite de Montreal }
projects / ape.git / commitdiff