# ordistep: Choose a Model by Permutation Tests in Constrained Ordination In vegan: Community Ecology Package

## Description

Automatic stepwise model building for constrained ordination methods (`cca`, `rda`, `capscale`). The function `ordistep` is modelled after `step` and can do forward, backward and stepwise model selection using permutation tests. Function `ordiR2step` performs forward model choice solely on adjusted R2 and P-value, for ordination objects created by `rda` or `capscale`.

## Usage

 ```1 2 3 4 5``` ```ordistep(object, scope, direction = c("both", "backward", "forward"), Pin = 0.05, Pout = 0.1, permutations = how(nperm = 199), steps = 50, trace = TRUE, ...) ordiR2step(object, scope, Pin = 0.05, R2scope = TRUE, permutations = how(nperm = 499), trace = TRUE, R2permutations = 1000, ...) ```

## Arguments

 `object` In `ordistep`, an ordination object inheriting from `cca` or `rda`. In `ordiR2step`, the object must inherit from `rda`, that is, it must have been computed using `rda` or `capscale`. `scope` Defines the range of models examined in the stepwise search. This can be a list containing components `upper` and `lower`, both formulae. If it is a single item, it is interpeted the target scope, depending on the `direction`. If `direction` is `"forward"`, a single item is interpeted as the `upper` scope and the formula of the input `object` as the `lower` scope. See `step` for details. In `ordiR2step`, this defines the upper scope; it can also be an ordination object from with the model is extracted. `direction` The mode of stepwise search, can be one of `"both"`, `"backward"`, or `"forward"`, with a default of `"both"`. If the `scope` argument is missing, the default for `direction` is `"backward"` in `ordistep` (and `ordiR2step` does not have this argument, but only works forward). `Pin, Pout` Limits of permutation P-values for adding (`Pin`) a term to the model, or dropping (`Pout`) from the model. Term is added if P <= `Pin`, and removed if P > `Pout`. `R2scope` Use adjusted R2 as the stopping criterion: only models with lower adjusted R2 than scope are accepted. `permutations` a list of control values for the permutations as returned by the function `how`, or the number of permutations required, or a permutation matrix where each row gives the permuted indices. This is passed to `anova.cca`: see there for details. `steps` Maximum number of iteration steps of dropping and adding terms. `trace` If positive, information is printed during the model building. Larger values may give more information. `R2permutations` Number of permutations used in the estimation of adjusted R2 for `cca` using `RsquareAdj`. `...` Any additional arguments to `add1.cca` and `drop1.cca`.

## Details

The basic functions for model choice in constrained ordination are `add1.cca` and `drop1.cca`. With these functions, ordination models can be chosen with standard R function `step` which bases the term choice on AIC. AIC-like statistics for ordination are provided by functions `deviance.cca` and `extractAIC.cca` (with similar functions for `rda`). Actually, constrained ordination methods do not have AIC, and therefore the `step` may not be trusted. This function provides an alternative using permutation P-values.

Function `ordistep` defines the model, `scope` of models considered, and `direction` of the procedure similarly as `step`. The function alternates with `drop` and `add` steps and stops when the model was not changed during one step. The `-` and `+` signs in the summary table indicate which stage is performed. It is often sensible to have `Pout` > `Pin` in stepwise models to avoid cyclic adds and drops of single terms.

Function `ordiR2step` builds model forward so that it maximizes adjusted R2 (function `RsquareAdj`) at every step, and stopping when the adjusted R2 starts to decrease, or the adjusted R2 of the `scope` is exceeded, or the selected permutation P-value is exceeded (Blanchet et al. 2008). The second criterion is ignored with option `R2scope = FALSE`, and the third criterion can be ignored setting `Pin = 1` (or higher). Adjusted R2 cannot be calculated if the number of predictors is higher than the number of observations, but such models can be analysed with `R2scope = FALSE`. The R2 of `cca` is based on simulations (see `RsquareAdj`) and different runs of `ordiR2step` can give different results.

Functions `ordistep` (based on P values) and `ordiR2step` (based on adjusted R2 and hence on eigenvalues) can select variables in different order.

## Value

Functions return the selected model with one additional component, `anova`, which contains brief information of steps taken. You can suppress voluminous output during model building by setting `trace = FALSE`, and find the summary of model history in the `anova` item.

Jari Oksanen

## References

Blanchet, F. G., Legendre, P. & Borcard, D. (2008) Forward selection of explanatory variables. Ecology 89, 2623–2632.

