The paper is:
van den Brink, P.J. & ter Braak, C.J.F. (1999). Principal response curves: Analysis of time-dependent multivariate responses of biological community to stress. Environmental Toxicology and Chemistry, 18, 138-148.
The example data used in the paper:
spdta <- structure(list(`Species 1` = c(100, 100, 100, 100, 100, 110,
110, 137.5, 157.143, 157.143, 120, 120, 150, 240, 240, 130, 130,
144.444, 216.667, 216.667), `Species 2` = c(100, 100, 100, 100,
100, 90, 90, 90, 90, 90, 80, 80, 80, 80, 80, 100, 100, 100, 100,
100), `Species 3` = c(100, 100, 100, 100, 100, 120, 120, 96,
84, 84, 140, 140, 112, 70, 70, 160, 160, 144, 96, 96), `Species 4` = c(100,
100, 100, 100, 100, 90, 90, 72, 63, 63, 80, 80, 64, 40, 40, 70,
70, 63, 42, 42), `Species 5` = c(200, 200, 200, 200, 200, 240,
240, 153.6, 117.6, 117.6, 240, 240, 153.6, 60, 60, 200, 200,
162, 72, 72), `Species 6` = c(100, 100, 100, 100, 100, 130, 130,
83.2, 63.7, 63.7, 70, 70, 44.8, 17.5, 17.5, 100, 100, 81, 36,
36)), .Names = c("Species 1", "Species 2", "Species 3", "Species 4",
"Species 5", "Species 6"), class = "data.frame", row.names = c(NA,
-20L))
# time
week <- gl(4, 5, labels=0:3)
# treatment
dose <- factor(rep(c("C","C","L","H","H"), 4), levels=c("C","L","H"),
ordered=FALSE)
Perform a principal response curve analysis:
require(vegan)
#the species data are natural log(x) transformed
mod <- prc(response = log(spdta), treatment = dose, time = week)
plot(mod)
However, the species scores are difference from the ones listed in the paper, but I can get it back by the function scores
# Species scores:
scores.sps.1999 <- scores (mod, choices=1, scaling=1,
const=-sqrt(nrow(spdta)), dis='sp')
scores.sps.1999
My question is:
The canonical coefficients are also different from the paper. I can manage to transform my results into the paper ones but don't know why.
sum_prc <- summary(mod,scaling=1)#, const=-sqrt(nrow(spdta))) sum_prc -(sum_prc$coefficients)/(sqrt(nrow(spdta))/2)
Also, the paper gives the standardized canonical coefficients
rdt
and standard deviationssdt
. How do I get this piece of information from my PRC results?rdt <- c(0,0,0,0,0.1805,0.3970,0,0.1805,0.7716,0,0.0852,0.5686) sdt <- c(1e-10,1e-10,1e-10,1e-10,0.2179,0.3,1e-10,0.2179,0.3, 1e-10,0.2179,0.3) Cdt <- 0.2*rdt/sdt