微生物组分数据分析的Logistic正态多项式回归模型

摘要

人类微生物组的变化与许多人类疾病有关。新一代测序技术使量化微生物组成成为可能，而无需实验室培养。微生物组数据分析的一个重要问题是确定与不同细菌分类群相关的环境/生物协变量。微生物组学研究中的分类群计数数据通常过于分散，包括许多零。为了解释这种过度分散，我们建议使用加法logistic正态多项式回归模型将协变量与细菌成分关联起来。该模型可以自然地解释采样变量和零观测值，并允许细菌分类群之间具有灵活的协方差结构。为了选择相关的协变量并估计相应的回归系数，我们提出了一组变量选择和估计的惩罚似然估计方法。我们开发了一种蒙特卡罗期望最大化算法来实现惩罚似然估计。我们的仿真结果表明，该方法的性能优于组变量选择中的惩罚多项式logistic回归和Dirichlet多项式回归模型。我们使用将人类肠道微生物组与微量营养素联系起来的数据集来演示方法，以确定与人类肠道微生物群肠型相关的营养素。

1.简介

人体微生物组是居住在人体不同部位的微生物的集合，包括肠道和肺气道。这些微生物组的基因组构成了我们遗传和代谢景观的一个组成部分，因此有助于我们的正常生理和疾病易感性(2012年人类微生物组项目联盟)。16S核糖体RNA的靶向扩增子测序通常用于细菌生态学研究。16S rRNA基因在细菌有机体中无处不在，具有允许扩增的缓慢进化区域和允许识别的快速进化区域(Kuczynski等人，2012年)。为了确定细菌组成，研究人员首先对16S rRNA基因某些可变区域的DNA链进行PCR扩增，然后使用条形码将目标扩增子原始序列分配给样本。处理后的序列通常以独立于分类学的方式在一定相似性水平上进一步聚类为操作分类单元（OTU）(Chaffron等人，2010年;Caporaso等人，2010年)。每个OTU都有一个代表性的DNA序列，通过与已知的16S rRNA数据库（例如核糖体数据库项目）进行比较，可以将其划分为分类谱系(Cole等人，2009年)。最后，在给定的分类水平（例如，属水平），通过将同一分类单元的OTU聚合，获得分类计数以表征微生物组的组成。

人类肠道中栖息着数量最多的微生物，其数量超过基因。最近的研究表明，人类肠道微生物组与肥胖和炎症性肠病有关（参见维珍和托德，2011年查看）。因此，了解与微生物组分相关的生物/环境协变量非常重要。肠道微生物组主要由两个细菌分支组成，即拟杆菌和厚壁菌。在属水平上，Arumugam等人（2011年）结论是，肠道微生物组分的变化也大体上是分层的，他们根据拟杆菌、普雷维菌和瘤胃球菌的不同组成将人类肠道微生物组划分为三个离散的肠型。Wu等人（2011）通过进行简单的相关性分析，将饮食摄入与人类肠道微生物群联系起来。我们的研究是由宾夕法尼亚大学正在进行的微生物组研究推动的。在这项研究中，收集了98名健康志愿者的粪便样本，并基于16S rRNA测序获得了不同分类水平的肠道细菌计数数据。此外，使用食物频率问卷收集这些受试者的饮食信息，从中获得214种微量营养素的测量值。该研究的目标之一是确定与细菌组成相关的微量营养素。

本文的重点是开发统计方法来识别与肠道细菌组成相关的协变量。已经为多元计数或成分数据开发了许多回归模型。其中，最常用的是多项式logistic（ML）回归模型。然而，ML模型不允许计数数据过度分散，这在基于16S的计数数据中常见。为了允许过度分散，可以应用Dirichlet多项式（DM）回归(陈和李，2013)。DM模型的一个缺点是，它的参数数量有限，无法充分建模成分的方差和协方差。特别是，Dirichlet变量之间的依赖结构不能独立于其平均值而确定；Dirichlet变量总是负相关的，这可能并不代表微生物组数据的性质。在本文中，我们建议应用加法logistic正态多项式（LNM）回归模型(艾奇逊，1982年;Billheimer、Guttorp和Fagan，2001年)将协变量与分类计数联系起来。在LNM模型中，观察到的计数由多项式分布建模，潜在的细菌组成被视为随机变量，并由逻辑正态（LN）分布建模，其中比例被转换为遵循多元正态分布。根据建议艾奇逊（1982）LN分布允许使用比Dirichlet模型更灵活的协方差结构。然而，这种LNM模型的估计并不简单，特别是当变量选择是研究的重点时。Billheimer等人（2001）开发了一种使用马尔可夫链蒙特卡罗（MCMC）的完全贝叶斯方法，以在协变量数较少的情况下拟合这样的LNM模型。此外，Billheimer等人（2001）仅在协变量数较少且不考虑变量选择时才考虑设置。

