使用一般线性模型进行多组差异比较，以及其与方差分析的联系

使用一般线性模型进行多组差异比较，以及其与方差分析的联系 by 小白鱼的生统笔记

使用一般线性模型进行多组差异比较，以及其与方差分析的联系

有时候我们会在一些文献中见到，作者使用某种回归模型去比较多组之间数值的差异显著性。是的，其实很多回归模型都可以用于差异分析，并不让人感到意外。

举个最简单的例子，一般线性模型（general linear model）和方差分析（Analysis of variance，ANOVA）。尽管二者具有独立起源和发展过程，但从函数形式上看，方差分析可以归类于回归模型的特例。在一般线性模型中，根据自变量的存在形式，可以分为简单和多元线性回归模型（自变量是连续变量时）和方差分析模型（自变量是类别或因子变量时），后者在某种程度上即等同于统计学中的方差分析。您不妨可以这样想，原则上二者最后都是通过计算F统计量获取对p值的估计，对于相同的数据而言，二者最后计算的F统计量肯定是相同的，所以不难想到二者进行差异分析时的效果也是相似的。

本篇就不妨通过一些简单的小例子，在R语言中演示通过一般线性模型进行多组差异比较的过程，并将它与方差分析的结果进行比较，以展示二者的联系。

下文中所使用的示例数据和R代码的百度盘链接（提取码，j0np）：

https://pan.baidu.com/s/1YGQPV0RODzzUMAb43jbkEw

若百度盘失效，也可在GitHub的备份中获取：

https://github.com/lyao222lll/sheng-xin-xiao-bai-yu-2021

以单因素方差分析和一元线性模型为例的比较

首先来看一个最简单的例子，只存在单因素的情况。在一元线性模型中，当自变量是类别（或因子）变量时，等效于单因素方差分析。以下展示其在R中的实现过程，并据此作详细解读。

示例数据集概要

已知啮齿动物神经损伤后，脊髓的背角神经元导致疼痛超敏反应，而参与神经疼痛信号的分子和信号传导级联是复杂的，其中可能涉及某些重要的circRNA。Zhang等（2019）通过脊神经结扎手术（SNL），在大鼠脊髓背角组织中诱导神经损伤，通过qPCR验证了脊髓特异性circAnks1a在神经损伤的组织中显著更高表达，并进而阐明该分子参与神经损伤引起的中枢敏化和疼痛行为。

网盘文件“circAnks1a_qPCR.txt”来自Zhang等（2019）的图1d源数据（您也可以在原文献的补充材料中获取，注意我只是借用该数据进行本推文的演示，不代表原文献也用了这些方法），通过qPCR定量了脊神经结扎手术（SNL）前后的大鼠脊髓背角组织中，circAnks1a在各个时间点的相对表达水平。Sham，SNL前的组织；SNL 3d、7d、10d、14d分别为SNL后的第3、7、10、14天的组织。

我们希望比较脊神经结扎手术（SNL）前后不同时间点的大鼠脊髓背角组织中， circAnks1a的相对表达水平高低情况，以观察circAnks1a是否在由SNL诱导的神经损伤的组织中存在更高的表达，以及是否随损伤时间的增加更明显。

通过方差分析比较circAnks1a在SNL前后各时间点的表达水平高低

对于这种单个响应变量在多组间的比较情形，我们不难下意识想到通过单因素方差分析结合多重比较，直接比较响应变量在多组间的均值差异，如下所示。

#读取 circAnks1a 的 qPCR 定量数据
circAnks1a <- read.delim('circAnks1a_qPCR.txt', stringsAsFactors = FALSE)
 
#预定义因子水平，也就是给分组按时间排个序
circAnks1a$Group <- factor(circAnks1a$Group, levels = c('Sham', 'SNL 3d', 'Day 7d', 'Day 10d', 'Day 14d'))
 
