理解分段式结构方程模型的构建方法

前言-写在前面
在进行一些论文研读的过程中，我发现针对模型复杂生态系统（比如农田土壤微生物的宏基因组学研究），研究者会将微生物基因丰度、环境因子等整合到一个SEM模型当中，由此来观察各变量之间直接或间接的相互影响。但是，细究某些论文后会发现，研究者的SEM计算并没有声明数据的类型，比如数据的一些基本要素，包括正态性及是否相互独立，以及是否存在时间自相关等情况；此外研究者也并未声明采用的是哪种SEM模型进行的计算。不同的SEM模型，对于数据的这些基本要素的要求不同，相同数据在不同SEM中计算得出的结论也许会有极大的差异。鉴于此，详细了解SEM的一些基本理论及适用范围，将有助于实现更加科学的模型计算，得出的结论也会更符合原始数据对应的生态学意义。
这篇文章中，我将了解到：
（1）SEM的基本概念
（2）传统SEM与分段式SEM的差别（应用条件、实现方法）
（3）有助于理论理解的2个实际例子

文中斜体为自己理解的注释。
非统计学专业翻译，译文仅供参考，欢迎批评指正。
转载须获得作者授权

原文标题：PIECEWISESEM: Piecewise structural equation modelling in R for ecology, evolution, and systematics, Jonathan S. Lefcheck, Methods in Ecology and Evolution 2015.

摘要

1，生态学家和进化生物学家依赖着逐渐复杂的统计工具来描述复杂的自然系统。其中有一个被称为结构功能模型（SEM）的工具已经在生态学领域获得了极大关注。SEM是一种路径分析工具，可用来分析解决生态系统中多种变量之间的复杂关系。
2，已有的关于SEM的计算是基于变量之间的协方差，而不是变量的真实值。虽然这种SEM可以进行各种模型的拟合，但是也限制了各种细节数据的整合。最新发展的SEM允许使用本地化估算方法同时实现正态-非正态分布、随机效应及不同相关性结构的计算，但是这种实现并不能自动化计算，这就限制了在复杂模型中SEM的计算。
3，在这里，我开发了一个基于R程序语言的开源R包PIECEWISESEM，可用来执行确证性路径分析。这个包将SEM方法扩展到了目前使用的所有线性模型（广义）、最小二乘模型（系统发生）以及混合效应模型。此外，我还提供了两个例子，第一个涉及到随机效应和时间自相关，第二个涉及到系统发育的独立性比较。
4，我的目标是提供一个用户友好且易于执行的SEM实现方法，且这种方法能够体现数据生成的生态学及方法学过程。
关键词：确证性路径分析，图形理论，混合模型，网络

引言

想要弄明白大自然的错综复杂性的欲望，是所有科学向前发展的唯一推动力。然而，直至上个世纪末，生态学家及进化生物学家也几乎只能弄清单个或少数几个因素对某个响应变量的影响。这种是因为计算能力的限制以及在严谨的实验中进行简化的必要性。然而，随着现代计算机的出现以及大尺度观测变得容易进行，使用复杂工具来处理复杂自然系统中多方面数据集的需求逐渐增加。结构方程模型（SEM）便提供了这样一种工具。

结构方程模型是一种概率模型，它在一个因果性网络中整合了多个预测因子和响应变量。模型通常使用路径图来进行表示，其中箭头代表观测变量之间的指向关系（见Fig.1和2）。这种指向关系可通过模型中对应各个路径的一系列结构方程获得。相较于传统的建模方法，SEMs具有以下两个主要特点。

1. 路径代表假设的因果关系。这其实违背了那句名言 "相关性不代表因果性"。事实上，相关性确实意味着因果性，但是因果的方向尚未确定，因为人们不知道是A导致B发生、还是B导致A发生或者A和B是第三个未知变量的结果。如果使用观察或者实验得到的先验知识，人们就能作出有依据的假设，从而判断A、B及其它可能影响它们之间关联的变量之间的因果结构。SEM允许直接检验这些假定的因果结构。从这个角度来讲，因为SEM明确地检验A导致B发生的这种假设，故其实际上与传统的线性模型相背离。关于以上这种思想的更深入的讨论超出了这篇文章的范围，但是在Pearl(2012)和Bollen&Pearl(2013)的文章中可以找到更多有关因果性及其与SEM关系的讨论。