由于我们研究的主要目标是选择与细菌组成相关的营养素，变量选择是拟合LNM模型的关键步骤。据我们所知，文献中尚未对LNM回归模型的变量选择进行研究。在本文中，我们通过引入一个组，为LNM模型开发了一个组化似然估计过程选择相关变量的惩罚函数。这种群惩罚函数已应用于其他多元回归设置(Meier van de Geer和Buhlmann，2008年;彭等，2010)。然而，由于LNM回归没有闭合形式的似然函数，标准坐标下降算法不能直接用于最大化惩罚似然函数。为了解决这个困难，我们开发了一种蒙特卡罗EM（MCEM）算法来实现惩罚似然估计，其中在期望（E）步骤中使用Metropolis–Hastings（MH）算法对未观察到的成分进行采样，在最大化（M）中使用组优化最小二乘估计'步骤选择相关变量。

文章结构如下。第2节介绍了将高维协变量与微生物组分数据联系起来的LNM回归模型。第3节提出了用于拟合模型和选择相关变量的MCEM算法。第4节给出了仿真结果，将LNM模型的变量选择性能与其他模型进行了比较。第5节将所提出的方法应用于人类肠道微生物组数据集，以确定影响决定肠道类型的细菌组成的营养物质。最后，对方法和结果进行了简要讨论第6节。

2.微生物组分数据的可加Logistic正态多项式回归模型

考虑一项人类微生物组研究，我们已经获得了关于分类群。让表示计数的随机向量细菌分类群。让是由测序深度确定的计数总数。我们请客米作为辅助统计，并对米.让是微生物分类群的基本组成，其中⁠。这意味着是中的随机向量K（K）-维单纯形⁠.调节开启⁠，我们对计数进行建模使用多项式分布，条件密度函数由

The expectation and variance of thek个第个组件是和⁠分别是。此外，让成为第页-测量协变量的维向量。在我们的数据分析中，这表示第页微量营养测量。

连接协变量的一种常见方法到作文是通过多项式（ML）回归模型。然而，微生物组研究中观察到的计数的实际变化通常高于ML模型所暗示的变化。这种过度分散来自于由于个体异质性或分类群之间的相关性而导致的潜在成分的随机变化。为了对这种过度分散的数据建模，我们考虑了一个附加的逻辑正态模型⁠此先验分布基于加法对数比率变换(⁠⁠)以映射从限制单纯形到K（K）-多维开放实空间(艾奇逊，1982年)。对数比率转换定义为

其中th分类单元被视为基团。逆算子，加法逻辑变换(⁠⁠)，哪个映射到⁠，带有k个th分量由给出

和

然后我们假设随机变量遵循多元正态分布用密度函数

哪里是平均向量是对数比率转换的成分数据的协方差矩阵。为了模拟协变量对潜在细菌组成的影响Billheimer等人（2001），我们假设

对于⁠，其中是截距和是第页-量纲回归系数，用于测量第页上的协变量k个这意味着。让是截距向量成为维系数矩阵。此外，我们表示我系数矩阵的第几行作为对于⁠。

协助解释回归系数矩阵⁠，以下是艾奇逊（1986）和Billheimer等人（2001），我们定义了一个扰动算子用于组成和另一个变量在里面通过

此运算符表示原始组成被“打扰”了形成新的构图⁠.何时⁠，的k个th分类单元增加了其相对丰度。根据转换后的数据，我们有也很容易验证哪里一是实际空间中的任意数字，并且是单纯形的中心⁠，表示相同的组成。

基于这个扰动算子，我们得到

哪里代表我第行系数矩阵⁠.矢量表示单纯形中的组件因此，在没有来自协变量的干扰的情况下，它被解释为“基线组成”。向量通过协变量中的单位变化测量成分相对于基线的变化我当其他协变量不变时。此外，当⁠，协变我与k个th分类单元和其他分类单元之间的关联为负。协变量中一个单位变化对细菌组成的干扰程度我是通过以下方式测量的

