Search code examples
rvegan

Explanation of the 'scope' part of ordiR2step

I'm trying to use the ordiR2step function from the vegan package. I'm able to get the similar ordistep function working fine with:

mrday1<-rda(y1bio~y1local + y1region + y1comp + (y1local * y1region)....)
Aiy1<-ordistep(mrday1,perm.max=200)
Aiy1$anova

But I don't think I fully understand how the 'scope' piece works, so I neither trust that ordistep is giving me what I'm looking for nor can I get ordiR2step to work (it requires scope).

In the documentation it says that the scope:

"Defines the range of models examined in the stepwise search. This should be either a single formula, or a list containing components upper and lower, both formulae."

with no example of it being used beyond

data(mite) 
data(mite.env) 
mite.hel = decostand(mite, "hel") 
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, mod1, perm.max = 200) 
step.res$anova  # Summary
#Note: this is a direct quote from the Vegan documentation

I'm confused on what the function of 'scope' is, and therefore how best to craft appropriate formulae for it. I tried:

mrday0<-yda(y1bio~1,newAbioy1)
 mrday1<-rda(y1bio~y1local + y1region + y1comp + (y1local * y1region)....)
Aiy1<-ordiR2step(mrday1,scope=mrday0, perm.max=200)
Aiy1$anova

but without fully understanding what that scope bounding function did, I can't begin to evaluate the results. Questions:

1) What does 'scope' actually do?

2) What kind of formula is it looking for?

UPDATE

The complete and functional code I used is: 

mrdayy02<-rda(y2bio ~ 1, datay2)
mrday2<-rda(y2bio~y2l + y2r 
         + y2c + y2lh + y2d
         + (y2l * y2r) + (y2l * y2c) + (y2l * y2lh) + (y2l * y2d)
         +               (y2r * y2c) + (y2r * y2lh) + (y2r * y2d)
         )
Aiy2<-ordiR2step(mrdayy02,scope=mrday2,direction="forward",R2scope= FALSE, perm.max=200)
Aiy2$anova

par(bg="transparent",new=FALSE)
plot(Aiy2,type="n",bty="n",main="RDAy2",
     xlab="RDA1",
     ylab="RDA2",
     col.main="black",col.lab="black", col.axis="white",
     xaxt="n",yaxt="n")
#abline(h=0,v=0,col="black",lwd=1)
points(Aiy2,display="species",col="gray",pch=20)
#text(rday2,display="species",col="gray")
points(Aiy2,display="cn",col="black",lwd=2)
text(Aiy2,display="cn",col="black",cex=0.5)
Solution

  • ordiR2step works only with forward selection (this is documented). It starts with the model given as the first argument. The second argument scope gives the model to which tries to advance: scope must give the formula of the largest possible model (maximum model). I think this answers your first question.

    The formula must be similar as the formula you use in your model. ordiR2step seems to be a kind function, and it will also extract the formula from the fitted ordination model.

    Your example is unreproducible (and four dots in formula would give a syntax error). However, it seems to me that myrday1 is your maximum model. So it should be used as the scope. Your myrday0 contains only the constant (~ 1) and it can be used as the model to start with. It can be used as the first argument in your function. The following should work:

    ordiR2step(myrday0, myrday1)
ordiR2step(myrday0, scope = formula(myrday1)) # same, but more explicit

    In your own example you reversed the order of these models, and gave the maximum model as the first argument (= initial model). This will not work, because ordiR2step cannot go backwards (the another alternative "both" means that after a step forward, it will try to take one step back, but it cannot take the first step from the maximum model).