广义线性模型应用举例之beta回归及R计算

比例数据、beta分布与beta回归

在生物学数据分析中，经常会涉及到比例型的数据类型。例如，试验成功的概率（试验成功次数/总次数），疾病发生频率（人群中患者数/总人数），物种相对丰度（物种个体数量/群落物种总数）等。比例数据的特点是分布在区间(0, 1)范围内的连续型变量，并且实现了对样本整体某种形式的归一化。    

常见的回归方法应用至比例数据时可能存在局限性

常规的线性回归假定响应变量服从正态分布。比例数据是连续型的数值变量，并且有时一些比例数据也表现出良好的正态分布，似乎可以通过线性回归进行分析，并将获得较合理参数估计以及显著性水平。但更多时候，数值分布大都存在较高的偏离正态性程度，使得线性回归分析比例数据的结果并不理想，回归系数和p值等比较糟糕，变量的生物学意义难以解释。并且基于正态分布的线性回归得出的预测值和置信区间很可能包含比例数据定义区间外的值，即拟合出区间(0, 1)范围外的偏离实际的值。

某些比例数据来自计数型数值的转换。如在群落研究中调查物种丰度时，观察到的物种个体数量是离散型的非负整数，对计数型物种丰度数据的建模，通常泊松回归或负二项回归是优先考虑的。但有时出于特定需要，可能会转换为物种相对丰度（物种个体数量/群落物种总数）用于后续统计分析，此时泊松或负二项回归的局限在于无法应用至小数数值的比例数据中。

比例数据可以为概率响应（如试验成功的概率），有时形似二项分布试验，容易联想到基于二项分布的logistic回归。但logistic回归要求响应变量为0-1的二分型或非数值的类别型，故不适合连续数值的比例数据分析。    

beta分布与beta回归

beta回归是一种分析比例数据的出色替代方案，其假定响应变量服从beta分布，对于分布在区间(0, 1)范围内的连续型响应变量的建模非常有用。

beta分布

首先来看beta分布，具有beta分布的变量的值是介于0到1之间的连续变量，beta分布的密度函数为：

式中p>0且q>0，Γ(·)是伽马函数（gamma function），beta分布密度函数的均值E(y)和方差var(y)分别为：

beta分布的密度曲线随均值和方差参数的变化而呈现出许多不同的形状。如下所示了几种代表样式，变量的所有值介于0到1之间，并且分布形状覆盖了各异的均值和离散度水平，而实际上很多生物学数据的变量分布都具有这种beta分布状态的特点。因此，beta分布覆盖面相当广，基于beta分布的beta回归在很多实际问题中有着广泛的应用。

也可以看到，beta分布中某些状态大致与正态分布相同，例如A样式。常规的线性回归假定响应变量服从正态分布，对于呈现A样式beta分布的响应变量，或许可以通过一般线性回归进行参数估计。但若是D、E、F这种，线性回归将会产生较大的偏差。

beta回归

因此，在回归分析中，若响应变量为比例数据，引进beta回归模型进行统计分析顺其自然，并且beta回归模型通常能够较好地对其进行拟合。

beta分布是指数分布族的一员，推广到广义线性模型中就是beta回归。若响应变量Y服从beta分布，则对于给定的一组自变量xiT（xt1, …, xti），Y的条件均值ut可写作广义线性回归模型：

βi为各自变量xti的回归系数。连接函数g(.)存在多种选择，一个最为常见的连接函数是logit连接（logit link），在这种情况下：

xiT即各自变量（xt1, …, xti），β为各自变量xti的回归系数（β1, …, βi）。

关于beta回归其它细节，可参考Ferrari和Cribari-Neto（2004）的阐述。

文献中使用beta回归的一个实例

上文提到，计数型变量在生物学数据中非常常见，例如物种丰度计数。在广义线性模型中，泊松回归或负二项回归是对这类计数型响应变量建模的常见方法。但有些情况下可能需要转换为物种相对丰度的比例数据进行统计分析，此时泊松回归或负二项回归不适用，就可以考虑beta回归。

因此，不妨首先来看一篇与植物根系微生物组有关研究，将文中涉及到beta回归的部分内容节选出来以帮助大家理解。

Edwards等（2018）研究了连续种植季水稻的根系微生物组。正文Result的第二部分内容中，作者通过特定门水平细菌类群的相对丰度，描述了根系相关微生物在水稻根系的富集模式，包括在土壤到根系区域的空间微尺度层面，以及随水稻根系生长发育的时间层面。

为了明确在细菌门水平上，是否存在特定微生物类群表现出被水稻根系富集或排斥的模式，使用beta回归对它们由土壤-根际-根表-根内（由根外到根内）的相对丰度进行了建模（图2B）。beta回归通过R包betareg执行，各细菌门的相对丰度在这里就作为响应变量对待，选择beta回归的原因就是考虑到各微生物的相对丰度都代表了其来自总微生物群落的相对比例，都是分布在 (0, 1) 区间内的比例数值，beta回归将更加有效。土壤、根际、根表、根内这4个空间微尺度区域是一组有序的类别型自变量，作者对它们分配了一组有序的离散值指代：土壤0、根际1、根表2、根内3。在55个细菌门中，共确定了从根的外部到内部空间分布明显不同的42个门（beta回归中的系数显著）。从土壤环境到根内部，相对丰度显著增加的有20个门（beta回归系数大于0），代表了它们被水稻根系富集。从土壤环境到根内部，相对丰度降低的有22个门（beta回归系数小于0），代表了它们被水稻根系排斥。根系排斥的细菌门的回归系数的绝对值总体上更大，表明只有少数比例的土壤微生物由根系富集。

作者还在4个空间微尺度区域内，分别进行了细菌门的相对丰度与水稻年龄的beta回归，以确定哪些微生物表现出明显的随水稻生长的响应（图2C）。观察到不同根系区域中的每个细菌门的丰度随时间表现出近乎一致的趋势。在根表和根内，大多数具有动态丰度模式的细菌类群在水稻生长时期中显著减少（分别为12/22和32/40），表明在根系发育初期存在大量土壤微生物随机定殖在根系，但最终在与其它类群竞争或者寄主植物选择中被排斥。

原文图2部分。Bulk Soil表示土壤，Rhizosphere表示根际，Rhizoplane表示根表，Endosphere表示根内；Plant Age表示水稻生长时间（天）。

B，beta回归探究不同的细菌门在土壤-根际-根表-根内（由根外到根内）的富集模式，棒棒糖图展示回归系数，回归系数>0为根系富集，0为根系富集，%”管道符，对于管道符这玩意，白鱼同学的感受一直就是，方便自己书写但不方便别人解读。所以下文就没直接套用作者的原代码，而是作了拆分和更改，虽然看着比先前啰嗦了点但是更方便理解分析过程，并且运行结果不变。

响应变量分布状态的探索

数据集中有几十类细菌门待进行趋势分析，提示要进行批处理。不过在批处理之前，先以其中某个类群为例，展示R函数的运行方法，以及对结果的解读。例如，期望探究Acidobacteria（酸杆菌门）相对丰度是否由土壤到根系发生富集或降低。