哪里是单位矩阵和是1s的矢量。

3.通过惩罚似然估计进行变量选择和参数估计

假设我们有n个随机向量的副本⁠，表示为对于⁠.让是未观察到的我第个样品和是其log-ratio变换，如(1)。我们考虑当协变量的数量第页很大，我们分析的主要目标是选择对细菌组成有影响的变量。即使在第页很小，基于似然的LNM模型推断并不简单，因为没有闭合形式的对数似然函数可用。

我们建议开发一种惩罚似然方法来选择与细菌组成相关的协变量。请注意，如果我这个协变量与细菌组成无关，我们有因此⁠这促使了对参数的估计=(⁠⁠)基于最大化以下惩罚对数似然函数，

哪里是的对数似然函数我基于LNM模型的第次观测，以及是调谐参数。小组惩罚函数(袁和林，2006)组合所有分类群的信息并诱导系数矩阵的行稀疏性⁠.由于log-likelihood函数没有闭合形式的表达式⁠，我们无法执行优化(2)直接。

3.1. LNM模型估计的MCEM算法

对于给定的调谐参数⁠，完整的数据组受惩罚的log-likelihood可以写成

要实现EM算法，请在t吨第E步，我们需要计算期望的完整数据惩罚对数似然(3),

其中期望是关于⁠、和是第M步。

对于我第th个样本，的条件分布给定观测数据和当前参数值为

哪里⁠为了计算条件期望，我们使用MH算法从该条件分布中进行采样。具体来说，对于第页第MH步，我们提出了一个新的向量从多元正态分布，⁠，然后计算Metropolis接受率

然后我们画一个随机数u个从均匀分布并接受建议的新值如果否则保留之前的值。在500次初始烧入后，我们选择下一次R（右）MH样本计算E步骤中的条件期望。

对于M步骤，很容易检查和表示为

对于⁠，其中⁠。

鉴于和⁠，我们下次更新通过最小化以下完整数据对数似然的负惩罚期望值，

哪里估计依据为⁠.自是一个正定矩阵，这是一个凸优化问题，具有唯一的解。我们可以使用块坐标下降算法⁠具体来说，我们写并定义⁠,和具有对于⁠.自是下对角线矩阵，是零向量当且仅当是零向量。然后是目标函数可以写为

最小化(6)然后将其简化为多元回归的标准群拉索问题，并在Peng等人（2010）可以直接应用。

由于在EM算法的E步中使用了MCMC，因此我们不期望参数值精确收敛。相反，我们通过检查参数估计图来检查以确保估计值稳定。对于我们的模拟和实际数据分析，该算法稳定得很快，通常在10 EM步内。由于E步骤中采样的不确定性，在估计值稳定后，我们运行额外的S公司EM逐步计算，并将估计值的中位数作为参数的最终估计值。我们使用用于仿真和用于实际数据分析。

3.2. 调整参数选择

我们使用五倍交叉验证来选择调谐参数⁠对于每个⁠，我们将样本划分为五个不重叠的子集。表示子集五通过以及通过以下方式收集剩余的个体子集⁠.基于⁠，我们执行MCEM算法以获得⁠，然后计算数据集的相对预测误差通过

哪里是折叠中的样本数五,⁠、和是针对我褶皱中的第个样品，带有和对应的总计数。交叉验证得分由⁠.调谐参数选择得分最小的。

4.仿真评估

4.1. 仿真模型和方法比较

我们用分类群计数数据模拟数据集和协变量数据⁠.给定⁠，我们生成我-第个个人

使用组件设置为运算符中的基本组件⁠.协变量向量是从第页-均值为零的多维正态分布和多项式衰减协方差矩阵由提供⁠。为了模拟稀疏协变量效应，我们选择相关协变量的数量第页与的协变量控制模型稀疏性。系数矩阵的非零元素由区间内的均匀分布生成⁠.误差向量模拟自具有⁠，其中用于控制平均信噪比（SNR）。信噪比的水平决定了模型中的过度分散程度，其中较小的信噪比表示协变量对成分的总体影响较小，因此在模型中产生了较高的过度分散水平。相反，较大的信噪比表示较低的过频散水平。在极端情况下，除了协变量的影响外，没有额外的变化，即固定为零，模拟数据集时不会出现过分散走向无限。最后，我们生成分类群阅读深度在区间内的均匀分布⁠，并模拟个体的观测计数向量我通过⁠。