##首先需要对响应变量进行正态性或方差齐性评估，满足这些假设才能进行方差分析
#可使用 car 包 QQ-plot 检查数据是否符合正态分布（所有的点都离直线很近，落在置信区间内说明正态性良好）
library(car)
qqPlot(lm(Relative_expression~Group, data = circAnks1a), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
 
#可使用 Bartlett 检验进行方差齐性检验（p 值大于 0.05 说明等方差）
bartlett.test(Relative_expression~Group, data = circAnks1a)
 
##上述两个假设通过，执行单因素方差分析（One Way ANOVA）
fit_aov <- aov(Relative_expression~Group, data = circAnks1a)
summary(fit_aov)  #查看方差分析概要
 
##方差分析后，使用 Tukey HSD 检验进行事后多重比较，继续探寻两两分组间的差异
#multcomp 包提供了直观的结果比较
library(multcomp)
 
tuk_aov <- glht(fit_aov, alternative = 'two.sided', linfct = mcp(Group = 'Tukey'))
summary(tuk_aov)  #查看 Tukey HSD 概要
 
tuk_aov.cld <- cld(tuk_aov, level = 0.05, decreasing = TRUE)
tuk_aov.cld  #按差异高低自动标识了 a、b、c 水平
plot(tuk_aov.cld, col = c('#7CBAD8', '#E1EBF4', '#D3E599', '#F9D2B8', '#E7CBE2'))

这里的单因素方差分析通过计算F统计量估计显著性p值。根据整体的p<0.05（F=14.53，p=4.35×10^-7），可以认为在通过脊神经结扎手术（SNL）诱导的大鼠神经损伤的组织中，与SNL前的非损伤组织相比，circAnks1a存在更高的表达。同时继续根据Tukey HSD检验的结果可知，circAnks1a在大鼠神经损伤组织中的表达随损伤时间的升高更明显。

通过线性模型比较circAnks1a在SNL前后各时间点的表达水平高低

正如本篇一开始时提到，很多回归方法也可用于组间差异分析。在这里，我们尝试将上述数据带入一般线性模型（general linear model），探讨circAnks1a表达随脊神经结扎手术（SNL）诱导的大鼠神经损伤时间的关系。由于数据中，自变量（包括SNL前的1个时间点以及SNL后的4个时间点）已定义为一组因子变量，此时将执行一个“类别或因子型自变量的线性模型”。

#读取 circAnks1a 的 qPCR 定量数据
circAnks1a <- read.delim('circAnks1a_qPCR.txt', stringsAsFactors = FALSE)
 
#预定义因子水平，也就是给分组按时间排个序
circAnks1a$Group <- factor(circAnks1a$Group, levels = c('Sham', 'SNL 3d', 'Day 7d', 'Day 10d', 'Day 14d'))
 
##一般来说，一般线性模型也对响应变量的正态性和方差齐性有要求，但在很多情况下经常被忽略
#如果您期望评估这两种假设条件，方法和上文方差分析一样
#可使用 car 包 QQ-plot 检查数据是否符合正态分布（所有的点都离直线很近，落在置信区间内说明正态性良好）
library(car)
qqPlot(lm(Relative_expression~Group, data = circAnks1a), simulate = TRUE, main = 'QQ Plot', labels = FALSE)
 
#可使用 Bartlett 检验进行方差齐性检验（p 值大于 0.05 说明等方差）
bartlett.test(Relative_expression~Group, data = circAnks1a)
 
##类别（或因子）型自变量的一元线性模型
fit_lm <- lm(Relative_expression~Group, data = circAnks1a)
summary(fit_lm)  #查看线性模型概要
 
##对于类别（或因子）型自变量的一元线性模型，仍可使用 Tukey HSD 检验进行事后多重比较，继续探寻两两分组间的差异
#仍然可使用 multcomp 包中的方法执行 Tukey HSD 检验
library(multcomp)
 
tuk_lm <- glht(fit_lm, alternative = 'two.sided', linfct = mcp(Group = 'Tukey'))
summary(tuk_lm)  #查看 Tukey HSD 概要
 
tuk_lm.cld <- cld(tuk_lm, level = 0.05, decreasing = TRUE)
tuk_lm.cld  #按差异高低自动标识了 a、b、c 水平
plot(tuk_lm.cld, col = c('#7CBAD8', '#E1EBF4', '#D3E599', '#F9D2B8', '#E7CBE2'))

在这个简单一元线性模型中，同样通过计算F统计量估计显著性p值。如果仍基于这两个统计量判断的话，F=14.53和p=4.35×10^-7，也即获得了和上文单因素方差分析一致的结论：在通过脊神经结扎手术（SNL）诱导的大鼠神经损伤的组织中，与SNL前的非损伤组织相比，circAnks1a存在更高的表达。同时继续根据Tukey HSD检验的结果可知，circAnks1a在大鼠神经损伤组织中的表达随损伤时间的升高更明显。

当线性模型的自变量为类别或因子型时，R函数lm()是如何处理的？

在R中，对于运行一般线性模型的函数lm()而言，大多数情况下我们输入的自变量通常都是数值型变量。这种情况也好理解，将依据普通最小二乘法（Ordinary Least Square，OLS）的原则执行回归。但是，如上这种情况，当自变量是类别或因子型变量时，lm()如何工作的呢？

在前文“类别或因子型自变量的线性回归”中提到，当lm()函数碰到类别或因子型自变量时，它会用一系列与因子水平相对应的数值型名义变量（dummy variables）来代替原有的因子。如果因子有k个水平，将会创建k–1个名义变量。在本示例中，因子变量trt有5个水平，即Sham、SNL 3d、Day 7d、Day 10d、Day 14d。对应于此，R构建了4个（等级数-1，即5-1=4）数值型名义变量代替原始的因子变量，以执行本该执行的回归操作。

#查看自变量预定义的因子水平和顺序
levels(circAnks1a$Group)
 
#通过 contrasts() 查看名义变量的编码过程
contrasts(circAnks1a$Group)
 
#查看上文线性模型的结果概要
summary(fit_lm)

在上述概要中，4个名义变量显示为GroupSNL 3d、GroupDay 7d、GroupDay 10d、GroupDay 14d，它们分别对应于contrasts()所示的编码表的4列。对于解读方式，详情也可见前文“类别或因子型自变量的线性回归”。例如，当关注通过脊神经结扎手术（SNL）诱导3天后的大鼠神经损伤的组织SNL 3d时，名义变量GroupSNL 3d等于1，其它名义变量GroupDay 7d、GroupDay 10d、GroupDay 14d都等于0，以此类推带入模型中获取参数估计。其余情况亦是如此。

其它类型的方差分析和线性模型的关系

上述是最简单的情形，单因素方差分析和一元线性模型，对应了变量之间的简单一元响应关系。对于其它情形也是类似的，如：协方差分析和多元线性模型，双（多）因素方差分析和带交互因素的线性模型，重复测量方差分析和线性混合效应模型，等等。这些就不再细说了。

其实如果只是仅基于F统计量获取对差异显著性p值的估计，无论使用方差分析或一般线性模型，获得的结果都是一致的。这种情况下貌似也没必要考虑那么多，直接使用方差分析就完全可以了，除非您有其它的重要参数获取。

参考资料

Robert I. Kabacoff. R语言实战（第二版）（王小宁刘撷芯黄俊文等译）. 人民邮电出版社, 2016.

Zhang SB, Lin SY, Liu M, et al. CircAnks1a in the spinal cord regulates hypersensitivity in a rodent model of neuropathic pain. Nat Commun. 2019;10(1):4119.

友情链接

数据预处理

常见的数据标准化或转化方法

常见的参数检验类型

两组比较：t检验（t-test）

两组多元响应：Hotelling's t-squared test

多组单因素比较：单因素方差分析（单因素ANOVA）+多重比较

多组协方差比较：单因素协方差分析（ANCOVA）

多组双因素比较：双因素方差分析（双因素ANOVA）

多组混合效应：重复测量方差分析

多组多元响应：多元方差分析（MANOVA）和稳健多元方差分析（稳健MANOVA）

常见的非参数检验类型