The function handles constrained ordination methods `cca`, `rda`, `dbrda` and `capscale`. The underlying functions are `add1.cca` and `drop1.cca`, and the function is modelled after standard `step` (which also can be used directly but uses AIC for model choice, see `extractAIC.cca`). Function `ordiR2step` builds upon `RsquareAdj`.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37``` ```## See add1.cca for another example ### Dune data data(dune) data(dune.env) mod0 <- rda(dune ~ 1, dune.env) # Model with intercept only mod1 <- rda(dune ~ ., dune.env) # Model with all explanatory variables ## With scope present, the default direction is "both" ordistep(mod0, scope = formula(mod1), perm.max = 200) ## Example without scope. Default direction is "backward" ordistep(mod1, perm.max = 200) ## Example of ordistep, forward ## Not run: ordistep(mod0, scope = formula(mod1), direction="forward", perm.max = 200) ## End(Not run) ### Mite data 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 ## Example of ordiR2step with default direction = "both" ## (This never goes "backward" but evaluates included terms.) step.res <- ordiR2step(mod0, mod1, perm.max = 200) step.res\$anova # Summary table ## Example of ordiR2step with direction = "forward" ## Not run: step.res <- ordiR2step(mod0, scope = formula(mod1), direction="forward") step.res\$anova # Summary table ## End(Not run) ```

### Example output