我们将LNM模型的变量选择性能与三个具有相似组的替代模型进行了比较处罚：

• 1

  logistic normal（LN）模型将对数比率转换视为多元响应，然后应用组变量选择程序Peng等人（2010）选择变量。

• 2

  带群的Dirichlet多项式（DM）模型惩罚函数由提出陈和李（2013）。

• 三。

  直接用组对计数数据建模的多项式logistic（ML）模型惩罚多项式logistic回归模型。

在LN模型中，将模拟数据集中的零计数替换为伪计数值0.05，以便于估计。参考组件在构建对数比率的组件。五倍交叉验证用于选择所有四个模型中的调谐参数。

我们考虑了几种不同的参数设置，以评估模型稀疏性、过分散性和读取深度水平对变量选择性能的影响。对于每个模型，进行50次重复，并通过以下测量的平均值评估变量选择性能

其中TP和FP分别是真阳性和假阳性的数量，以及给出了变量选择性能的总体度量。

4.2. 仿真结果

对于第一组模拟，我们检验了相关协变量比例对变量选择性能的影响。我们选择⁠,⁠,⁠，并考虑三个级别的模型稀疏性，⁠.其他参数选择为SNR=4.5，⁠,⁠,⁠,⁠图的顶部面板1给出了召回率、精确度和具有不同稀疏度级别的模型的值。当模型稀疏时(⁠⁠)，我们观察到这四种方法具有相似的性能，所提出的LNM略优于其他三种方法。作为增加并且模型变得不那么稀疏，所有四种方法的选择性能都会降低。然而，LNM模型显示出比其他方法更高的召回率。例如，当⁠LNM平均识别约75%的相关协变量，而LN、DM和ML模型仅分别识别约60%、28%和20%的相关变量。在精度方面，LN和LNM方法的精度随着增加而降低⁠然而，就整体衡量得分，LNM仍然优于其他三种方法。

接下来，我们评估过分散对模型选择性能的影响。我们通过将信噪比固定为1.5、4.5和⁠分别是。我们选择⁠,⁠,⁠,⁠,⁠,⁠,⁠,⁠结果显示在图的中间面板中1对于所有方法，较小的过分散会导致更好的变量选择性能。LNM方法的召回率高于其他准确率较低的方法。然而，LNM的总体变量选择性能由得分高于其他三种车型。这在很大程度上是由于使用交叉验证选择了调整参数，因为当过分散非常大时，CV往往会选择更大的模型来实现更好的预测。这导致在模型中选择更多的协变量，从而降低精确度和提高召回率。如果使用较小的调整参数，所提出的模型可以获得与我们比较的其他模型相似的精度。

最后，我们检验了不同测序深度数据的变量选择性能。请注意，测序深度决定了我们观察到的每个细菌分类群的计数数量。我们选择⁠,⁠,⁠,⁠，信噪比=4.5，⁠、和⁠.我们认为⁠结果显示在图的底部面板中1我们可以看到，LNM方法优于召回率较高的其他方法，在所有设置中，F1得分仅略低于准确率。此外，ML和DM方法的性能在很大程度上取决于分类群计数的总数。当测序深度较大时，LM和LNM方法的性能类似。

4.3. 计算复杂性、对MH样本的敏感性和MCEM的收敛性

为了评估MCMC实现中参数估计值对MH样本数的敏感性，我们用⁠,⁠,⁠、和⁠.系数矩阵的非零元素固定为2。此外，我们将信噪比设置为1.5⁠,⁠该模型对应于高水平的过度分散和大的多项式采样变化。我们应用五倍交叉验证来选择调谐参数⁠为了获得30–40%左右的平均接受率，将MH步长设置为1.5。图2显示了使用不同老化次数和不同MH运行次数对MCEM算法参数的估计，以接近EM算法中的预期值。我们观察到，参数估计值在几个EM迭代内保持稳定，并且估计值对E步骤中使用的烧入数量或MH样本数量不太敏感。我们考虑的其他模型也观察到了类似的结果。最后，我们还检查了初始值设置为不同的值。参数估计值对初始值不敏感（见图2).

