ordistep: Choose a Model by Permutation Tests in Constrained Ordination

Description Usage Arguments Details Value Author(s) References See Also Examples

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.

Author(s)

Jari Oksanen

References

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

See Also

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
Loading required package: lattice
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 

Step: R2.adj= 0 
Call: mite.hel ~ 1 
 
                R2.adjusted
<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

Step: R2.adj= 0.2608453 
Call: mite.hel ~ WatrCont 
 
                R2.adjusted
<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

Step: R2.adj= 0.3177536 
Call: mite.hel ~ WatrCont + Shrub 
 
                R2.adjusted
<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

Step: R2.adj= 0.3653551 
Call: mite.hel ~ WatrCont + Shrub + Substrate 
 
                R2.adjusted
<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

Step: R2.adj= 0.4004249 
Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo 
 
                R2.adjusted
<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

Step: R2.adj= 0.4367038 
Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo + SubsDens 
 
                 R2.adj Df     AIC       F Pr(>F)   
+ 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
Step: R2.adj= 0 
Call: mite.hel ~ 1 
 
                R2.adjusted
<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

Step: R2.adj= 0.2608453 
Call: mite.hel ~ WatrCont 
 
                R2.adjusted
<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

Step: R2.adj= 0.3177536 
Call: mite.hel ~ WatrCont + Shrub 
 
                R2.adjusted
<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

Step: R2.adj= 0.3653551 
Call: mite.hel ~ WatrCont + Shrub + Substrate 
 
                R2.adjusted
<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

Step: R2.adj= 0.4004249 
Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo 
 
                R2.adjusted
<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

Step: R2.adj= 0.4367038 
Call: mite.hel ~ WatrCont + Shrub + Substrate + Topo + SubsDens 
 
                 R2.adj Df     AIC       F Pr(>F)   
+ 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.004 **
+ 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

vegan documentation built on Dec. 2, 2017, 1:06 a.m.