R包vegan的主响应曲线（PRC）

主响应曲线（principal response curve，PRC）是重复测量设计中多元响应冗余分析的一个特例，最初被提出用于生态群落的研究中，比传统的约束排序更容易解释（van den Brink and ter Braak, 1999）。经过一系列拓展，目前已经广泛用于研究多元响应试验中A因素的效应对B因素的依赖水平，即“A+ A:B”。例如，通过跟踪试验中对照组和处理组之间的差异随时间的变化，分析处理产生的影响。

本篇简介通过R包vegan执行PRC的方法。

数据集

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”效应下，物种丰度的响应水平。

执行PRC

首先设定样本的分组，如上所述，包含了沟渠来源、采样时间以及化学处理剂量。

详情参阅：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)

同RDA，约束轴需经过检验，通过后才能说明模型是可信的。上述示例均观测了第一RDA轴的物种响应特征，如果期望使用它，则需要检验该轴。

#所选观测轴的置换检验，确定显著性，以第 1 轴为例，99 次置换
ctrl <- how(plots = Plots(strata = ditch,type = 'free'), within = Within(type = 'series'), nperm = 99)
anova(mod, permutations = ctrl, first = TRUE)

结果是显著的，表明模型可用。

参考文献

van den Brink P J, ter Braak C J F. Principal response curves: Analysis of time-dependent multivariate responses of biological community to stress. Environmental Toxicology & Chemistry, 1999, 18(2):138-148.  

友情链接