```Loading required package: permute
This is vegan 2.4-3

Start: dune ~ 1

Df    AIC      F Pr(>F)
+ Management  3 87.082 2.8400  0.005 **
+ Moisture    3 87.707 2.5883  0.005 **
+ Manure      4 89.232 1.9539  0.010 **
+ A1          1 89.591 1.9217  0.070 .
+ Use         2 91.032 1.1741  0.280
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step: dune ~ Management

Df   AIC    F Pr(>F)
- Management  3 89.62 2.84  0.005 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Df    AIC      F Pr(>F)
+ Moisture  3 85.567 1.9764  0.005 **
+ Manure    3 87.517 1.3902  0.055 .
+ A1        1 87.424 1.2965  0.265
+ Use       2 88.284 1.0510  0.440
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step: dune ~ Management + Moisture

Df    AIC      F Pr(>F)
- Moisture    3 87.082 1.9764  0.010 **
- Management  3 87.707 2.1769  0.005 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Df    AIC      F Pr(>F)
+ Manure  3 85.762 1.1225  0.285
+ A1      1 86.220 0.8359  0.570
+ Use     2 86.842 0.8027  0.765

Call: rda(formula = dune ~ Management + Moisture, data = dune.env)

Inertia Proportion Rank
Total         84.1237     1.0000
Constrained   46.4249     0.5519    6
Unconstrained 37.6988     0.4481   13
Inertia is variance

Eigenvalues for constrained axes:
RDA1   RDA2   RDA3   RDA4   RDA5   RDA6
21.588 14.075  4.123  3.163  2.369  1.107

Eigenvalues for unconstrained axes:
PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12  PC13
8.241 7.138 5.355 4.409 3.143 2.770 1.878 1.741 0.952 0.909 0.627 0.311 0.227

Start: dune ~ A1 + Moisture + Management + Use + Manure

Df    AIC      F Pr(>F)
- A1          1 85.933 0.7933  0.660
- Use         2 86.056 0.8330  0.630
- Manure      3 87.357 1.0737  0.365
- Management  2 87.672 1.1976  0.225
- Moisture    3 88.818 1.3320  0.170

Step: dune ~ Moisture + Management + Use + Manure

Df    AIC      F Pr(>F)
- Use         2 85.762 0.8441  0.635
- Manure      3 86.842 1.1004  0.355
- Management  2 87.201 1.2054  0.315
- Moisture    3 88.546 1.4355  0.120

Step: dune ~ Moisture + Management + Manure

Df    AIC      F Pr(>F)
- Manure      3 85.567 1.1225  0.325
- Management  2 86.060 1.1986  0.255
- Moisture    3 87.517 1.5788  0.055 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step: dune ~ Moisture + Management

Df    AIC      F Pr(>F)
- Moisture    3 87.082 1.9764  0.015 *
- Management  3 87.707 2.1769  0.005 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: rda(formula = dune ~ Moisture + Management, data = dune.env)

Inertia Proportion Rank
Total         84.1237     1.0000
Constrained   46.4249     0.5519    6
Unconstrained 37.6988     0.4481   13
Inertia is variance

Eigenvalues for constrained axes:
RDA1   RDA2   RDA3   RDA4   RDA5   RDA6
21.588 14.075  4.123  3.163  2.369  1.107

Eigenvalues for unconstrained axes:
PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12  PC13
8.241 7.138 5.355 4.409 3.143 2.770 1.878 1.741 0.952 0.909 0.627 0.311 0.227

Start: dune ~ 1

Df    AIC      F Pr(>F)
+ Management  3 87.082 2.8400  0.005 **
+ Moisture    3 87.707 2.5883  0.005 **
+ Manure      4 89.232 1.9539  0.015 *
+ A1          1 89.591 1.9217  0.040 *
+ Use         2 91.032 1.1741  0.305
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step: dune ~ Management

Df    AIC      F Pr(>F)
+ Moisture  3 85.567 1.9764  0.005 **
+ Manure    3 87.517 1.3902  0.065 .
+ A1        1 87.424 1.2965  0.210
+ Use       2 88.284 1.0510  0.370
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Step: dune ~ Management + Moisture

Df    AIC      F Pr(>F)
+ Manure  3 85.762 1.1225   0.29
+ A1      1 86.220 0.8359   0.62
+ Use     2 86.842 0.8027   0.70

Call: rda(formula = dune ~ Management + Moisture, data = dune.env)

Inertia Proportion Rank
Total         84.1237     1.0000
Constrained   46.4249     0.5519    6
Unconstrained 37.6988     0.4481   13
Inertia is variance

Eigenvalues for constrained axes:
RDA1   RDA2   RDA3   RDA4   RDA5   RDA6
21.588 14.075  4.123  3.163  2.369  1.107

Eigenvalues for unconstrained axes:
PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12  PC13
8.241 7.138 5.355 4.409 3.143 2.770 1.878 1.741 0.952 0.909 0.627 0.311 0.227

Call: mite.hel ~ 1

<All variables>  0.43670383
+ WatrCont       0.26084533
+ Shrub          0.20716190
+ Topo           0.15205437
+ Substrate      0.07718348
+ SubsDens       0.02632468
<none>           0.00000000

Df     AIC     F Pr(>F)
+ WatrCont  1 -84.336 25.35  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont

<All variables>   0.4367038
+ Shrub           0.3177536
+ Topo            0.3120057
+ Substrate       0.3091579
+ SubsDens        0.3066715
<none>            0.2608453

Df     AIC     F Pr(>F)
+ Shrub  2 -88.034 3.836  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub

<All variables>   0.4367038
+ Substrate       0.3653551
+ Topo            0.3525851
+ SubsDens        0.3446967
<none>            0.3177536

Df     AIC      F Pr(>F)
+ Substrate  6 -87.768 1.8251  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub + Substrate

<All variables>   0.4367038
+ Topo            0.4004249
+ SubsDens        0.3901844
<none>            0.3653551

Df     AIC      F Pr(>F)
+ Topo  1 -90.924 4.5095  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo

<All variables>   0.4367038
+ SubsDens        0.4367038
<none>            0.4004249

Df     AIC      F Pr(>F)
+ SubsDens  1 -94.489 4.7999  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo + SubsDens

+ WatrCont      0.26085  1 -84.336 25.3499  0.002 **
+ Shrub         0.31775  2 -88.034  3.8360  0.002 **
+ Substrate     0.36536  6 -87.768  1.8251  0.002 **
+ Topo          0.40042  1 -90.924  4.5095  0.002 **
+ SubsDens      0.43670  1 -94.489  4.7999  0.002 **
<All variables> 0.43670
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Call: mite.hel ~ 1

<All variables>  0.43670383
+ WatrCont       0.26084533
+ Shrub          0.20716190
+ Topo           0.15205437
+ Substrate      0.07718348
+ SubsDens       0.02632468
<none>           0.00000000

Df     AIC     F Pr(>F)
+ WatrCont  1 -84.336 25.35  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont

<All variables>   0.4367038
+ Shrub           0.3177536
+ Topo            0.3120057
+ Substrate       0.3091579
+ SubsDens        0.3066715
<none>            0.2608453

Df     AIC     F Pr(>F)
+ Shrub  2 -88.034 3.836  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub

<All variables>   0.4367038
+ Substrate       0.3653551
+ Topo            0.3525851
+ SubsDens        0.3446967
<none>            0.3177536

Df     AIC      F Pr(>F)
+ Substrate  6 -87.768 1.8251  0.004 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub + Substrate

<All variables>   0.4367038
+ Topo            0.4004249
+ SubsDens        0.3901844
<none>            0.3653551

Df     AIC      F Pr(>F)
+ Topo  1 -90.924 4.5095  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo

<All variables>   0.4367038
+ SubsDens        0.4367038
<none>            0.4004249

Df     AIC      F Pr(>F)
+ SubsDens  1 -94.489 4.7999  0.002 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo + SubsDens