[Vegan-commits] r640 - pkg/vegan/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 18 09:25:39 CET 2008
Author: jarioksa
Date: 2008-12-18 09:25:38 +0100 (Thu, 18 Dec 2008)
New Revision: 640
Modified:
pkg/vegan/R/kendall.post.R
Log:
kendall.post: latest version from G. Blanchet; removed useless system.time reporting
Modified: pkg/vegan/R/kendall.post.R
===================================================================
--- pkg/vegan/R/kendall.post.R 2008-12-17 18:45:13 UTC (rev 639)
+++ pkg/vegan/R/kendall.post.R 2008-12-18 08:25:38 UTC (rev 640)
@@ -9,134 +9,129 @@
mult <- match.arg(mult, c("sidak", "holm", "bonferroni"))
- a <- system.time({
-
- ##CC# Make sure Y is a matrix and find number of rows and columns of Y
- Y <- as.matrix(Y)
- n <- nrow(Y)
- p <- ncol(Y)
- if(p < 2) stop("There is a single variable in the data matrix")
+ ##CC# Make sure Y is a matrix and find number of rows and columns of Y
+ Y <- as.matrix(Y)
+ n <- nrow(Y)
+ p <- ncol(Y)
+ if(p < 2) stop("There is a single variable in the data matrix")
- ##CC# Transform the species abundances to ranks, by column
- R <- apply(Y,2,rank)
+ ##CC# Transform the species abundances to ranks, by column
+ R <- apply(Y,2,rank)
- if(missing(group)) group <- rep(1,p)
- if(length(group) != p){
- stop("The number of species in the vector differs from the total number of species")
- }
+ if(missing(group)) group <- rep(1,p)
+ if(length(group) != p){
+ stop("The number of species in the vector differs from the total number of species")
+ }
- ##CC# Separate tests for the variables in each group
- group <- as.factor(group)
- gr.lev <- levels(group)
- ngr <- nlevels(group)
+ ##CC# Separate tests for the variables in each group
+ group <- as.factor(group)
+ gr.lev <- levels(group)
+ ngr <- nlevels(group)
- gr <- as.list(1:ngr)
+ gr <- as.list(1:ngr)
- n.per.gr <- vector(length=ngr)
- for(i in 1:ngr){
- gr[[i]] <- which(group==gr.lev[i])
- n.per.gr[i] <- length(gr[[i]])
- ## Vector with the number of variables per group
- }
+ n.per.gr <- vector(length=ngr)
+ for(i in 1:ngr){
+ gr[[i]] <- which(group==gr.lev[i])
+ n.per.gr[i] <- length(gr[[i]]) # Vector with the number of
+ # variables per group
+ }
- ##===============================
- ##CC# start permutation procedure
- ##===============================
- ##CC# Initialize the counters
- counter <- as.list(1:ngr)
- for(i in 1:ngr){
- counter[[i]] <- vector(l=n.per.gr[i])
- }
- W.gr <- vector("list",ngr)
- if(ngr > 1) spear.gr <- vector("list",ngr)
+ ##===============================
+ ##CC# start permutation procedure
+ ##===============================
+ ##CC# Initialize the counters
+ counter <- as.list(1:ngr)
+ for(i in 1:ngr){
+ counter[[i]] <- vector(l=n.per.gr[i])
+ }
+ W.gr <- vector("list",ngr)
+ if(ngr > 1) spear.gr <- vector("list",ngr)
- for(i in 1:ngr){
- p.i <- n.per.gr[i]
- if(p.i < 2) stop("There is a single variable in group ",gr.lev[i])
- ##CC# Extract variables part of group i
- R.gr <- R[,gr[[i]]]
- ## Table with species of group 'i' only as columns CC#
- ##Calculate the mean of the Spearman correlations for each
- ##species in a group
- spear.mat <- cor(R.gr)
- diag(spear.mat) <- NA
- spear.mean <- apply(spear.mat, 1, mean, na.rm=TRUE)
- ##CC# Calculate Kendall's W for each variable
- W.var <- ((p.i-1)*spear.mean+1)/p.i
+ for(i in 1:ngr){
+ p.i <- n.per.gr[i]
+ if(p.i < 2) stop("There is a single variable in group ",gr.lev[i])
+ #CC# Extract variables part of
+ #group i
+ R.gr <- R[,gr[[i]]] # Table with species of group 'i' only
+ #as columns CC# Calculate the
+ #mean of the Spearman
+ #correlations for each species
+ #in a group
+ spear.mat <- cor(R.gr)
+ diag(spear.mat) <- NA
+ spear.mean <- apply(spear.mat, 1, mean, na.rm=TRUE)
+ #CC# Calculate Kendall's W for each variable
+ W.var <- ((p.i-1)*spear.mean+1)/p.i
#for(j in 1:n.per.gr[i]){
- for(j in 1:p.i){ # Test each species in turn
+ for(j in 1:p.i){ # Test each species in turn
#CC# Create a new R where the
#permuted variable has been
#removed
- R.gr.mod <- R.gr[,-j]
- #CC# Set counter
- counter[[i]][j] <- 1
- #CC# The actual permutation
- #procedure
- for(k in 1:nperm){
- var.perm <- sample(R.gr[,j])
- spear.mat.perm <- cor(cbind(R.gr.mod, var.perm))
- diag(spear.mat.perm) <- NA
- spear.mean.j.perm <- mean(spear.mat.perm[p.i,], na.rm=TRUE)
- W.perm <- ((p.i-1)*spear.mean.j.perm+1)/p.i
- if(W.perm >= W.var[j]) counter[[i]][j] <- counter[[i]][j]+1
- }
+ R.gr.mod <- R.gr[,-j]
+ ##CC# Set counter
+ counter[[i]][j] <- 1
+ ##CC# The actual permutation procedure
+ for(k in 1:nperm){
+ var.perm <- sample(R.gr[,j])
+ spear.mat.perm <- cor(cbind(R.gr.mod, var.perm))
+ diag(spear.mat.perm) <- NA
+ spear.mean.j.perm <- mean(spear.mat.perm[p.i,], na.rm=TRUE)
+ W.perm <- ((p.i-1)*spear.mean.j.perm+1)/p.i
+ if(W.perm >= W.var[j]) counter[[i]][j] <- counter[[i]][j]+1
}
- W.gr[[i]] <- W.var
- if(ngr > 1) spear.gr[[i]] <- spear.mean
}
- ##CC# Calculate P-values
- for(i in 1:ngr) {
- counter[[i]] <- counter[[i]]/(nperm+1)
- }
- ## Correction to P-values for multiple testing
- ## Write all P-values to a long vector 'vec'
- vec <- counter[[1]]
- if(ngr > 1) {
- for(i in 2:ngr) vec = c(vec, counter[[i]])
- }
- if(length(vec) != p) stop("Error in putting together vector 'vec'")
+ W.gr[[i]] <- W.var
+ if(ngr > 1) spear.gr[[i]] <- spear.mean
+ }
+ ##CC# Calculate P-values
+ for(i in 1:ngr) {
+ counter[[i]] <- counter[[i]]/(nperm+1)
+ }
+
+ ## Correction to P-values for multiple testing
+ ## Write all P-values to a long vector 'vec'
+ vec <- counter[[1]]
+ if(ngr > 1) {
+ for(i in 2:ngr) vec = c(vec, counter[[i]])
+ }
+ if(length(vec) != p) stop("Error in putting together vector 'vec'")
- if(mult == "sidak") {
- vec.corr = NA
- for(i in 1:p) vec.corr = c(vec.corr, (1-(1-vec[i])^p))
- vec.corr <- vec.corr[-1]
- }
- if(mult == "holm") vec.corr <- p.adjust(vec, method="holm")
- if(mult == "bonferroni") vec.corr <- p.adjust(vec, method="bonferroni")
+ if(mult == "sidak") {
+ vec.corr = NA
+ for(i in 1:p) vec.corr = c(vec.corr, (1-(1-vec[i])^p))
+ vec.corr <- vec.corr[-1]
+ }
+ if(mult == "holm") vec.corr <- p.adjust(vec, method="holm")
+ if(mult == "bonferroni") vec.corr <- p.adjust(vec, method="bonferroni")
- if(ngr > 1) {
- vec.gr <- vector("list",ngr)
- end <- 0
- for(i in 1:ngr){
- beg <- end+1
- end <- end + n.per.gr[i]
- vec.gr[[i]] <- vec.corr[beg:end]
- }
+ if(ngr > 1) {
+ vec.gr <- vector("list",ngr)
+ end <- 0
+ for(i in 1:ngr){
+ beg <- end+1
+ end <- end + n.per.gr[i]
+ vec.gr[[i]] <- vec.corr[beg:end]
}
- ## Create data frames containing the results
- if(ngr == 1) {
- table <- rbind(spear.mean, W.var, counter[[1]], vec.corr)
- rownames(table) <- c("Spearman.mean", "W.per.species", "Prob", "Corrected prob")
- } else {
- table <- as.list(1:ngr)
- for(i in 1:ngr) {
- table[[i]] <- rbind(spear.gr[[i]], W.gr[[i]], counter[[i]], vec.gr[[i]])
- rownames(table[[i]]) <- c("Spearman.mean", "W.per.species", "Prob", "Corrected prob")
- }
+ }
+ ## Create data frames containing the results
+ if(ngr == 1) {
+ table <- rbind(spear.mean, W.var, counter[[1]], vec.corr)
+ rownames(table) <- c("Spearman.mean", "W.per.species", "Prob", "Corrected prob")
+ } else {
+ table <- as.list(1:ngr)
+ for(i in 1:ngr) {
+ table[[i]] <- rbind(spear.gr[[i]], W.gr[[i]], counter[[i]], vec.gr[[i]])
+ rownames(table[[i]]) <- c("Spearman.mean", "W.per.species", "Prob", "Corrected prob")
}
-
- })
- a[3] <- sprintf("%2f",a[3])
- cat("Time to compute a posteriori tests (per species) =",a[3]," sec",'\n')
-
+ }
if(ngr == 1) {
out <- list(A_posteriori_tests=table, Correction.type=mult)
} else {
out <- list(A_posteriori_tests_Group=table, Correction.type=mult)
}
-
+ #
class(out) <- "kendall.post"
out
}
More information about the Vegan-commits
mailing list