vegan包的內建資料集“pyrifos”，包含了殺蟲劑處理前後進行的11次研究中，12個溝渠中水生無脊椎動物的對數轉換後的丰度。
library(vegan)
data(pyrifos)
head(pyrifos[ ,1:10])
#行是樣本，列是物種，交叉區域為物種在樣本中的丰度值
試驗中隨機分配了12個觀測組（“溝渠隔離的微生態系統”），其中4個用作對照，其餘8個用殺蟲劑處理：劑量水平分別為0.1、0.9、6和44 mug/L，各劑量下對應2個溝渠。
從處理前的第4周到處理後的第24週期間，共進行了11次取樣，觀測水生無脊椎動物的丰度，並作了對數轉化。即整個“pyrifos”資料集包含了12×11=132個樣本的水生無脊椎動物的對數轉換後的丰度數值。
詳情參閱：https://www.rdocumentation.org/packages/vegan/versions/2.4-2/topics/pyrifos
可知，這裡所涉及的“重複測量設計”，即指在整個試驗期間對同一溝渠反覆進行了11次觀測。試驗中的主因素（A因素）對應處理劑量，次因素（B因素）對應處理時間。
然後期望透過PRC，觀測“A+ A:B”效應下，物種丰度的響應水平。
首先設定樣本的分組，如上所述，包含了溝渠來源、取樣時間以及化學處理劑量。
詳情參閱：https://www.rdocumentation.org/packages/vegan/versions/2.4-2/topics/pyrifos
#分組資訊：溝渠編號
ditch <- gl(12, 1, length = 132)
#分組資訊：取樣時間，處理前的第 4 周開始，至處理後的第 24 周
week <- gl(11, 12, labels = c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24))
#分組資訊：化學劑量濃度
dose <- factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11))
vegan中，執行PRC的函式為prc()。響應變數為物種丰度資料集，treatment對應處理，time對應時間。
#PRC 執行，詳情 ?prc
mod <- prc(response = pyrifos, treatment = dose, time = week)
如上所述，由於PRC是RDA的一種特例，因此可以首先檢視約束軸的解釋程度。在這裡，化學劑量濃度（主因素）在RDA中作為解釋變數對待，取樣時間（次因素）在RDA中作為解釋變數的同時，還作為了化學劑量濃度的協變數。對於R包vegan的rda()方法的細節部分，可參考前文。
#這裡實質上執行了一種附帶偏 RDA 的過程
mod
#對應於 RDA 函式 rda()，則上式 prc() 等同於下式 rda()
rda(pyrifos ~ dose * week + Condition(week))
既然執行了RDA，那麼會計算物種得分以及典範係數，可透過summary()函式檢視細節。注意一次只能檢視一個軸的得分。
#檢視 PRC 的物種得分及典範係數，詳情 ?prc
#需指定觀測的軸（即 RDA 的第幾約束軸），以及使用的標尺縮放型別
#這裡以 RDA 第 1 主軸為例，'symmetric' 縮放為例
prc_summary <- summary(mod, axis = 1, scaling = 'symmetric')
#提取物種得分
sp <- prc_summary$sp
#提取典範係數
coeff <- prc_summary$coefficients
繪製觀測物種響應處理和時間的PRC曲線。曲線代表了整體響應特徵，或者說處理在群落和物種水平產生的效應；圖右側顯示了物種得分。
#物種響應曲線，詳情 ?prc
#同樣地，需指定觀測的軸（即 RDA 的第幾約束軸），以及使用的標尺縮放型別
#物種太多時，不妨透過 select 選擇特定物種展示，例如這裡選擇展示一些較高丰度的物種
plot(mod, axis = 1, scaling = 'symmetric', select = colSums(pyrifos) > 100)
#所選觀測軸的置換檢驗，確定顯著性，以第 1 軸為例，99 次置換
ctrl <- how(plots = Plots(strata = ditch,type = 'free'), within = Within(type = 'series'), nperm = 99)
anova(mod, permutations = ctrl, first = TRUE)
結果是顯著的，表明模型可用。
