I want to execute the same rda analysis sequence (fitting a model, testing significance of the model, the axis, and the term, plotting the data) on subsets of the same datasets. Therefore I wrote a function. The problem now is that the call to the anova.cca function does not work well within a function when I want to test the axis. It cannot find the Y.sub dataset
Error in eval(expr, envir, enclos) : object 'RV.sub' not found
Minimal working example:
library(vegan)
data(dune)
data(dune.env)
rda.subsetfunc <- function(RV, Y){
#RV.sub <- subset(RV, !Y$Use%in%c("BF"))
#Y.sub <- subset(Y, !Y$Use%in%c("BF"))
RV.sub <- RV; Y.sub <- Y
rda.mod <- rda(RV.sub ~ Manure, Y.sub)
axis.test <- anova(rda.mod, by = "axis")
return(list(rda.mod, axis.test))
}
rda.subsetfunc(RV = dune, Y = dune.env)
I found some other related questions like here but that seems be a lot more complicated than what I am doing. I tried to implement the do.call approach as is mentioned here but I couldn't get it to work. If it is really not possible to do this without really digging deep into functions I will find a way of programming around it. But to me, it feels like I'm trying to do something that makes total sense. So it is probably more likely that I am doing something wrong than that I am doing something impossible.
This a scoping issue in anova.cca(..., by="axis") which should find items from several different environments when it is updating the formula (I won't go into technical details). You really cannot embed the function to analyse the significances of axes. This is a known issue. We have solved that in the development version of vegan. The re-designed function in https://github.com/vegandevs/vegan seemed to work with this example. All ordination and significance functions are radically changed there, and they are not yet completely finished. We plan to release them in vegan 2.5-0 in the last quarter of 2017, but they are not finished yet.
The problem is that anova.cca(..., by = "axis") must find items that it builds within the function, and in addition it can find items that were available when the original model was built, but it cannot find items that you generate in functions that embed the function. You must circumvent this by making your embedding function to write its objects somewhere that these can be found. The easiest (but dirty) solution is to write them into parent environment with <<-. The following version of your function adds this <<- and seems to work in vegan 2.4-3
rda.subsetfunc <- function(RV, Y){
RV.sub <<- RV; Y.sub <- Y
rda.mod <- rda(RV.sub ~ Manure, Y.sub)
axis.test <- anova(rda.mod, by = "axis")
list(rda.mod, axis.test)
}
rda.subsetfunc(RV = dune, Y = dune.env)