2. 变量可同时作为预测变量和响应变量。通过让一个变量在一条路径中充当响应变量，并在另一条路径当中充当预测变量，这样SEM就可以检验并量化其它单一模型无法实现的间接或级联效应(e.g. Grace et al. 2007)。

结构方程模型中需要将生态学或进化学问题进行结构化并进行检验，且在一个单一网络中同时实现评估多个因果假设。

传统的SEMs使用最大似然方法来估计参数值，这些参数能最好地重现观察变量-协变量矩阵的整体。然后可以使用卡方检验来评估SEM的拟合优度，并于观察到的协方差矩阵进行比较。然而，这种方法假定所有的观测独立且所有变量服从（多重）正态分布。此外这种方法还限制了最小观测次数，因为这种传统SEM拟合中需要有足够的*度来估计整个变量-协变量矩阵(Grace 2006)。

这些限制催生了一种基于图论应用的有向无环的，或者说分段式的SEMs的并行发展。首先，在分段式SEM中，将各条路径图转换为（结构化的）线性方程集合，并且进行独立估算。这种将所有方程进行同时计算的全局估算转换为每个方程分开计算的局部估算的处理，会有利于各种分布（正态或非正态）及取样设计的数据的拟合 (Shipley 2000a, 2009)。此外，理论上来说分段式SEMs还可进行更小的数据集的拟合，因为拟合特定组分模型（可理解为单个局部估算方程模型）只需要对应足够的*度便可 (Shipley 2000a)。最后，它可以结合从分类学或系统发育中获得的进化距离来解决共同进化史的潜在混淆效应(Von Hardenberg & Gonzalez-Voyer 2013)。因为分段式SEM并没有合并潜在变量或复合变量，故我通常把它称为验证性路径分析。尽管有这三点不同，但鉴于Grace et al. (2012)将局部估算纳入了所谓“第三代SEM”，我还是会继续把这种分段式SEM称为SEM。

由于分段式SM并不产生有效的全局协方差矩阵，因此需要进行替代性的拟合优度检验。替代的典型方法是使用Shipley的有向分离检验。这种方法用来检验所有变量是条件独立的这一假设。说的直白一些，条件独立意味着这些未连接的变量间不存在遗失的关系（指定路径就存在关系，不指定的就不存在关系，Fig.1中的前一年夏季的海藻冠层面积与第二年的波浪扰动间就不存在关系） (Shipley 2000a)。定向分离检验的第一步，是需要获得与假设路径图有关且满足条件独立要求的最小集合，我们称之为基集。基集可以转化为一个线性方程组集合，且每一个线性方程都可以像其它任何线性模型一样被处理。对于任意一个给定的独立要求的变量的显著性，也就是p值，都可以估计且提取。可通过结合方程1所示的统计检验（Fisher’s C检验）的基集中的所有p值来进行定向分离检验。

其中
Pi表示包含k个变量的基集中的第i个独立变量。
基于这个公式，C就可以与具有2k个*度的卡方分布相比较。当有微弱的证据支持条件独立变量之和时，我们假设的关系就可看作与数据是一致的。公式中得出的C所代表的这些关系的集合有时是随机的，在这种情况下卡方检验的P值可能大于我们选择的显著性检验阈值（一般为0.05）。我们可以在Shipley(2000a, 2009)的这些文章中找到关于基集的推到过程的一些例子。

Shipley (2013)的文章中指出了Fisher’s C统计检验可以通过赤池信息准则的一个值来获得，公式如下
AIC=C+2K---公式2
其中
C指公式1中的C，K是*度的似然值（不要与基集中的独立变量数小k相混淆）。因为这个估计值并不是从最大似然数中获得，所以这个K值有时代表C统计量信息准则（CIC，参照 Cardon et al. 2011）。公式2可被扩展适用到少样本（变量）（AICc）的情况中，特别是当参数的个数超过了总样本数的n/40时，修正的AICc=C+2K(n/(n-K-1))。

分段式SEM的应用受限于基集的正确规范化及估算结果，且这种应用不能手动实现，对复杂模型则更不能实现了。为解决此问题，我提供了基于R程序语言的完全文档化且开源的程序包，并命名为PIECEWISESEM。这个包可实现分段式SEM的各种计算，包括构建基集、执行全体或各组分模型的拟合优度检验、计算AIC值、返回（标准化的）估计参数、绘制局部相关性以及产生预测信息。SEMs的构建基于一系列的结构方程，这些方程可由R中的大多数常见的线性建模函数指定，因此可以实现正态-非正态分布、层次结构以及不同的估计参数的计算。本文中，我给出了两个实际的例子。第一个结合了混合效应模型及时间相关观测，第二个涉及到由系统发育独立对比得到的非独立性数据。执行所有这些分析的数据和R代码都可以在文章附件当中获得。

