[Vegan-commits] r580 - in pkg: lmodel2/R vegan/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Nov 20 16:34:15 CET 2008
Author: jarioksa
Date: 2008-11-20 16:34:14 +0100 (Thu, 20 Nov 2008)
New Revision: 580
Modified:
pkg/lmodel2/R/lmodel2.R
pkg/vegan/R/oecosimu.R
Log:
lmodel2: slopes Inf, -Inf or NA instead of arbitrary numeric values
Modified: pkg/lmodel2/R/lmodel2.R
===================================================================
--- pkg/lmodel2/R/lmodel2.R 2008-11-20 15:26:43 UTC (rev 579)
+++ pkg/lmodel2/R/lmodel2.R 2008-11-20 15:34:14 UTC (rev 580)
@@ -315,9 +315,9 @@
x$theta," degrees",sep="",'\n')
##if((x$nperm > 0) || (x$info.slope == 1) || (x$info.CI == 1)) cat("")
if(x$info.slope == 1)
- cat("MA, SMA, RMA slopes = 555.5555 when the correlation is zero\n")
+ cat("MA, SMA, RMA slopes = NA when the correlation is zero\n")
if(x$info.CI == 1)
- cat("999.9999 or -999.9999 when the slope is infinite (90 deg. angle)\n")
+ cat("Inf or -Inf when the slope is infinite (90 deg. angle)\n")
if(x$nperm > 0) {
cat("Permutation tests of OLS, MA, RMA slopes: 1-tailed, tail corresponding to sign\n")
cat("A permutation test of r is equivalent to a permutation test of the OLS slope\n")
@@ -352,8 +352,8 @@
b0 <- ybar - b.ma*xbar
angle <- atan(b.ma)*180/pi
} else {
- b0 <- 555.5555
- b.ma <- 555.5555
+ b0 <- NA
+ b.ma <- NA
angle <- NA
}
MA.res <- c(b0, b.ma, angle)
@@ -367,8 +367,8 @@
b0 <- ybar - b.sma*xbar
angle <- atan(b.sma)*180/pi
} else {
- b0 <- 555.5555
- b.sma <- 555.5555
+ b0 <- NA
+ b.sma <- NA
angle <- NA
}
SMA.res <- c(b0, b.sma, angle)
@@ -405,13 +405,13 @@
} else {
A <- sqrt(H/(1-H))
if((b.ma*A) == -1) { # angle = -90 deg., b1inf = -infinity
- b1inf <- -999.9999
+ b1inf <- -Inf
info.CI <- 1
} else {
b1inf <- ratio * (b.ma-A) / (1+b.ma*A)
}
if((b.ma*A) == 1) { # angle = +90 deg., b1sup = +infinity
- b1sup <- 999.9999
+ b1sup <- Inf
info.CI <- 1
} else {
b1sup <- ratio * (b.ma+A) / (1-b.ma*A)
Modified: pkg/vegan/R/oecosimu.R
===================================================================
--- pkg/vegan/R/oecosimu.R 2008-11-20 15:26:43 UTC (rev 579)
+++ pkg/vegan/R/oecosimu.R 2008-11-20 15:34:14 UTC (rev 580)
@@ -1,12 +1,24 @@
`oecosimu` <-
function(comm, nestfun, method, nsimul=99,
- burnin=0, thin=1, statistic = "statistic", ...)
+ burnin=0, thin=1, statistic = "statistic",
+ control, ...)
{
nestfun <- match.fun(nestfun)
method <- match.arg(method, c("r00", "r0", "r1", "r2", "c0",
"swap", "tswap", "backtrack",
- "quasiswap"))
- comm <- ifelse(comm > 0, 1, 0)
+ "quasiswap", "permat")) # "permat" method added
+ ## eveluate according to method
+ if (method == "permat") {
+ quant <- TRUE
+ if (missing(control))
+ control <- permat.control()
+ pfull <- control$ptype == "full"
+ } else quant <- FALSE
+
+ ## conditional on quant value
+ if (!quant)
+ comm <- ifelse(comm > 0, 1, 0)
+
ind <- nestfun(comm, ...)
if (is.list(ind))
indstat <- ind[[statistic]]
@@ -14,43 +26,92 @@
indstat <- ind
n <- length(indstat)
simind <- matrix(0, nrow=n, ncol=nsimul)
- if (method %in% c("swap", "tswap")){
- checkbrd <- 1
- if (method == "tswap") {
- checkbrd <- sum(designdist(comm, "(J-A)*(J-B)", "binary"))
- M <- ncol(comm)
- N <- nrow(comm)
- checkbrd <- M*(M-1)*N*(N-1)/4/checkbrd
- thin <- round(thin*checkbrd)
+
+ ## conditional on quant value quant = FALSE
+ if (!quant) {
+ if (method %in% c("swap", "tswap")){
+ checkbrd <- 1
+ if (method == "tswap") {
+ checkbrd <- sum(designdist(comm, "(J-A)*(J-B)", "binary"))
+ M <- ncol(comm)
+ N <- nrow(comm)
+ checkbrd <- M*(M-1)*N*(N-1)/4/checkbrd
+ thin <- round(thin*checkbrd)
+ }
+ attr(simind, "thin") <- thin
+ attr(simind, "burnin") <- burnin
+ x <- comm
+ if (burnin > 0)
+ for(i in 1:burnin)
+ x <- commsimulator(x, method= method, thin = round(checkbrd))
+ for(i in 1:nsimul) {
+ x <- commsimulator(x, method = method, thin = thin)
+ tmp <- nestfun(x, ...)
+ if (is.list(tmp))
+ simind[,i] <- tmp[[statistic]]
+ else
+ simind[,i] <- tmp
+ }
}
- attr(simind, "thin") <- thin
- attr(simind, "burnin") <- burnin
- x <- comm
- if (burnin > 0)
- for(i in 1:burnin)
- x <- commsimulator(x, method= method, thin = round(checkbrd))
- for(i in 1:nsimul) {
- x <- commsimulator(x, method = method, thin = thin)
- tmp <- nestfun(x, ...)
- if (is.list(tmp))
- simind[,i] <- tmp[[statistic]]
- else
- simind[,i] <- tmp
+ else {
+ for (i in 1:nsimul) {
+ x <- commsimulator(comm, method=method)
+ tmp <- nestfun(x,...)
+ if (is.list(tmp))
+ simind[,i] <- tmp[[statistic]]
+ else
+ simind[,i] <- tmp
+ }
}
- }
- else {
- for (i in 1:nsimul) {
- x <- commsimulator(comm, method=method)
- tmp <- nestfun(x,...)
- if (is.list(tmp))
- simind[,i] <- tmp[[statistic]]
- else
- simind[,i] <- tmp
+ ## this is new addition for quantitative null model simulations
+ } else { # quant = TRUE
+ if (pfull) {
+ for (i in 1:nsimul) {
+ x <- permatfull(comm, fixedmar=control$fixedmar,
+ shuffle=control$shuffle,
+ reg=control$reg, hab=control$hab,
+ mtype=control$mtype, times=1)
+ tmp <- nestfun(x$perm[[1]])
+ if (is.list(tmp)) {
+ simind[, i] <- tmp[[statistic]]
+ } else simind[, i] <- tmp
+ }
+ attr(simind, "thin") <- NULL
+ attr(simind, "burnin") <- NULL
+ } else {
+ if (control$method %in% c("swap", "tswap")) {
+ if (burnin > 0) {
+ m <- permatswap(comm, method=control$method,
+ reg=control$reg, hab=control$hab,
+ mtype=control$mtype, times=1,
+ burnin=control$burnin, thin=0)$perm[[1]]
+ } else m <- comm
+ }
+ for (i in 1:nsimul) {
+ x <- permatswap(m, method=control$method,
+ reg=control$reg, hab=control$hab,
+ mtype=control$mtype, times=1,
+ burnin=0, thin=thin)
+ tmp <- nestfun(x$perm[[1]], ...)
+ if (is.list(tmp))
+ simind[, i] <- tmp[[statistic]]
+ else simind[, i] <- tmp
+ }
+ attr(simind, "thin") <- control$thin
+ attr(simind, "burnin") <- control$burnin
}
}
+ ## end of addition
+
z <- (indstat - rowMeans(simind))/apply(simind, 1, sd)
p <- 2*pmin(rowSums(indstat > simind), rowSums(indstat < simind))
p <- (p + 1)/(nsimul + 1)
+
+ ## ADDITION: if z is NA then it is not correct to calculate p
+ ## values try e.g. oecosimu(dune, sum, "permat")
+ if (any(is.na(z)))
+ p[is.na(z)] <- NA
+
if (is.null(names(indstat)))
names(indstat) <- statistic
if (!is.list(ind))
More information about the Vegan-commits
mailing list