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

Question

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
   

Answer 1

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