例子1：风暴频率与海藻森林食物网

这个例子使用的数据来自于Byrnes et al. (2011)，他们研究了风暴事件对美国加利福利亚州海藻森林的多样性及食物网结构的影响。研究者整合了针对海藻森林的各种生态学调查数据，包括35个样点及八年的基于文献数据的潜在食物网联系、基于物理检测站的波浪高度及频率数据、基于卫星图像的海藻冠层面积数据。基于这些数据，他们使用关于这一海藻系统的先验知识以及实验数据（先验知识，也就是已经知道特定的因果性是存在的），在一个单一因果网络中汇总了这些变量。然后，他们使用传统方差-协方差SEM来估算了这个模型。

Byrnes et al. 假设冬季风暴产生的波浪扰动可能影响已有海藻的数量，这种数量会对春季冠层面积产生交互影响。春季冠层反过来会影响夏季冠层面积，夏季冠层面积又会受到物理作用力（指风暴）的影响。冠层面积的大小、春季或夏季的季节变化，都将为诸如藻类、固着底栖无脊椎动物以及它们的消费者（捕食者）等各种生物提供特殊结构的栖息地。栖息地中总物种丰度，将最终决定我们所观察到的食物网中潜在营养连接的数量（或称连接强度，指每一个观测物种对应到的捕食关系的平均数量）。

上述例子最原始的分析是使用LAVAAN包（Rosseel 2012）执行的。模型卡方拟合优度检验的结果（ka square = 8.784, P = 0.118）显示数据进行了充分的拟合。Byrnes et al. (2011)的研究表明，春季冠层面积的大小受到波浪扰动及前期冠层面积交互的强烈影响：如果前一年的冠层面积增加，波浪扰动对次年春季冠层面积的影响变得更负相关。春季冠层面积对物种丰度有着直接的影响，且通过影响夏季冠层面积间接影响着物种丰度。物种丰度又反过来增加了食物网的复杂性。然而他们注意到，路径系数的相乘结果（0.38×0.29=0.11）显示，比起间接影响，春季冠层面积对于物种丰度的直接负向影响的量级更大（标准化β=-0.23）。因此，他们认为冬季风清除的春季冠层事实上能够增加物种丰度（通过减少强烈的直接负向影响），并最终在短期内增加食物网的复杂性。然而，如果考虑到海藻面积冠层面积降低带来的影响，那珊瑚礁在受到连续多年波浪的冲击后，总物种丰度应该降低才对。

Byrnes et al. (2011)的分析中将各观测数据看作是相互独立的。然而事实上，距离较近的观测点可能具有相同的特征，且在一个观测点中，观测时间间隔较近的数据的相似程度可能高于观测间隔较远的数据。为了解决这两个问题，我使用分段式SEM重新拟合了他们的原始模型。重新分析的第一步，我将使用NLME包 (Pinheiro et al. 2013)将观测变量拟合到一般线性混合效应模型当中，由此解决样点之间的非独立性。虽然也可以使用分段式SEM计算，当我还是将变量进行log转化（而不是进行整数拟合）到泊松分布中，以有助于将与原始分析进行直接比较。对于每一个组分模型，我都拟合了样点的随机影响并且只让模型截距发生变化。然后我将组分模型整合到一个串行中，并将串行传递给sem.fit函数。这个函数可以返回SEM的定向分离检验的结果、Fisher’s C 统计量和AIC值。接下来，使用sem.coefs函数来获得标准化回归系数（通过均值和标准误进行标准化，参照Byrnes et al.）。

Fisher’s C统计量与卡方分布的数据比较（C10=15.64，P=0.11）显示，混合模型的分段式SEM生成的数据与 LAVAAN包产生的一致。重新分析的结果可见于fig1b当中。一般情况下，使用sem.model.fits函数计算的基于混合和随机效应的R^2值，分段式SEM模型解释的方差一般要比传统SEM的要大。