5.将营养素与人体肠道微生物组分联系起来

我们应用LNM模型对饮食和粪便微生物组分之间的关系进行了横断面研究(Wu等人，2011年)。本研究的目的是研究长期饮食对人体肠道微生物组分的影响。在本研究中，收集了98名健康志愿者的粪便样本，并通过V1–V2区16S rRNA基因片段的454/Roche焦测序对DNA样本进行了分析。QIIME管线分析了火山岩序列(Caporaso等人，2010年)其中定义了3608个OTU。这些OTU可进一步合并为11门127属。此外，还使用食物频率问卷收集了这98名受试者的长期习惯性饮食信息，该问卷对214种微量营养素进行了定量测量。微量营养素测量值高度相关，如果相关系数大于0.90，则进一步分组。对于每个营养素组，选择了一种具有代表性的营养素，总共我们分析中的协变量。将这些微量营养素分组也可以提高结果的可解释性。例如，我们将钾和不添加补充剂的钾分为一个钾组，将原花青素三聚体、原花色素4-6聚合体和原花色素7-10聚合体合为一个原花色素组。

我们重点研究了微量营养素对拟杆菌属、普雷维特菌属和瘤胃球菌属三个细菌属组成的影响。这三个细菌属显示为Arumugam等人（2011年）将单个样本分为三类，因此被用于定义所谓的肠型。我们同样使用多维缩放和聚类分析在数据集中观察到了这样的聚类(Wu等人，2011年)。在这三个属中，在所有98个高丰度样本中观察到拟杆菌，仅在36个样本中观察到普雷氏菌，在73个样本中观察到瘤胃球菌。在含有拟杆菌的样本中，计数范围为216至11890，平均计数为3687。对于36个普氏菌样本，计数范围为1至9623，平均计数为1786。对于73个带有瘤胃球菌的样本，计数范围为2至728，平均计数为109。显然，这些计数因样品和属而异，差异很大。

我们申请了小组对该数据集进行惩罚LNM回归，将Prevetella作为对数比率转换中的基本成分。初始烧入后，在每个E步骤中使用2000 MH样品。对于每个调谐参数，参数估计在10次EM迭代内稳定。我们使用五次交叉验证来选择最终的调整参数。最终模型中总共选择了10种微量营养素（见表1这些营养素的名称）。为了评估这些营养素对肠型的影响，我们应用了加性logistic变换非零回归系数的估计值，以获得选定营养素扰动的估计值并将结果显示在表中1。这些结果与Wu等人（2011）使用简单的Spearman相关分析。例如，我们观察到普氏菌与麦芽糖、蔗糖、添加的小麦胚芽和维生素E/食品强化呈正相关，而拟杆菌与麦芽糖类、蔗糖、维生素E/食物强化呈负相关。相反，脂肪和氨基酸/胆碱营养素（包括棕榈酸反式脂肪酸、脯氨酸和总胆碱）与脯氨酸和胆碱的影响程度并不太大，但两者之间存在反向关联。

除了上述相关营养素外，LNM模型还确定了其他三种营养素，由于它们对这三个属的组成有较大的扰动影响，因此值得进一步研究。我们观察到丙烯酰胺和β-隐黄素与普氏菌呈正相关，而维生素C与拟杆菌呈正相关。丙烯酰胺是富含碳水化合物食物的天然副产品，它会显著增加普氏菌的含量。在某些类型的高温烹饪中，如油炸、烘焙和烘焙，糖和氨基酸（天冬酰胺）会形成丙烯酰胺。β-隐黄质和维生素C是肠道愈合补充剂中最常见的营养素(Claesson等人，2012年).

为了进一步证明营养素对肠型的影响，我们使用了图中所示的三元图三通过单纯形上的回归曲线来可视化营养物质的扰动模式。三元图长期以来一直用于显示组成数据，作为三部分组成的绘图区域（例如，拟杆菌、普雷维氏菌和瘤胃球菌的比例）。它可以被视为满足三部分组成的求和约束的正值平面(艾奇逊，1986年)。为了清楚地证明，我们假设三种肠型的基线组成位于单纯形的中心，并且只绘制了扰动程度较大的前九种营养素。该图清楚地表明，瘤胃球菌的组成相对于营养成分的变化是相当稳定的，而拟杆菌和普氏菌的相对丰度则随着微量营养物质的变化而急剧变化。

