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?
UPDATEThe 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)
ordiR2step
works only with forward selection (this is documented). It starts with the model given as the first argument. The second argumentscope
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 thescope
. Yourmyrday0
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: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 gobackwards
(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).