Simulate Responses with Gaussian Error or Permuted Residuals for Constrained Ordination

Function simulates a response data frame so that it adds Gaussian error to the fitted responses of Redundancy Analysis (rda), Constrained Correspondence Analysis (cca) or distance-based RDA (capscale). The function is a special case of generic simulate, and works similarly as simulate.lm.

Usage

# S3 method for class 'rda'
simulate(object, nsim = 1, seed = NULL, indx = NULL,
    rank = "full", correlated = FALSE, ...)

Arguments

object

an object representing a fitted rda, cca or capscale model.

nsim

number of response matrices to be simulated. Only one dissimilarity matrix is returned for capscale, and larger nsim is an error.

seed

an object specifying if and how the random number generator should be initialized (‘seeded’). See simulate for details.

indx

Index of residuals added to the fitted values, such as produced by shuffleSet or sample. The index can have duplicate entries so that bootstrapping is allowed. If nsim \(>1\), the output should be compliant with shuffleSet with one line for each simulation. If nsim is missing, the number of rows of indx is used to define the number of simulations, but if nsim is given, it should match number of rows in indx. If null, parametric simulation is used and Gaussian error is added to the fitted values.

rank

The rank of the constrained component: passed to predict.rda or predict.cca.

correlated

Are species regarded as correlated in parametric simulation or when indx is not given? If correlated = TRUE, multivariate Gaussian random error is generated, and if FALSE, Gaussian random error is generated separately for each species. The argument has no effect in capscale which has no information on species.

...

additional optional arguments (ignored).

Details

The implementation follows "lm" method of simulate, and adds Gaussian (Normal) error to the fitted values (fitted.rda) using function rnorm if correlated = FALSE or mvrnorm if correlated = TRUE. The standard deviations (rnorm) or covariance matrices for species (mvrnorm) are estimated from the residuals after fitting the constraints. Alternatively, the function can take a permutation index that is used to add permuted residuals (unconstrained component) to the fitted values. Raw data are used in rda. Internal Chi-square transformed data are used in cca within the function, but the returned matrix is similar to the original input data. The simulation is performed on internal metric scaling data in capscale, but the function returns the Euclidean distances calculated from the simulated data. The simulation uses only the real components, and the imaginary dimensions are ignored.

Value

If nsim = 1, returns a matrix or dissimilarities (in capscale) with similar additional arguments on random number seed as simulate. If nsim > 1, returns a similar array as returned by simulate.nullmodel with similar attributes.

Author

Jari Oksanen

See also

simulate for the generic case and for lm objects, and simulate.nullmodel for community null model simulation. Functions fitted.rda and fitted.cca return fitted values without the error component. See rnorm and mvrnorm (MASS package) for simulating Gaussian random error.

Examples

data(dune)
data(dune.env)
mod <- rda(dune ~  Moisture + Management, dune.env)
## One simulation
update(mod, simulate(mod) ~  .)
#> Call: rda(formula = simulate(mod) ~ Moisture + Management, data = dune.env)
#> 
#> -- Model Summary --
#> 
#>               Inertia Proportion Rank
#> Total         86.0878     1.0000     
#> Constrained   57.3442     0.6661    6
#> Unconstrained 28.7436     0.3339   13
#> 
#> Inertia is variance
#> 
#> -- Eigenvalues --
#> 
#> Eigenvalues for constrained axes:
#>   RDA1   RDA2   RDA3   RDA4   RDA5   RDA6 
#> 21.698 17.716  7.226  4.866  3.814  2.024 
#> 
#> Eigenvalues for unconstrained axes:
#>   PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10  PC11  PC12  PC13 
#> 7.014 4.665 4.266 2.746 2.628 2.236 1.700 1.192 0.795 0.532 0.492 0.328 0.150 
#> 
## An impression of confidence regions of site scores
plot(mod, display="sites")
for (i in 1:5) lines(procrustes(mod, update(mod, simulate(mod) ~ .)), col="blue")

## Simulate a set of null communities with permutation of residuals
simulate(mod, indx = shuffleSet(nrow(dune), 99))
#> An object of class “simulate.rda” 
#> ‘simulate index’ method (abundance, non-sequential)
#> 20 x 30 matrix
#> Number of permuted matrices = 99 
#>