[R-sig-eco] Issue with BiodiversityR::nested.npmanova
Jari Oksanen
jari.oksanen at oulu.fi
Mon Feb 25 09:55:40 CET 2013
Yes sure, it is a version problem in the sense that the problem first appeared in BiodiversityR version 2.0. The previous BiodiversityR version (1.6) was OK with this function. However, that version was not OK with the current R: The 1.6 version of BiodiversityR triggers funny warnings, and these had to be fixed. When fixing these warnings, one line was overlooked and this triggers the error later in code.
Below is a fix for the current R against BiodiversityR 2.1
diff --git a/R/nested.npmanova.R b/R/nested.npmanova.R
index 70d1f8c..9db2142 100755
--- a/R/nested.npmanova.R
+++ b/R/nested.npmanova.R
@@ -45,7 +45,7 @@ function(formula,data,method="euc",permutations=100,warnings=FALSE){
resp <- all.vars(formula)[1]
toplev <- all.vars(formula)[2]
lowlev <- all.vars(formula)[3]
- environment(formula) <- .GlobalEnv
+ environment(formula) <- .BiodiversityR
data1 <- data
assign("data1",data,envir=.BiodiversityR)
model <- capscale(formula,data1,distance="euclidean")
Cheers, Jari Oksanen
On 23/02/2013, at 17:29 PM, Baldwin, Jim -FS wrote:
> I suspect it must be a version problem. Running Windows XP (yes, still XP) with R version 2.13.1 (2011-07-08) and just now loading vegan and BiodiversityR I get the following:
>
>> library(BiodiversityR)
> Loading required package: vegan
> Loading required package: permute
> This is vegan 2.0-2
> Warning messages:
> 1: package 'BiodiversityR' was built under R version 2.13.2
> 2: package 'vegan' was built under R version 2.13.2
> 3: package 'permute' was built under R version 2.13.2
>> library(vegan)
>> data(mite)
>>
>> mite <- mite[1:60,]
>> env <- data.frame(Fac = rep(LETTERS[1:2], each = 30), NFac = rep(letters[1:4], each = 15))
>>
>> nested.npmanova(mite ~ Fac + NFac, data = env, method = "jac", permutations = 1000)
> Total sum of squares of distance matrix: 16.77718
> Total sum of squares for non-parametric manova: 16.7771827228089
>
> Nested anova for NFac nested within Fac
>
> Df SumsofSquares F N.Perm Pr(>F)
> Fac 1 3.6153 7.7348 1000 0.347652
> NFac 2 0.9348 2.1408 1000 0.001998 **
> Residual 56 12.2270
> ---
> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> Jim
>
>
> -----Original Message-----
> From: r-sig-ecology-bounces at r-project.org [mailto:r-sig-ecology-bounces at r-project.org] On Behalf Of Kay Cichini
> Sent: Saturday, February 23, 2013 5:56 AM
> To: Jari Oksanen; r-sig-ecology at r-project.org
> Subject: Re: [R-sig-eco] Issue with BiodiversityR::nested.npmanova
>
> I also tried the following:
>
> library(BiodiversityR)
> library(vegan)
> data(mite)
>
> mite <- mite[1:60,]
> env <- data.frame(Fac = rep(LETTERS[1:2], each = 30), NFac = rep(letters[1:4], each = 15))
>
> nested.npmanova(mite ~ Fac + NFac, data = env, method = "jac", permutations = 1000)
>
> throwing this error:
>
> Total sum of squares of distance matrix: 16.77718 Fehler in model.frame(data = data1, na.action = function (object, ...) :
> Objekt 'data1' nicht gefunden
>
> I inspected the source of the function but can not figure out why data1 is not assigned correctly:
>
> function (formula, data, method = "euc", permutations = 100,
> warnings = FALSE)
> {
> randomize = function(data, toplev, lowlev) {
> newdata <- data
> orig.levs <- levels(data[, lowlev])
> nl <- length(orig.levs)
> new.levs <- orig.levs[sample(nl)]
> for (i in 1:nl) {
> subs1 <- data[, lowlev] == orig.levs[i]
> subs2 <- data[, lowlev] == new.levs[i]
> newtops <- data[subs2, toplev]
> newtops <- newtops[1]
> newtops <- rep(newtops, sum(subs1))
> newdata[subs1, toplev] <- newtops
> }
> return(newdata)
> }
> randomize2 = function(data, strata) {
> newdata <- data
> orig.levs <- levels(data[, strata])
> nl <- length(orig.levs)
> for (i in 1:nl) {
> subs <- data[, strata] == orig.levs[i]
> nsub <- sum(subs == T)
> subdata <- data[subs, ]
> newdata[subs, ] <- subdata[sample(nsub), ]
> }
> return(newdata)
> }
> ow <- options("warn")
> if (warnings == FALSE) {
> options(warn = 0)
> }
> formula <- as.formula(formula)
> .BiodiversityR <- new.env()
> environment(formula) <- .BiodiversityR
> if (length(all.vars(formula)) > 3)
> stop(paste("function only works with one main and one nested
> factor"))
> x <- eval(as.name((all.vars(formula)[1])))
> if (inherits(x, "dist")) {
> distmatrix <- as.matrix(x)
> }
> else {
> distmatrix <- as.matrix(vegdist(x, method = method))
> }
> SStot <- sum(distmatrix^2)/(2 * nrow(distmatrix))
> cat("Total sum of squares of distance matrix:", SStot, "\n")
> resp <- all.vars(formula)[1]
> toplev <- all.vars(formula)[2]
> lowlev <- all.vars(formula)[3]
> environment(formula) <- .GlobalEnv
> data1 <- data
> assign("data1", data, envir = .BiodiversityR)
> model <- capscale(formula, data1, distance = "euclidean")
> anovares <- anova(model, perm.max = 1, by = "terms")
> anovadat <- data.frame(anovares)
> df1 <- anovadat[1, 1]
> df2 <- anovadat[2, 1]
> df3 <- nrow(distmatrix) - df1 - df2 - 1
> adonis1 <- adonis(formula, data1, permutations = 2, method = method)
> adonis1 <- data.frame(adonis1$aov.tab)
> sstop <- anovadat[1, 2] <- adonis1[1, 2]
> sslow <- anovadat[2, 2] <- adonis1[2, 2]
> ssres <- anovadat[3, 2] <- adonis1[3, 2]
> vartot <- adonis1[4, 2]
> Ftop <- anovadat[1, 3] <- (sstop/df1)/(sslow/df2)
> Flow <- anovadat[2, 3] <- (sslow/df2)/(ssres/df3)
> counter <- 1
> for (i in 1:permutations) {
> data2 <- randomize(data, toplev, lowlev)
> assign("data2", data2, envir = .BiodiversityR)
> adonis2 <- adonis(formula, data = data2, method = method,
> permutations = 2)
> adonis2 <- data.frame(adonis2$aov.tab)
> Frand <- (adonis2[1, 2]/df1)/(adonis2[2, 2]/df2)
> if (Frand >= Ftop) {
> counter <- counter + 1
> }
> }
> signi <- counter/(permutations + 1)
> anovadat[1, 4] <- anovadat[2, 4] <- permutations
> anovadat[1, 5] <- signi
> counter <- 1
> for (i in 1:permutations) {
> data2 <- randomize2(data, toplev)
> assign("data2", data2, envir = .BiodiversityR)
> adonis2 <- adonis(formula, data = data2, method = method,
> permutations = 2)
> adonis2 <- data.frame(adonis2$aov.tab)
> Frand <- (adonis2[2, 2]/df2)/(adonis2[3, 2]/df3)
> if (Frand >= Flow) {
> counter <- counter + 1
> }
> }
> remove("data1", envir = .BiodiversityR)
> remove("data2", envir = .BiodiversityR)
> signi <- counter/(permutations + 1)
> anovadat[2, 5] <- signi
> colnames(anovadat) <- c("Df", "SumsofSquares", "F", "N.Perm",
> "Pr(>F)")
> mod <- paste("Nested anova for", lowlev, "nested within",
> toplev, "\n")
> head <- paste("Total sum of squares for non-parametric manova:",
> vartot, "\n")
> options(ow)
> structure(anovadat, heading = c(head, mod), Random.seed = NA,
> class = c("anova.cca", "anova", "data.frame")) }
> <environment: namespace:BiodiversityR>
>
>
> 2013/2/23 Jari Oksanen <jari.oksanen at oulu.fi>
>
>> Kay,
>>
>> This should not work if the function is correctly written. You say
>> that the terms in your formula are in data=env, but there are no
>> variables env$NFac and env$Fac in env -- but there are NFac and Fac.
>>
>> Cheers, Jari
>>
>> Sent from my iPad
>>
>> On 23.2.2013, at 12.15, "Kay Cichini" <kay.cichini at gmail.com> wrote:
>>
>>> Hi all,
>>>
>>> Does anyone have a clue why this is not working:
>>>
>>> library(BiodiversityR)
>>> library(vegan)
>>>
>>> data(mite)
>>> mite <- mite[1:60,]
>>>
>>> env <- data.frame(Fac = rep(LETTERS[1:2], each = 30), NFac =
>>> rep(letters[1:4], each = 15))
>>>
>>> nested.npmanova(mite ~ env$Fac + env$NFac, data = env, method =
>>> "jac", permutations = 1000)
>>>
>>> Regards,
>>> Kay
>>>
>>> --
>>>
>>> Kay Cichini, MSc Biol
>>>
>>> Grubenweg 22, 6071 Aldrans
>>>
>>> Tel.: 0650 9359101
>>>
>>> E-Mail: kay.cichini at gmail.com
>>>
>>> Web: www.theBioBucket.blogspot.co.at<
>> http://www.thebiobucket.blogspot.co.at/><
>> http://www.theBioBucket.blogspot.co.at>
>>> --
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-ecology mailing list
>>> R-sig-ecology at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>>
>
>
>
> --
>
> Kay Cichini, MSc Biol
>
> Grubenweg 22, 6071 Aldrans
>
> Tel.: 0650 9359101
>
> E-Mail: kay.cichini at gmail.com
>
> Web: www.theBioBucket.blogspot.co.at<http://www.thebiobucket.blogspot.co.at/><http://www.theBioBucket.blogspot.co.at>
> --
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
>
>
>
> This electronic message contains information generated by the USDA solely for the intended recipients. Any unauthorized interception of this message or the use or disclosure of the information it contains may violate the law and subject the violator to civil or criminal penalties. If you believe you have received this message in error, please notify the sender and delete the email immediately.
>
--
Jari Oksanen, Dept Biology, Univ Oulu, 90014 Finland
jari.oksanen at oulu.fi, Ph. +358 400 408593, http://cc.oulu.fi/~jarioksa
More information about the R-sig-ecology
mailing list