写在前面
在生态学研究中，CWM方法已经被广泛应用几十年。当计算出CWM后，我们就可以进行进一步分析（相关性分析、方差分析、多元线性回归、线性混合效应模型等），以此来研究整个群落的性状模式，而不仅仅是单个物种的变化。然而，在性状-环境分析中，CWM方法容易出现 I 类误差膨胀，可能导致过于乐观的性状—环境关系……本文将对这些内容进行简述。
1 CWM的计算
1.1 概念
植物特征（性状）的群落加权平均值，定义为群落内植物特征（性状）的加权平均值。
从最广泛的意义上讲，功能多样性或功能结构可以定义为群落里中性状值的分布[1]。因此，性状的群落加权均值也是表征功能多样性的一种指数[2]。此外，CWM通常被理解为定义群落中的主要特征，并且与Grime的质量比假说[3]直接相关，该假说认为最丰富的物种的性状在很大程度上决定了生态系统过程。
1.2 计算公式
其计算公式如下：
其中 pi 是物种 i 对群落的贡献（通常用物种相对丰度衡量），tratii 是物种 i 的性状值。
物种相对丰度通常可以使用群落内物种的相对丰富度、相对盖度或相对生物量来代表[4]。
1.3 R语言中CWM的计算
尽管研究者使用自己编写的一段代码来计算CWM是可行的[5]，但一些R包简化了这一流程，这里我主要介绍2个最常用的R包，即FD包和BAT包。在后面的内容中，我还将提到名为weimea的包。
（1）数据准备
为了便于讲解。这里先让我们生成两个模拟数据集。
comm <- matrix(c(2,5,0,0,0,1,1,0,0,0,0,1,2,0,0,0,0,0,10,1), nrow = 4, ncol = 5, byrow = TRUE)rownames(comm) = c("Site1","Site2","Site3","Site4")colnames(comm) = c("Sp1","Sp2","Sp3","Sp4","Sp5")trait <- data.frame(Trait1 = c(1,0,0,2,0), Trait2 = c(rep("A",2), rep("B",3)))rownames(trait) = colnames(comm)
comm是群落的物种丰富度矩阵，它只能是matrix格式。
comm数据集
trait是物种的性状矩阵，包括一种定量性状和一种定性性状（分类变量）。它可以是matrix或者data.frame格式。
trait数据集
注意，性状矩阵的行名应该与物种丰度矩阵的列名完全一致。这里物种丰度未被转换为相对丰度，因为后面介绍的包默认会执行转换。
（2）BAT包计算CWM
BAT包被开发来评估群落在所有维度（分类学、系统发育和功能）上的alpha 和 beta 多样性[6]。这里我们只涉及使用其中的cwm()函数。
首先加载该包：
if(!require(BAT)) install.packages("BAT")然后使用cwm()函数即可计算CWM：
cwm(comm, trait)可以看出，定性性状被视为哑变量转换成了2个变量，它们的CWM都被计算了出来。
如果我们的数据中有NA值，该函数还支持我们通过设置参数na.rm = TRUE来处理：
cwm(comm, trait, na.rm = TRUE)（3）FD包计算CWM
FD包可以计算不同的功能多样性指标。
首先加载该包：
if(!require(FD)) install.packages("FD")然后使用dbFD()函数一次性计算多个功能多样性指标：
fd <- dbFD(x = trait, a = comm)警告信息提醒我们于群落中功能性独特物种的数量少于3个，某些功能性多样性指数无法计算，这不影响CWM的计算。
我们使用下面的代码查看CWM的计算结果：
fd$CWM这里，Trait1的CWM值和使用BAT包得到的结果一致，而Trait2的CWM值并没有被FD包计算出来。
也就是说，如果我们使用FD包计算包括CWM在内的多种功能多样性指数，我们需要提前将分类变量转换为哑变量。
2 CWM导致的争议
我之所以写这节内容，是为了让读者了解最近几年学界对CWM的质疑，以便我们对CWM有更全面的认识。
在性状—环境关系的分析中，已经有很多种方法被提出。常见的方法包括CWM分析、第四角分析和物种生态位质心（SNC）分析。这个三种方法在数学上非常相关，但微小的差异也会导致它们在统计功效上截然不同的表现。
2.1 CWM方法可能夸大了 I 型错误
2016年，Peres-Neto等人[7]的研究指出，与第四角相关性相比，使用 CWM 和 SNC 估计性状—环境相关性的精度非常低，尤其是在相关性真正为零的情况下。原文写得很详细，这里我们只是简单回顾一下。
Peres-Neto等人模拟了一个例子，其中 100 个群落中 50 个物种的丰度分布是根据单个环境梯度构建的，但单个性状的性状值是跨物种随机生成的（图1A）。另一个例子是，50 个物种的分布是根据单个性状构建的，但环境值是跨群落随机生成的（图1B）。
在这些例子中，人们应该期望性状和环境之间没有相关性，但CWM和SNC与环境因子间的相关性非常高，且表现出统计学上的显著意义。研究指出，CWM方法可能会返回过于乐观的结果， I 型错误率被夸大，错误地表明物种和样本属性之间存在联系，而实际上没有这种关系。相比之下，第四角相关性接近于零，这表明第四角方法是唯一一种能够对真实相关值产生一致估计的方法。
图1. 将性状-环境变异联系起来的三种方法。
2.2 真实案例中的CWM分析
CWM与跨群落的环境因素相关（或回归）是否真的产生高度夸大的 I 型错误？接下来，我们使用一个真实的案例来进行探索[8]。
让我们想象一片茂密的森林，林下生长着草本植物。穿过树冠的光是有限的，因此草本植物可能需要通过长高来争夺光资源。你可能想知道：林下植物的高度与穿过树冠的光量有关吗？这样的问题是群落生态学中一项常见研究任务。
图2. 森林插画（由百度AI生成）
这里我们使用 Vltava 数据集来分析一个叶片性状（比叶面积，SLA）和其中一个环境变量（林下光照强度）的群落加权平均值之间的关系。相关理论预测，耐阴物种将产生具有更高 SLA 的叶片，以优化在遮阴条件下的光合作用增益。让我们使用真实数据来查看这是否有效。
首先，让我们下载该数据集：
load (url ('https://raw.githubusercontent.com/zdealveindy/anadat-r/master/data/vltava.RData'))# 如果通过该网站手动下载则使用下面的代码
这里我们使用林冠覆盖度作为光照强度的代替，较高的覆盖度表明林下光照较少。
comm <- vltava$herbs$specover <- vltava$env$COVERE32trait <- data.frame(SLA=vltava$herbs$traits$SLA)rownames(trait) = colnames(comm)trait <- sqrt(trait) #使用平方根变换改善数据的正态分布
接着，使用BAT包的cwm()函数计算CWM：
CWM_SLA <-cwm(comm, trait, na.rm = TRUE)然后构建一元线性回归模型：
LM <- lm (CWM_SLA ~ cover)我们使用F检验来看cover对CWM_SLA的影响是否显著：
anova (LM)结果是显著的。让我们再看下两者的关系图：
plot (CWM_SLA ~ cover, xlab = 'Canopy cover [%]', ylab = 'CWM SLA')abline (lm (CWM_SLA ~ cover))
SLA与cover的关系
可见，结果表明这种关系是正向且显著的（P < 0.001）。
然而，问题在于，即使性状是随机生成的，我们也有相当大的机会获得 CWM 和环境变量之间的显著关系。
让我们重新排列物种的SLA 值，通过500次随机重排计算的CWM测试CWM和原始环境变量的回归关系（这可能会运行较长的时间）。
# 定义一个函数来执行单次模拟simulate_cwm <- function() {# 随机打乱SLA值SLA_shuffled <- sample(trait$SLA, replace = TRUE)# 创建一个临时数据框，其中包含打乱后的SLA值temp_trait <- data.frame(SLA = SLA_shuffled)rownames(temp_trait) <- rownames(trait)# 使用cwm函数计算群落加权平均SLACWM_SLA_rand <- cwm(comm, temp_trait, na.rm = TRUE)# 进行线性回归LM_rand_col <- lm(CWM_SLA_rand ~ cover)# 获取 p 值P_value <- anova(LM_rand_col)$`Pr(>F)`[1]# 返回 p 值return(P_value)}# 进行 500 次模拟P_rand <- replicate(500, simulate_cwm())# 计算小于 0.05 的 p 值的次数rejected <- sum(P_rand < 0.05)rejected
结果令人惊讶，500 个回归中大约有130多个是显著的 ，但如果测试具有正确的 I 类错误率，则它应该在 5% = 25个（总共 500 个）显著结果之内。
2.3 CWM—环境关系的max检验
解决这个问题的可能途径之一是使用“max”排列检验，该检验结合了基于行和基于列的排列的结果，并选择较高的 P 值作为参考值。
“max”检验可以通过weimea包的test_cwm()函数计算。首先让我们从github上安装它。
# install.packages("devtools")devtools::install_github("zdealveindy/weimea")
类似于BAT包，我们可以使用cwm()函数计算CWM_SLA，我们无需删除NA值，函数可以处理它们（详见函数说明文档）。
library (weimea)CWM_SLA_w <- cwm (com = comm, traits = trait)
使用test_cwm()函数可以直接测试 CWM 和环境变量的关系。你需要选择将 CWM 与环境变量关联的方法（目前仅适用于线性回归和 Pearson 相关性）。
CWM_cover <- test_cwm (cwm = CWM_SLA_w, env = cover, method = 'lm', test = 'all')CWM_cover
这个函数可以很快计算出结果。我们看到，最后的 P_max 值，这是显著的 （P = 0.02）。
我们还可以绘制出两者的关系：
plot (CWM_cover)SLA与cover的关系
最终，我们可以下结论：即使通过“max”检验，森林林下草本物种的CWM_SLA与林下光照之间的关系也显著。这有力地证明了之前的科学假设，即更强的光照可能会将物种过滤到 叶片SLA 较低的群落中。
2.4 有关CWM方法的最新讨论
上述只是校正CWM方法结果的一种途径，我们仍然有别的选择来实现我们的目的。
2018年，Zelený在一篇综述中指出[9]，传统的CWM 方法检验的假设可以分为三类，它们在原假设的表述上有所不同（图3），相应的 I 类错误率也不同。具体可见原文章，这里不再展开描述。
图3. 使用 CWM 方法测试的三类假设背后的假设差异。
粗体链接表示假设明确或隐含地假设了物种组成矩阵 （L） 与物种属性向量 （t） 或向量样本属性 （e） 之间的联系，因此此链接未经过检验。另一方面，问号表示此链接未明确或隐含确认，可以进行检验。
2023年，Lepš和de Bello提出了一种解读CWM值的新观点[10]，即如果环境梯度上的群落主要是由某个优势物种的变化导致的，或者环境梯度上的相似环境下的群落内物种的差异很小时（或者说物种的生态位重叠度很高时），那么通过CWM分析观察到的性状变化可能是准确的。
图4. 基于随机性状的CWM可能沿环境梯度产生同样显著的结果。
情景 （a） 是一个极端的例子，即沿梯度仅由一个物种（ sp3）引起 CWM 变化，而沿梯度的样地中没有其他物种的丰度变化（sp1、sp2、sp4 和 sp5）。sp3 具有的任何特征都将影响 CWM。情景 （b） 反映了物种经常同时出现的情况，这些物种具有重叠的物种生态位和具有相似环境条件、相似组成（在正文中定义为“结构化”）的样地。
此外，该研究也肯定了零模型通过随机化跨物种的性状值是可以纠正CWM带来的误差的，这种方法被称为CWM-sp方法。
因为这些文章的内容十分详实，所以本推送不可能完全概括其中的讨论，最后还是希望大家去多多阅读相关文献。
参考资料
[1] Díaz, S. & Cabid, M. (2001) Vive la différence: plant functional diversity matters to ecosystem processes. Trends in Ecology and Evolution, 16, 646–655.[2] Lavorel, S., K. Grigulis, S. McIntyre, N. S. G. Williams, D. Garden, J. Dorrough, S. Berman, F. Quétier, A. Thebault and A. Bonis (2008) Assessing functional diversity in the field - methodology matters! Functional Ecology 22:134-147.[3] Grime JP (1998) Benefits of plant diversity to ecosystems: immediate, filter and founder effects. J Ecol 86:902–910[4] Lepš J. et al. 2011. Community trait response to environment: disentangling species turnover vs intraspecific trait variability effects. Ecography 34: 856–863.[5] Community Weighted Mean Calculation. https://rpubs.com/CPEL/cwm[6] Cardoso, P., Rigal, F., & Carvalho, J. C. (2015). BAT–Biodiversity Assessment Tools, an R package for the measurement and estimation of alpha and beta taxon, phylogenetic and functional diversity. Methods in Ecology and Evolution, 6(2), 232-236.[7] Peres‐Neto, P. R., Dray, S., & ter Braak, C. J. (2017). Linking trait variation to the environment: critical issues with community‐weighted mean correlation resolved by the fourth‐corner approach. Ecography, 40(7), 806-816.[8] Community-weighted mean approach (CWM) and the fourth-corner approach. https://www.davidzeleny.net/anadat-r/doku.php/en:cwm_fc_examples[9] Zelený, D. (2018). Which results of the standard test for community‐weighted mean approach are too optimistic?. Journal of Vegetation Science, 29(6), 953-966.[10] Lep J , Bello F .Differences in trait–environment relationships: Implications for community weighted means tests[J].The Journal of Ecology, 2023, 111(11):2328-2341.DOI:10.1111/1365-2745.14172.
