Principle Response Curves with R (II)
OK, so in the first post we saw that it is quite easy to calculate Principle Response Curves in R.
How are these interpreted?
This plot shows on the x-Axis the time and on the y-Axis the difference from the control treatments. The farther apart from the x-Axis the more different are the communities compared to the control.
We see a clearly treatment-related effect: After application at time 0, the treated communities rapidly change treatment dependent. However to the end of the experiment the treated and the control get similar again, which we may call ‘recovery’.
On the right side we see the species names. This shows the contribution of the species to this plot.
For example, ceanhora (=Caenis horaria) has the lowest weight (~-3), so this species decreases the most after application. binitent(=Bithynia tentaculata) has a positive weight and abundance increases slightly.
Let´s move back to the paper and the tables therein: Table 1 contains the proportions of variance that can be attributed to Time and Treatment. We can extract this information if we just print the prc-object:
Since a PRC is a special case of a RDA model
so Conditional refers in this output to time and Constrained to treatment and its interaction with time.
Looking at inertia we see that time explains 22% and treatment (+ the interaction with time) 33% of variance.
The second part of this table is a little harder: Either we refit the model with rda() as above and look at the summary where this information is displayed, or we dig a little in the object returned by prc().
This object is quite big with a lot of information in it, for more information see ?cca.object. The eigenvalues are stored in pyr_prc$CCA$eig and to get the explained variance per axis we need to divide these eigenvalues by their sum:
So the first axis explains 26% and the second 8.6% of variance.
As a side note:
would return the explained variances by treatment and treatment from above.
The first (displayed) axis can also be tested for ‘significance’ using a permutation test. However, observations from an experimental ditch are not independent since the same ditch was measured repeatedly during the experiment. We have to take this into account: each ditch represents a time-series. We will permute the whole series of one ditch, keeping the temporal order.
To setup such a permutation scheme we use the permute package, which is automatically loaded with vegan:
Note: I use the devel-version of the permute package here. You can install it from R-Forge via
With this setup we can create a permutation matrix. Each row therein is one permutation, the values are the row numbers of the original data set.
This can be passed to permutest, testing the first eigenvalue of our model.
which gives the significance of the first PRC-axis (cf. Table 2).
Table 2 gives also the period of significant influence of the treatment: Van den Brink and Ter Braak did this using Monte Carlo permutation tests per sampling date, using the natural log-transformed nominal dose as the explanatory variable.
Of course, this can also be done in R ;)
First we need to transform the treatment dose:
And then we run a permutation test (via anova.cca) per sampling date:
I am looping through time and for every sampling week I run a permutation test on an RDA. The results are in accordance with Table 2: there is a statistically significant effect of treatment from week 0.1 till week 19.
Besides the overall significance of treatment, they also looked which treatments differed from control in order to get a no-observed-effect Concentration (NOEC) [=the concentration below the lowest significant concentration].
Testing by permutation fails here, because there are not enough unique permutations (we have only 2 treated and 4 controls per sampling date). Therefore they applied a Williams Test  on the first principle component of a PCA on each sampling date.
This is a little bit more challenging since we need to compute a PCA, extract the scores and the run a Williams-Test on these.
Instead of using a trend-test we could use Dunnett-Contrasts - Comparing every treatment to the control. However, this approach has less power then Williams-Test.
I R we could do something like this:
So we run through time, via a for-loop. For every time-point we compute a PCA and extract the scores. Then we run the Dunnett contrasts with the multcomp package. Finally, we extract our results.
As we see, we get a different result then in the paper. Here our NOEC would be 0.9 $\mu g/L$. But I think this can be attributed to the lower power of Dunnetts test. To run a test with Williams contrasts simply change in the code above
OK, so now we are finished with Principle Response Curves. Of course, there are a lot more possibilities to analyze such data sets. My next posts will probably show some alternatives to PRC: SPEAR, mvabund and control charts.
Of course, I can´t (and don´t want) give any warranty on the correctness of these blog entries. So be critical!