Simultaneous Analysis of a Sequence of Paired Ecological Tables:
A Comparison of Several Methods.

1. Summary

A pair of ecological tables is made of one table containing environmental variables (in columns), and another table containing species data (in columns). The rows of these two tables are identical and correspond to the sites where environmental variables and species data have been measured. Such data are used to analyse the relationships between species and their environment. If sampling is repeated over time for both tables, one obtains a sequence of pairs of ecological tables. Analysing this type of data is a way to assess changes in species-environment relationships, which can be important for conservation Ecology or for global change studies. We present the COSTATIS analysis, a new data analysis method adapted to the study of this type of data, and we compare it to two other methods on the same data set, STATICO and BGCOIA. All three methods are implemented in the ade4 package for the R environment.

2. BGCOIA

This section can be used to redo the computations and graphical displays of Between-Group Co-Inertia Analysis (BGCOIA). First the ade4 library and the "meau" data set are loaded. A correlation matrix PCA (acpmil) is done on the environmental parameters table (meau$mil), and a covariance matrix PCA (acpfau) is done on the Ephemeroptera species table (meau$fau) with the dudi.pca function. Co-inertia analysis (coinfm) is then computed on these two analyses with the coinertia function. The BGCOIA (bgcoinfm) is then computed with the betweencoinertia function, using the output of co-inertia analysis (coinfm) and the site factor (meau$plan$sta).

Graphical displays are plotted with the s.arrow function, using species and environmental variable loadings (bgcoinfm$l1 and bgcoinfm$c1). The graph of sites is plotted with the s.class function, using the row scores of the two tables (bgcoinfm$msX and bgcoinfm$msY) and the site factor (meau$plan$sta).

Note that the sign of axes is arbitrary and has been changed here to get figures comparable to the figures in the AoAS paper.

library(ade4) data(meau) acpmil <- dudi.pca(meau$mil, scannf=FALSE) acpfau <- dudi.pca(meau$fau, scale=FALSE, scannf=FALSE) coinfm <- coinertia(acpmil, acpfau, scannf=FALSE) bgcoinfm <- betweencoinertia(coinfm, meau$plan$sta, scannf=FALSE) s.arrow(-bgcoinfm$l1) s.arrow(-bgcoinfm$c1) s.class(-bgcoinfm$msX, meau$plan$sta, cell=0, xlim=range(-bgcoinfm$msX[,1], -bgcoinfm$msY[,1]), ylim=range(-bgcoinfm$msX[,2], -bgcoinfm$msY[,2]), pch=21) s.class(-bgcoinfm$msY, meau$plan$sta, cell=0, xlim=range(-bgcoinfm$msX[,1], -bgcoinfm$msY[,1]), ylim=range(-bgcoinfm$msX[,2], -bgcoinfm$msY[,2]), add.plot=TRUE)

3. STATICO

This section can be used to redo the computations and graphical displays of STATICO.

3.1 Interstructure and compromise

Two Within-Group PCA are computed. The first one (wit1) is a normed Within-Group PCA on the environmental parameters table (meau$mil) with the site factor and using the withinpca function. The second one (wit2) is a Within-Group PCA on species data (meau$fau) with the site factor and using the within function.

Two K-tables are computed (kta1 and kta2), starting from the two Within-Group analyses, with the ktab.within function. The STATICO analysis (statico1) is then obtained with the statico function on these two ktabs.

Graphs are plotted using the s.corcircle (interstructure, statico1$RV.coo), s.arrow (compromise rows, statico1$l1 and compromise columns, statico1$c1), and plot (squared cosines, statico1$cos2 vs. statico1$tabw) functions.

library(ade4) data(meau) wit1 <- withinpca(meau$mil, meau$plan$dat, scan = FALSE, scal = "total") pcafau <- dudi.pca(meau$fau, scale = FALSE, scan = FALSE, nf = 2) wit2 <- within(pcafau, meau$plan$dat, scan = FALSE, nf = 2) kta1 <- ktab.within(wit1, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) kta2 <- ktab.within(wit2, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) statico1 <- statico(kta1, kta2, scan = FALSE) par(mfrow=c(2,2)) s.corcircle(statico1$RV.coo, sub = "Interstructure", possub="topleft") s.arrow(statico1$l1, sub = "Compromise", possub="topleft") s.arrow(statico1$c1, sub = "Compromise", possub="topleft") plot(statico1$tabw, statico1$cos2, xlab = "Table weights", ylab = "Cos 2", type="n") text(statico1$tabw, statico1$cos2, label=row.names(statico1$RV.coo))

3.2 Intrastructure (Variables and species)

