# Replicating tables in van den Brink, P.J. & ter Braak, C.J.F. (1999) with vegan's prc

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:

1. 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)  2. Also, the paper gives the standardized canonical coefficients rdt and standard deviations sdt. 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  ## 1 Answer I'm not sure whether I follow what you tried to do. As this question is old, I'll just show the way to extract the correct species scores. require(vegan) data(pyrifos) # time week <- gl(11, 12, labels = c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24)) # treatment dose <- factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11)) ditch <- gl(12, 1, length = 132) mod <- prc(response = pyrifos, treatment = dose, time = week) #see the species scores summary(mod) #make new object containing the species scores sum_mod <- summary(mod) sum_mod$sp


These two blog posts by Eduard Szöcs are really usefull when doing PRC with Vegan in R.

PRC 1

PRC 2