作为比较，我们还将ML模型和DM模型应用于该数据集。ML模型选择了五种营养素，蔗糖、脯氨酸、小麦添加胚芽、胆固醇和维生素B12，而DM模型选择了六种营养素，包括蔗糖、脯氨酸、小麦添加胚芽、胆碱（磷脂酰胆碱）、胆固醇和蛋氨酸。除了胆固醇、维生素B12和蛋氨酸外，LNM模型选择了所有这些营养素。然而，对肠型组成干扰最大的三种营养素，即丙烯酰胺、麦芽糖和棕榈酸反式脂肪酸，均未通过ML或DM模型进行选择。图4显示了LMN、ML和DM模型对拟杆菌和瘤胃球菌的拟合计数以及观察到的计数。除大约10个样本外，所有三个模型都很好地拟合了大多数数据。LMN模型似乎更符合拟杆菌数量。对于极低丰度的瘤胃球菌，由于观察到的计数较少，很难从这三种模型中的任何一种获得精确的拟合。

最后，为了评估结果对用于营养成分分组的相关截止值的敏感性，我们还通过将微量营养成分分组为营养组的相关系数为0.95。结果与使用0.90的相关截止值时的结果一致。LNM模型选择了七个营养组，其中六个营养组在使用0.90的相关截断值时也与之前一样被选择，包括棕榈酸反式脂肪酸、丙烯酰胺、维生素C、麦芽糖、蔗糖、添加小麦宝石。LNM还选择了营养组甜菜碱和胆碱，当使用0.90的相关截值时，这与LNM选择的总胆碱营养组重叠。

6.讨论

在本文中，我们提出了一个组LNM回归模型的惩罚估计，以选择与细菌组成相关的协变量。将每个协变量对应的系数作为一组处理，该模型结合分类群数据中的信息来选择干扰细菌组成在单纯形上整体位置的相关协变量。与计数数据的常用Dirichlet多项式回归模型相比，LNM模型提供了一种更灵活的建模细菌组成相关性的方法。我们选择LNM模型的一个特殊原因是，它允许逻辑正态分布的丰富结构来描述基团间协方差。更一般地说，我们的建模框架可以被视为多项式数据的随机效应模型。这种方法允许微生物组研究中典型数据的多项式外变异。我们使用MCEM算法估计模型，并在M步中进行惩罚估计。已经进行了广泛的模拟，将所提出的模型与其他回归模型进行比较，以选择相关的协变量。结果表明，该方法具有较好的变量选择性能。我们应用该模型分析了人类肠道微生物组数据，并确定了与定义人类肠道微生物群肠型的三个属的组成相关的营养素。

过多的零值会给分析微生物组分数据带来一些困难。这些零可以表示社区中真正不存在的组件（称为结构化零），也可以表示从社区中抽取的样本中不存在的罕见组件。如果没有额外的知识，这两种可能性是无法区分的。由于多项式分布允许零计数，我们的LMN模型可以处理零观测值，正如我们在肠道微生物组数据分析中所看到的那样。根据分析的目标，研究人员应该决定如何解释这些零计数，并相应地选择分析方法。在我们对微生物组数据的分析中，由于只考虑了三个属，因此零观测值的问题并没有那么严重。然而，当考虑数百或数千个细菌分类群时，我们预计数据集中会出现过多的零，因为许多细菌分类群非常罕见和稀疏。虽然所提出的MCEM算法原则上可以用于分析此类数据集，但在计算上更具挑战性。此外，如何更好地处理这种稀疏的成分数据还需要进一步研究。要处理所考虑的大量分类群，一个可能的解决方案是考虑稀疏群在总体回归系数矩阵中引入稀疏性的惩罚函数。该方法已用于陈和李（2013）。

7.补充资料

我们已经将该方法实现为R包“PenLNM_1.0”，该包现在可以在CRAN上使用(http://cran.r-project.org/web/packages/PenLNM/index.html)。本文附带的zip文件PenLNM_1.1.tr.gz可在Wiley Online Library的生物统计学网站上找到。

致谢

这项研究得到了美国国立卫生研究院CA127334和GM097505拨款以及香港研资局研究拨款（766511M）的支持。我们感谢评审人员和AE提供了许多非常有帮助的意见。

工具书类