Computations are identical to section 3.1. Graphical displays are obtained using environmental variables and species loadings (statico1$Tli and statico1$Tco) with the s.arrow function.

library(ade4) data(meau) wit1 <- withinpca(meau$mil, meau$plan$dat, scan = FALSE, scal = "total") pcafau <- dudi.pca(meau$fau, scale = FALSE, scan = FALSE, nf = 2) wit2 <- within(pcafau, meau$plan$dat, scan = FALSE, nf = 2) kta1 <- ktab.within(wit1, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) kta2 <- ktab.within(wit2, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) statico1 <- statico(kta1, kta2, scan = FALSE) layout(matrix(c(1:16), 4, 4, byrow = TRUE)) for(i in 1:4) { s.arrow(statico1$Tli[statico1$TL[,1]==i,], xlim=range(statico1$Tli[,1]), ylim=range(statico1$Tli[,2]), boxes=FALSE, clab=1.5) s.arrow(statico1$Tco[statico1$TC[,1]==i,], xlim=range(statico1$Tco[,1]), ylim=range(statico1$Tco[,2]), boxes=FALSE, clab=1.5) }

3.3 Intrastructure (Sites)

Computations are identical to section 3.1. Graphical displays are obtained using site scores (statico1$supIX and statico1$supIY) with the s.label and s.traject functions.

library(ade4) data(meau) wit1 <- withinpca(meau$mil, meau$plan$dat, scan = FALSE, scal = "total") pcafau <- dudi.pca(meau$fau, scale = FALSE, scan = FALSE, nf = 2) wit2 <- within(pcafau, meau$plan$dat, scan = FALSE, nf = 2) kta1 <- ktab.within(wit1, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) kta2 <- ktab.within(wit2, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) statico1 <- statico(kta1, kta2, scan = FALSE) sites=as.factor(rep(c(1:6),4)) par(mfrow=c(2,2)) for (i in levels(reorder(meau$plan$dat, rep(c(1, 2, 3, 4), each = 6)))) { s.traject(statico1$supIX[meau$plan$dat==i,], label=i, , xlim=c(-3.5, 2), ylim=c(-2, 3)) s.label(statico1$supIX[meau$plan$dat==i,], label=sites, , xlim=c(-3.5, 2), ylim=c(-2, 3), add.plot=T) } for (i in levels(reorder(meau$plan$dat, rep(c(1, 2, 3, 4), each = 6)))) { s.traject(statico1$supIY[meau$plan$dat==i,], label=i, , xlim=c(-3.5, 2), ylim=c(-2, 3)) s.label(statico1$supIY[meau$plan$dat==i,], label=sites, , xlim=c(-3.5, 2), ylim=c(-2, 3), add.plot=T) }

4. COSTATIS

This section can be used to redo the computations and graphical displays of COSTATIS. Two Within-Group PCA are computed. The first one (wit1) is a normed Within-Group PCA on the environmental parameters table (meau$mil) with the site factor and using the withinpca function. The second one (wit2) is a Within-Group PCA on species data (meau$fau) with the site factor and using the within function.

Two K-tables are computed (kta1 and kta2), starting from the two Within-Group analyses, with the ktab.within function. The COSTATIS analysis (costatis1) is then obtained with the costatis function on these two ktabs.

Graphical displays are obtained by superimposition of site scores (costatis1$supIX) and environmental variable loadings (costatis1$c1) or of site scores (costatis1$supIY) and species loadings (costatis1$l1), using the s.arrow and s.class functions.

data(meau) wit1 <- withinpca(meau$mil, meau$plan$dat, scan = FALSE, scal = "total") pcafau <- dudi.pca(meau$fau, scale = FALSE, scan = FALSE, nf = 2) wit2 <- within(pcafau, meau$plan$dat, scan = FALSE, nf = 2) kta1 <- ktab.within(wit1, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) kta2 <- ktab.within(wit2, colnames = rep(c("S1","S2","S3","S4","S5","S6"), 4)) costatis1 <- costatis(kta1, kta2, scan = FALSE) s.arrow(costatis1$c1*4, label=row.names(costatis1$c1), xlim=c(-2.2,1.6), ylim=c(-2,2.9), boxes=FALSE) s.class(costatis1$supIX, meau$plan$sta, cell=0, , xlim=c(-2.2,1.6), ylim=c(-2,2.9), add.plot=TRUE) s.arrow(costatis1$l1*3, label=row.names(costatis1$li), xlim=c(-2,1.5), ylim=c(-2.1,1.4), boxes=FALSE) s.class(costatis1$supIY, meau$plan$sta, cell=0, , xlim=c(-2,1.5), ylim=c(-2.1,1.4), add.plot=TRUE)