# The "trichoptera" data set contains the abundances of 56 Trichoptera species captured
# during 100 consecutive days in 22 emergence traps along a stream. The 22 traps
# (sites) form a regular transect, with geographic positions 1 to 22. The original
# daily data collected at each site were pooled into 10 survey periods for the study
# of Legendre et al. (2010) in order to reduce the very high proportion of zeros in the
# original data matrix. Order of the observations in the data set: the 22 traps (sites)
# are nested within the survey weeks, as required by the 'stimodels' and 'quicksti'
# functions.
data(trichoptera)
# log-transform the species data, excluding the Site and Date colums
trich.log <- log1p(trichoptera[,-c(1,2)])
# A log-chord transformation (Legendre & Borcard 2018) would also be appropriate for
# these data: trich.tr <- decostand(log1p(trichoptera[,-c(1,2)]), method="norm")
# Example 1. Compute the space-time interaction test using model 5. By specifying the
# number of sites (traps), the sofware assumes that they form a regular transect with
# equispaced sites. Important note to users – In real analyses, use more than 99
# permutations.
out.1 <- stimodels(trich.log, S=22, Ti=10, model="5", nperm=99)
# The interaction is significant. Because of that, test results for the main effects,
# space and time, obtained with model 5, cannot be interpreted. Tests of spatial
# variation can be done for individual times using simple RDA against dbMEM.S
# variables. Likewise, tests of temporal variation can be done for individual sites
# using simple RDA against dbMEM.T variables. A global test of the hypothesis that none
# of the times shows a significant spatial structure can be done with model 6a. For a
# global test of temporal structure at the various sites, use model 6b.
# Code not run during CRAN software tests
# Example 2. Run space-time analysis with global tests for main effects after testing
# the interaction, which is significant in this example
out.2 <- quicksti(trich.log, S=22, Ti=10, nperm=999)
# Since the interaction is significant, function 'quicksti' will carry out the
# tests of existence of a spatial (at least in one of the time periods) and temporal
# (at least at one of the sites) structures using models 6a and 6b, respectively.
# 3. Run space-time analysis for two time blocks only, i.e. time periods 6 and 7,
# then time periods 8 and 9.
# Example 3.1. Time periods 6 and 7. The interaction is not significant. In that case,
# function 'quicksti' carries out the tests of the main effects using model 5.
out.3 <- quicksti(trich.log[111:154,], S=22, Ti=2, nperm=999)
# Example 3.2. Time periods 8 and 9. The interaction is significant. In that case,
# 'quicksti' carries out the tests of the spatial effects using model 6a. Model 6b
# cannot proceed with the test of the temporal effect because Ti=2. An explanation is
# printed in the output list.
out.4 <- quicksti(trich.log[155:198,], S=22, Ti=2, nperm=999)
# 4. Illustrations of the use of 'COD.S' and 'COD.T' in STI analysis
# The following examples illustrate how to use other representations of the spatial or
# temporal relationships among observations, through arguments 'COD.S' and
# 'COD.T' of functions 'stimodels' and 'quicksti'. The trichoptera data
# are used again.
# Example 4.1. Explicitly compute dbMEMs for the spatial structure along the regular
# transect, using function 'dbmem' of adespatial, and provide it to 'stimodels'
# or 'quicksti' as argument 'COD.S'. The dbMEMs must first be computed on the
# transect, then repeated (Ti-1) times to provide Ti repeats in total.
dbMEM.S1 <- as.matrix(dbmem(1:22))
dbMEM.S10 <- dbMEM.S1
for(j in 2:10) dbMEM.S10 <- rbind(dbMEM.S10, dbMEM.S1)
out.5 <- stimodels(trich.log, S=22, Ti=10, model="5", COD.S=dbMEM.S10, nperm=999)
# Results should be identical to those in output file out.1 of Example 1, except for
# P-values which can vary slightly.
# Example 4.2. Assume now that the sampling sites have irregular positions, as
# described by the following matrix of geographic coordinates 'xy.trich'. Provide
# this matrix to argument S of 'stimodels'
xy.trich = matrix(c(1:5,11:15,21:25,31:35,41,42,rep(c(1,2),11)),22,2)
plot(xy.trich, asp=1) # Plot a quick map of the site positions
out.6 <- stimodels(trich.log, S=xy.trich, Ti=10, model="5", nperm=999)
# Example 4.3. Compute a matrix of dbMEMs for the sites. The coding matrix provided to
# argument 'COD.S' must contain repeated dbMEM.S codes because that matrix must have
# the same number of rows as matrix Y. Construct coding matrix dbMEM.S10 containing the
# dbMEM.S codes repeated 10 times.
dbMEM.S1 <- as.matrix(dbmem(xy.trich))
dbMEM.S10 = dbMEM.S1
for(i in 1:9) dbMEM.S10 <- rbind(dbMEM.S10, dbMEM.S1)
out.7 <- stimodels(trich.log, S=22, Ti=10, model="5", COD.S=dbMEM.S10, nperm=999)
# Compare the results with those obtained in the output file out6, example 4.2.
# Careful: If an analysis requires a dbMEM coding matrix for 'COD.T', the dbMEM.T
# codes must follow the required data arrangement: sites must be nested within times.
# The following function can be used to construct a dbMEM.T matrix.
MEM.T <- function(s, tt, coord.T=NULL)
# Documentation of function MEM.T –
# Generate a matrix of temporal eigenfunctions for input into stimodels,
# with sites nested within times.
# Arguments –
# s : number of space points (sites)
# tt : number of time points
# coord.T : optional matrix or vector giving the time point coordinates
{
n <- s*tt
if(is.null(coord.T)) coord.T <- as.matrix(1:tt)
MEM.TT <- as.matrix(dbmem(coord.T))
dbMEM.T <- matrix(0,n,ncol(MEM.TT)) # Empty matrix to house dbMEM.T
beg.x <- seq(1, n, by=s)
for(i in 1:tt) { # Fill tt blocks of rows with identical MEM.TT values
for(j in 1:s) dbMEM.T[(beg.x[i]+j-1),] <- MEM.TT[i,]
}
dbMEM.T
}
# Example of use of function MEM.T
dbMEM.T <- MEM.T(s=6, tt=5)
# Check the size of the dbMEM.T output matrix
dim(dbMEM.T)
# End of code not run during CRAN software tests
Run the code above in your browser using DataLab