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

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

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

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

首先来看一个最简单的例子，只存在单因素的情况。在一元线性模型中，当自变量是类别（或因子）变量时，等效于单因素方差分析。以下展示其在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.

小结

本文原创作者：小白鱼，请支持原创！

感谢大家耐心看完，自己的文章都写的很细，代码都在原文中，希望大家都可以自己做一做，请关注后私信回复“数据链接”获取所有数据和本人收集的学习资料。如果对您有用请先收藏，再点赞转发。

也欢迎大家的意见和建议，大家想了解什么统计方法都可以在文章下留言，说不定我看见了就会给你写教程哦。

如果你是一个大学本科生或研究生，如果你正在因为你的统计作业、数据分析、论文、报告、考试等发愁，如果你在使用SPSS,R，Python，Mplus, Excel中遇到任何问题，都可以联系我。因为我可以给您提供好的，详细和耐心的数据分析服务。

如果你对Z检验，t检验，方差分析，多元方差分析，回归，卡方检验，相关，多水平模型，结构方程模型，中介调节，量表信效度等等统计技巧有任何问题，请私信我，获取详细和耐心的指导。

If you are a student and you are worried about you statistical #Assignments, #Data #Analysis, #Thesis, #reports, #composing, #Quizzes, Exams.. And if you are facing problem in #SPSS, #R-Programming, #Excel, Mplus, then contact me. Because I could provide you the best services for your Data Analysis.

Are you confused with statistical Techniques like z-test, t-test, ANOVA, MANOVA, Regression, Logistic Regression, Chi-Square, Correlation, Association, SEM, multilevel model, mediation and moderation etc. for your Data Analysis...??

Then Contact Me. I will solve your Problem...

加油吧，打工人！

猜你喜欢

R数据分析：逻辑斯蒂回归与泊松回归

R数据分析：多分类逻辑回归

R数据分析：逐步回归的做法和原理，案例剖析

R数据分析：一般线性回归的做法和解释

R数据分析：多元逻辑斯蒂回归的做法