Fig.1a,b显示的两个模型存在着诸多不同点。首先，即使它们的影响方向不变，前一年海藻冠层面积以及其与波浪扰动间的交互作用对于春季冠层面积的影响均减少了2/3。Fig.1b中最重要的结论是春季冠层面积与物种丰度之间的负相关关系不显著。基于层次结构的嵌套观测结果，之前认为由冠层面积所导致的变异被重新分配到了随机（空间）变异当中。如此一来，基于分段式SEM的分析结果便显示了波浪扰动能直接或间接地减少春季冠层面积，从而间接地减少食物网的复杂性（波浪扰动破坏了一连串的正向干扰过程：首先减少了春季和夏季的冠层面积，进而减少了物种丰度，并最终降低了食物网的复杂性）。

第二个重分析中，我解决了样点间的非独立性及潜在的时间自相关问题。步骤为：首先保持与之前一致的随机结构；然后从NLME包(Pinheiro et al. 2013)的CAR1函数中，使用一个连续自回归的自相关结构对样点年份之间的相关性进行了建模。通过对Fisher’s C统计量与卡方分布进行比较(C8=7.84, P=0.45)，我发现这种分析可以很好地再现数据。重分析结果可见Fig.1c。

与没有自相关结构的分段式SEM相比，每个组分模型的方差解释量稍微更大一些。但是，两种分段式模型间的显著差异较少。春季冠层面积与物种丰度之间的路径仍然不显著。夏季冠层面积与波扰动之间出现了显著正相关，而且前一年的冠层面积与物种丰度间的相关性由显著变成了不显著。然而，使用AIC对两种分段式SEMs进行比较后发现，加入CAR1自相关结构的模型比只加入层次随机结构的模型的拟合效果要差 (AICc=97.69 cf. 81.44)（一般而言AIC越小，拟合效果越好）。

总的来说，这种重分析揭示了对于数据层次结构的建模，会使相同原始数据出现不同的解释结果：波浪扰动主要通过清除物种栖息地来降低食物网复杂性。虽然我的这一结论支持(Byrnes et al, 2011)的整体结论：他们的研究认为风暴事件（比如说波浪扰动）应该会降低食物网复杂性，具体过程为重复不断的波浪扰动事件降低了物种丰度，然后降低了食物网复杂性；然而我所展示出来的影响过程为，第一次风暴扰动清除了栖息地，从而降低了食物网的复杂性（也就是说，两种模型的主要结论相似，但是过程解释不相同）。此外，AIC模型比较显示，对潜在时间自相关进行建模并不能增加我们理解这个相互作用系统的能力。

Byrnes等人(2011)对模型的进一步探索发现，将总物种丰度分解为营养成分后，冠层面积的降低显著减少了藻类数量，但并不是他们一开始分析的那样减少了固着无脊椎动物或移动型消费者物种丰度（参照附件数据）。对样点的随机效应建模有可能稀释了富藻与贫藻样点间的差异，因此在简单的分段式SEM中更难观察到藻类丰度对总物种丰度的贡献 (Fig. 1b,c)。这些进一步的分析证实了对Byrnes等人 (2011)研究工作的深入探索是有必要的，因为这能使统计学输出的结论能更加符合系统生物学过程。

海绵虾的真社会性及生态成功学

在第二个例子当中，我使用来自海绵虾的一个属（Synalpheus）的种群及生态学数据去探索生态成功学（我将ecological success翻译为生态成功学，可理解为种群成功占据了优势生态位）的驱动因素。这个属的物种展示了一系列的社会结构，从成对形成到真社会性，每个群体有一个繁殖的雌性。这些Synalpheus种表现出的复杂的社会结构已经被假设为可以有生态学上的优势去获得更强的竞争能力及资源获取能力（换句话说，真社会性的物种群落结构可获得更有利的生态位，从而从各方面有利于种群的发展）。为了验证这一假设，Duffy & Macdonald(2010) 整理了Belize地区20种Synalpheus的雌性体重、寄主种类的数量(寄主范围)和单位面积丰度（proportional regional abundance）的数据。他们进一步计算了每个物种的真社会性指数。他们假设群落的真社会学物种越多（例如单个雌性生成更多的后代），群落就会占据更广泛的寄主，并由此会导致群落会更容易捍卫它们的寄主（类似于领土）（例如在研究区域内达到更高的群落物种丰度）。此外研究者还假设寄主范围可能与物种体型有关，因为大多数真社会性物种的体型都比较小。

