1
$\begingroup$

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
    
$\endgroup$
0
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.