3

I am trying to build a function which wraps around the ordiR2step() function from the vegan package. This function is based on the step() function.

This is the code to that works perfectly well outside of a function:

    install.packages("vegan")
    require(vegan)
    data(mite) 
    data(mite.env) #explanatory variables
    mite.hel = decostand(mite, "hel") #response variable
    mod0 <- rda(mite.hel ~ 1, mite.env)  # Model with intercept only 
    mod1 <- rda(mite.hel ~ ., mite.env)  # Model with all explanatory variables
    step.res <- ordiR2step(mod0, scope = formula(mod1), direction="forward") 
    step.res$anova

However, if I try and wrap this into a function:

    forward <- function(X, Y) {
    intercept <- rda(X ~ 1, data = Y)  # Model with intercept only
    model <- rda(X ~ ., data = Y)  # Model with all explanatory variables
    # this is where debugging is stuck
    forStep <- ordiR2step(object=intercept, scope = formula(model), 
                    direction = "forward", trace = FALSE)
    list(forStep$anova)
    }
    forward(mite.hel, mite)

I get the following error: Error in eval(expr, envir, enclos) : object 'X' not found

I know some problems can arise when trying to wrap these types of functions, and these issues have been been discussed many times on stackoverflow and other places, such as here, here and here. However, none of these solution seem to have worked properly.

From my understanding of this behaviour, the step() function, and by extension the ordiR2step(), cannot read from the environment defined within the function, but only from the work space environment, because if I define X before calling the function, everything works well. However this defeats the purpose, so I tried some proposed solutions such as:

    forward2 <- function(X, Y) {
    intercept <- do.call("rda",list(X ~ 1, data = Y))  
    model <- do.call("rda", list(X ~ ., data = Y))  
    forStep <- ordiR2step(object=intercept, scope = formula(model), 
                    direction = "forward", trace = FALSE)
    list(forStep$anova)
    }
    forward2(mite.hel, mite.env)

Same error... No dice... Any suggestions? Thanks for your time!

Community
  • 1
  • 1
Xavier GB
  • 97
  • 1
  • 7
  • 1
    Without unpacking the innards of `ordiR2step()`, I don't think this is going to be possible. The error is being raised by `ordiParseFormula()` and the call to which is buried quite deeply in the call stack at that point. What I mean is, these functions were written to work interactively and bets are off when you then start embedding these in further functions. I will take a look at the offending function and see if we can't improve it to help. The issues is actually when `ordiParseFormula()` starts to parse the formula for the next model in the forward selection chain... – Gavin Simpson May 15 '14 at 14:44
  • ...it is looking for an `X` in a few locations but not in the execution frame of your function `forward()`. – Gavin Simpson May 15 '14 at 14:45
  • 1
    The call stack is pretty deep, indeed. It seems that `ordiParseformula(..., envdepth = 6)` will get the data (the default is `envdepth = 2`). The error is triggered when `update()` tries to change the model. There seems to be no way of changing the `envdepth` value. Perhaps you have some secret agenda, but for this particular task, there is no need to wrap these three commands within `forward()`: you can issue those directly and then the analysis runs with no problem. – Jari Oksanen May 16 '14 at 13:09
  • It seems that with 2 vegan developers on the case, I am in good hands! Thanks for your comments here. Indeed, I do have a secret agenda, here I have only provided the simplest example to reproduce the error I am trying to circumvent. I am trying many different rda models that incorporate dbMEMs, and so its useful to automate the exploration and visualization of selected spatial structure. While I also tracked the error to `update()`, I loose track and can't even tell where/how `ordiParseformula()` is called... so I've reached the limit of my own debugging skills (which are quite humble). – Xavier GB May 16 '14 at 18:26
  • Anyhow, I will take your word that I can't change the `envdepth` value, so for now I will just define the object in the global environment before running my function. If you do figure this one out though, I suspect there may be other grateful explorers that may try and wrap `ordistep()` for testing some cool applications of model selection with some other functions bundled up. Cheers and thanks for the help! – Xavier GB May 16 '14 at 18:34

0 Answers0