第一步，我使用LAVAAN包的sem函数来拟合传统SEM，并且假定20各数据点（物种）是相互独立的。模型拟合的数据结果良好（ka^2=0.653，P=0.419），结果见Fig.2a。有两条显著的路径值得关注，真社会性对于基于物种体重的寄主范围强烈正向影响 (标准化β=0.58)，以及寄主范围对相对丰度的强烈正向影响(0.47)。但是，真社会性与相对丰度之间并没有显著的直接相关关系。由此可见，真社会性物种的成功（占据生态位）在极大程度上是由于它们占据有大范围的寄主。正因为这些物种占据了广泛的栖息地（寄主），因此它们组成了总群落丰度的很大一部分。但是这个模型并不支持这种假设：真社会性赋予了这些物种直接的优势以帮助物种维持并获取特殊的生存资源。

当然，Duffy & Macdonald (2010)文章中正确指出了这些物种并不独立，因为一些物种与其它物种之间相关性密切。为解决这一问题，我在Fig.1a（原文1a，事实上应该为2a）重新拟合了SEM。此外我还根据这地区Synalpheus物种 (Hultgren & Duffy 2012)的系统发育距离，进一步地固定了模型相关性矩阵。使用APE包(Paradis, Claude & Strimmer 2004)的corBrownian函数，计算得到来自系统发育树的模型相关性，并且使用NLME包(Pinheiro et al. 2013)中的gls函数拟合了组分模型。将组分模型存储在一个串行数据中，然后使用sem.fit函数计算了SEM。最终模型计算得到的数据很好（C8=0.57, P=0.751）,结果见Fig. 2b。

Fig.2所示的两个SEM有着明显的差异，其中基于系统发育的SEM出现了生物量与寄主范围的显著负相关影响（-0.32），这支持物种尺寸确实有重要影响这一结论。即使存在物种尺寸效应，真社会性对于寄主范围仍然有显著正相关影响（而且事实上相关性指数更大了，达到0.80）。和Fig.2a所示的SEM一样，真社会性对于单位区域丰度（% Abundance）并没有直接影响。同样地，这种关系也是通过寄主范围来调节的。使用CAPER包中的pgls函数重复以上分析，这个函数可以估计一个新增的标度参数λ，并产生几乎一样的结果（见附件）。

在Duffy & Macdonald (2010) 的原始文章中，他们使用多元线性回归来探索这些变量之间的关系。在考虑了体型和进化差异后，他们发现真社会性与相对丰度及寄主范围均有强烈的正相关关系。这里，使用SEM对他们的数据的进行了重分析，我发现真社会性与相对丰度之间的关系并不直接，而是占据了更广泛的寄主后的间接影响，且这个结论并不能从单个单独的多元回归中推断得出。系统发育方法在SEM中的扩充使用，有利于进化生态学中更加复杂、多变量假设数据的检验，并且能得到一些重要的观点。

讨论

在本文中，我简要介绍了分段式SEM的基本概念，并且将分段式SEM应用到了两个实际例子当中。在两个实例中，我通过引入随机变异或进化距离，确认了数据存在非独立性，并由此得出了迥异于多重回归或传统方差-协方差SEM的推论。此外我还展示了如何将一个R包-PIECEWISESEM快速且简便地应用于局部估计。事实上，这个包已被用于探索大叶藻床生态系统功能的驱动因子 (Duffy et al. 2015)，以及研究实验性河流口岸中尺度生态系统中不同营养级中功能多样性的影响因素(Lefcheck & Duffy 2015)，以及量化草原功能多样性的生物及非生物驱动力(Jing et al. 2015)。

------我是分割线------

译文总结

1，模型的适用条件上
传统SEM要求所有变量间相互独立，且均为正态分布；piecewiseSEM则对这些要求不作限制，甚至变量间存在时间相关也没问题。传统SEM对观测次数有最低要求，以使*度能达到检验的最低阈值。
2，模型的实现手段上
传统SEM对模型进行全局估计，也就是将各变量与其它变量间同时建立一个模型，并使用卡方检验计算模型拟合参数；
piecewiseSEM需要先构建一个基集，基集中的各组分方程的结合能扩展为原始变量（这类似于线性代数中的基的概念，基中的各向量是线性无关的，原矩阵的向量可由基中的部分向量组合得到）。
3，模型的评价参数上
传统SEM采用卡方检验，piecewiseSEM采用Fisher's C统计量。
4，原始数据及代码链接
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12512