### Example 1: St. Marguerite River Salmon Transect
data(salmon)
## Converting the data from data frames to to matrices:
Abundance <- log1p(as.matrix(salmon[,"Abundance",drop = FALSE]))
Environ <- as.matrix(salmon[,3L:5])
## Creating a spatial eigenvector map:
map1 <- eigenmap(x = salmon[,"Position"], weighting = wf.binary,
boundaries = c(0,20))
## Case of a single descriptor:
mca1 <- MCA(Y = Abundance, X = Environ[,"Substrate",drop = FALSE],
emobj = map1)
mca1
mca1_partest <- test.cdp(mca1)
mca1_partest
summary(mca1_partest)
par(mar = c(6,4,2,4))
plot(mca1_partest, las = 3, lwd=2)
mca1_pertest <- permute.cdp(mca1)
if (FALSE) {
## or:
mca1_pertest <- parPermute.cdp(mca1, permute = 999999)
}
mca1_pertest
summary(mca1_pertest)
plot(mca1_pertest, las = 3)
mca1_pertest$UpYXcb$C # Array containing the codependence coefficients
## With all descriptors at once:
mca2 <- MCA(Y = log1p(as.matrix(salmon[,"Abundance",drop = FALSE])),
X = as.matrix(salmon[,3L:5]), emobj = map1)
mca2
mca2_partest <- test.cdp(mca2)
mca2_partest
summary(mca2_partest)
par(mar = c(6,4,2,4))
plot(mca2_partest, las = 3, lwd=2)
mca2_pertest <- permute.cdp(mca2)
if (FALSE) {
## or:
mca2_pertest <- parPermute.cdp(mca2, permute = 999999)
}
mca2_pertest
summary(mca2_pertest)
plot(mca2_pertest, las = 3, lwd=2)
mca2_pertest$UpYXcb$C # Array containing the codependence coefficients
mca2_pertest$UpYXcb$C[,1L,] # now turned into a matrix.
### Example 2: Doubs Fish Community Transect
data(Doubs)
## Sites with no fish observed are excluded:
excl <- which(rowSums(Doubs.fish) == 0)
## Creating a spatial eigenvector map:
map2 <- eigenmap(x = Doubs.geo[-excl,"DFS"])
## The eigenvalues are in map2$lambda, the MEM eigenvectors in matrix map2$U
## MCA with multivariate response data analyzed on the basis of the Hellinger
## distance:
Y <- LGTransforms(Doubs.fish[-excl,],"Hellinger")
mca3 <- MCA(Y = Y, X=Doubs.env[-excl,], emobj = map2)
mca3_pertest <- permute.cdp(mca3)
if (FALSE) {
## or:
mca3_pertest <- parPermute.cdp(mca3, permute = 999999)
}
mca3_pertest
summary(mca3_pertest)
par(mar = c(6,4,2,4))
plot(mca3_pertest, las = 2, lwd=2)
## Array containing all the codependence coefficients:
mca3_pertest$UpYXcb$C
## Display the results along the transect
spmeans <- colMeans(Y)
pca1 <- svd(Y - rep(spmeans, each=nrow(Y)))
par(mar = c(5,5,2,5) + 0.1)
plot(y = pca1$u[,1L], x = Doubs.geo[-excl,"DFS"], pch = 21L, bg = "red",
ylab = "PCA1 loadings", xlab = "Distance from river source (km)")
## A regular transect of sites from 0 to 450 (km) spaced by 1 km intervals
## (451 sites in total). It is used for plotting spatially-explicit
## predictions.
x <- seq(0,450,1)
newdists <- matrix(NA, length(x), nrow(Doubs.geo[-excl,]))
for(i in 1L:nrow(newdists))
newdists[i,] <- abs(Doubs.geo[-excl,"DFS"] - x[i])
## Calculating predictions for the regular transect under the same set of
## environmental conditions from which the codependence model was built.
prd1 <- predict(mca3_pertest,
newdata = list(target = eigenmap.score(map2, newdists)))
## Projection of the predicted species abundance on pca1:
Uprd1 <-
(prd1 - rep(spmeans, each = nrow(prd1))) %*%
pca1$v %*% diag(pca1$d^-1)
lines(y = Uprd1[,1L], x = x, col=2, lty = 1)
## Projection of the predicted species abundance on pca2:
plot(y = pca1$u[,2L], x = Doubs.geo[-excl,"DFS"], pch = 21L, bg = "red",
ylab = "PCA2 loadings", xlab = "Distance from river source (km)")
lines(y = Uprd1[,2L], x = x, col = 2, lty = 1)
## Displaying only the observed and predicted abundance for Brown Trout.
par(new = TRUE)
plot(y = Y[,"TRU"], Doubs.geo[-excl,"DFS"], pch = 21L,
bg = "green", ylab = "", xlab = "", new = FALSE, axes = FALSE)
axis(4)
lines(y = prd1[,"TRU"], x = x, col = 3)
mtext(side = 4, "sqrt(Brown trout rel. abundance)", line = 2.5)
### Example 3: Borcard et al. Oribatid Mite
## Testing the (2-dimensional) spatial codependence between the Oribatid Mite
## community structure and environmental variables, while displaying the
## total effects of the significant variables on the community structure
## (i.e., its first principal component).
data(mite)
map3 <- eigenmap(x = mite.geo)
Y <- LGTransforms(mite.species, "Hellinger")
## Organize the environmental variables
mca4 <- MCA(Y = Y, X = mite.env, emobj = map3)
mca4_partest <- test.cdp(mca4, response.tests = FALSE)
summary(mca4_partest)
plot(mca4_partest, las = 2, lwd = 2)
plot(mca4_partest, col = rainbow(1200)[1L:1000], las = 3, lwd = 4,
main = "Codependence diagram", col.signif = "white")
## Making a regular point grid for plotting the spatially-explicit
## predictions:
rng <- list(
x = seq(min(mite.geo[,"x"]) - 0.1, max(mite.geo[,"x"]) + 0.1, 0.05),
y = seq(min(mite.geo[,"y"]) - 0.1, max(mite.geo[,"y"]) + 0.1, 0.05))
grid <- cbind(x = rep(rng[["x"]], length(rng[["y"]])),
y = rep(rng[["y"]], each = length(rng[["x"]])))
newdists <- matrix(NA, nrow(grid), nrow(mite.geo))
for(i in 1L:nrow(grid)) {
newdists[i,] <- ((mite.geo[,"x"] - grid[i,"x"])^2 +
(mite.geo[,"y"] - grid[i,"y"])^2)^0.5
}
spmeans <- colMeans(Y)
pca2 <- svd(Y - rep(spmeans, each = nrow(Y)))
prd2 <- predict(mca4_partest,
newdata = list(target = eigenmap.score(map3, newdists)))
Uprd2 <-
(prd2 - rep(spmeans, each = nrow(prd2))) %*%
pca2$v %*% diag(pca2$d^-1)
## Printing the response variable (first principal component of the mite
## community structure).
prmat <- Uprd2[,1L]
dim(prmat) <- c(length(rng$x), length(rng$y))
zlim <- c(min(min(prmat), min(pca2$u[,1L])), max(max(prmat),
max(pca2$u[,1L])))
image(z = prmat, x = rng$x, y = rng$y, asp = 1, zlim = zlim,
col = rainbow(1200L)[1L:1000], ylab = "y", xlab = "x")
points(
x=mite.geo[,"x"], y=mite.geo[,"y"], pch=21L,
bg = rainbow(1200L)[round(1+(999*(pca2$u[,1L]-zlim[1L])/(zlim[2L]-zlim[1L])),0)])
Run the code above in your browser using DataLab