R语言中简单线性回归和多项式回归分析简介：结合实例加深常见线性回归应用认识

一个简单线性回归和多项式回归在R中的实现示例

常见线性回归的原理无需多言，大家都知晓。它就是普通最小二乘法（Ordinary Least Square，OLS）回归法，此方法主要被应用于简单线性回归、多项式回归以及多元线性回归。

在 R 语言里，用于拟合线性回归的函数是 lm()，这想必大家都不陌生。实际上，线性回归方程本身似乎没什么可讲的，它的执行式也很容易理解。关键在于如何使用它，它适用于哪些场合，能拟合哪些变量间的关系，以及如何根据变量的结构特点选择合适的回归模型等，这些才是最实际的需求。所以，白鱼同学会尽可能找一些有趣的实例，以帮助大家加深对这些的认识。

本篇首先从简单线性回归开始，接着介绍如何在 R 中执行回归分析。然后从多项式回归开始，继续介绍在 R 中执行回归分析的相关内容。

某文献中的回归示例（部分节选）

来看文献“Bacterioplankton Biogeography of the Atlantic Ocean: A Case Study of the Distance-Decay Relationship”中的部分内容。

描述大西洋表层浮游细菌群落相似度的距离衰减模式

如下图为原文图 3B 所示，横轴是采样地之间的地理距离，纵轴代表着采样地间 40m 水深处 FL 浮游细菌群落的组成相似度。其中，South 代表单独考虑南半球时的群落相似度距离 - 衰减趋势，North 代表单独考虑北半球时的群落相似度距离 - 衰减趋势，Both 代表综合考虑南北半球时的群落相似度距离 - 衰减趋势。

备注：一些样本与边界地区（大西洋中靠近最南和最北地带）有关，这些样本用红色标注。作者在分析中会对它们进行特殊对待，下文会提及。

接下来以该部分分析为例，以采样地点间的群落相似度与地理距离的关系为例，展示在 R 中如何拟合线性回归方程。

示例数据和预处理过程

原文献附表“Table S3”里，作者给出了研究中采样地点的地理位置（经纬度）以及浮游细菌物种组成等数据，这些数据可通过搜索原文来找到并进行下载。

白鱼同学整理了一小部分，他挑选出了“40m 深度的 FL 浮游细菌物种组成丰度表”以及“对应采样地点的经纬度坐标”，这些是用于重现上述 3 图中回归分析的内容。

整理后的数据可见网盘附件（提取码，3m3l）：

节选的部分示例数据

网盘附件“site.txt”中记录了采样地点的地理位置。样本名称在第一列。40 这个数字代表了水深为 40m 处的水体样本。Lat 指的是纬度，Lon 指的是经度。Province 表示采样地所属的地区。需要特别注意的是，NADR 和 FKLD 分别是最北和最南部地区。如上文所提到的，作者考虑到它们可能会对整体趋势产生影响，所以在执行回归时特别进行了考虑，具体详情可详见下文过程。

“FL.depth40.otu_table.txt”是这些采样地点在 40m 深度的水体样本所对应的 FL 浮游细菌物种组成的丰度表。其中，每一行为一种微生物的 OTU，而每一列代表着一个样本，列的名称就是采样地的名称。

计算样点的地理距离及群落组成相似度

为了拟合采样地间 FL 浮游细菌群落组成相似度与地理距离的线性方程，首先得：

根据各地点的经纬度坐标记录，将其用于计算采样地之间的地理距离。

根据各地点的 FL 浮游细菌物种组成丰度信息，来计算采样地之间 FL 浮游细菌群落的组成相似度。

##根据经纬度计算采样点的地理距离
#读取采样点的地理位置数据
site <- read.delim('site.txt', sep = '\t', row.names = 1, check.names = FALSE)
#site <- site[which(site$'Lat (°N)' > 0), ]  #只选择北半球
#site <- site[which(site$'Lat (°N)' < 0), ]  #只选择南半球
 
#计算采样点间的地理距离
geosphere 包中的 distm()函数能够依据经纬度来计算地理距离，其默认的距离单位是米。
#distm() 要求两列数据，第一列是经度，第二列是纬度
site_dis <- geosphere::distm(site[c('Lon (°E)', 'Lat (°N)')])
rownames(site_dis) <- rownames(site)
colnames(site_dis) <- rownames(site)
 
#将采样点地理距离矩阵转换为两两对应数值的数据框结构
site_dis <- reshape2::melt(site_dis)
site_dis <- subset(site_dis, value != 0)
head(site_dis)
 
计算群落间的物种组成相似度，依据原文的描述，运用 Bray-curtis 相似度。
#浮游细菌群落物种组成数据
spe <- read.delim('FL.depth40.otu_table.txt', sep = '\t', row.names = 1, check.names = FALSE)
spe <- data.frame(t(spe))
spe <- spe[rownames(site), ]  #与 site 中预先选择的南北半球数据对应
 
vegan 包中的 vegdist()函数用于计算群落间物种组成的 Bray-curtis 相异度矩阵
通过 1 - Bray - curtis 相异度来获取 Bray - curtis 相似度。
comm_sim <- 1 - as.matrix(vegan::vegdist(spe, method = 'bray'))
 
#将矩阵转换为两两群落对应数值的数据框结构
diag(comm_sim) <- 0  #去除群落相似度矩阵中的对角线值，它们是样本的自相似度
comm_sim[upper.tri(comm_sim)] <- 0  #群落相似度矩阵是对称的，因此只选择半三角（如下三角）区域的数值即可
comm_sim <- reshape2::melt(comm_sim)
comm_sim <- subset(comm_sim, value != 0)
head(comm_sim)
 
##采样点距离和群落相似度数据合并
comm_dis <- merge(comm_sim, site_dis, by = c('Var1', 'Var2'))
names(comm_dis) <- c('site1', 'site2', 'comm_sim', 'site_dis')
 
标识与最北的 NADR 地区采样地点有关的距离对，以及标识与最南的 FKLD 地区采样地点有关的距离对
site_edge <- rownames(subset(site, Province %in% c('NADR', 'FKLD')))
comm_dis 中，当 as.character(comm_dis$site1) 属于 site_edge 或者 as.character(comm_dis$site2) 属于 site_edge 时，其对应的 comm_dis[which(as.character(comm_dis$site1) %in% site_edge | as.character(comm_dis$site2) %in% site_edge),'edge'] 这部分内容。<- '0'
在 comm_dis 中，找出那些 edge 列的值不属于 '0' 的行，然后提取出这些行的 edge 列的值。<- '1'
head(comm_dis)
 
使用 write.table 函数将 comm_dis 写入到 'comm_dis.depth40.Both.txt' 文件中，分隔符为 '\t'，不显示行名，不使用引号进行包裹
使用 write.table 函数将 comm_dis 写入到 'comm_dis.depth40.North.txt' 文件中，以制表符 '\t' 作为分隔符，不包含行名，不使用引号进行包裹
使用 write.table 函数将 comm_dis 写入到 'comm_dis.depth40.South.txt' 文件中，以制表符 '\t' 作为分隔符，不包含行名，不进行引号引用。

如果选择南半球，就输出文件“comm_dis.depth40.South.txt”，内容如上；如果选择北半球，就输出文件“comm_dis.depth40.North.txt”，内容如上；如果选择南北半球整体，就输出文件“comm_dis.depth40.Both.txt”，内容如上。

其中，一个采样地点名为 site1，另一个采样地点名为 site2；二者的 FL 浮游细菌群落的组成相似度为 comm_sim，是通过 Bray-curtis 相似度指数来衡量的；二者的地理距离是 site_dis，单位是米（m）；如果 site1 和 site2 中有一个采样地点来自地区 NADR 或 FKLD，就标记为 1，否则标记为 0。

这里使用，所得结果和原文中计算的地理距离有一些微小差异，数值不完全相同，可能是在考虑地球椭球形状时出现了估计偏差。所以，下文线性回归的拟合 R2 也不完全一样，其实差别很小，整体趋势是一致的，无需见怪……

一个简单线性回归示例

OLS 回归表达式的形式大致为“Y =β1X+β0”，这是最简单的一元一次线性回归结构。其中，β1 代表斜率，β0 代表截距。

接下来，在 R 中进行操作，拟合采样地之间 FL 浮游细菌群落的组成相似度与采样地之间的地理距离的关系，以理解其用法。

作者在文中把南北半球进行了分开讨论，并且将其归因于不同的生物地理学方面的驱动力。

北大西洋FL浮游细菌群落的距离-衰减关系

首先考虑北半球的情形。

作者观察到，来自最北部地区 NADR 的采样地的群落结构没有明显的特殊性。并且，这些采样地与其他地点的相似度也没有明显的偏离趋势（图中红色的点）。因此，作者使用文件中所有数据直接进行了拟合。

在R中，直接在lm()中以lm(Y~X)的方程式指定。

#读取上述输出文件，首先以北半球的为例
comm_dis <- read.delim('comm_dis.depth40.North.txt', sep = '\t')
 
注：site_dis 列的单位为米（m），在这种大尺度的范围当中，将千米（km）作为度量会更为适宜一些。
comm_dis$site_dis_km <- comm_dis$site_dis/1000
 
#lm() 用于拟合群落组成相似度（comm_sim）与地理距离（site_dis_km）之间的一元线性关系
fit <- lm(comm_sim~site_dis_km, data = comm_dis)
summary(fit)  #展示拟合模型的简单统计

这是lm()拟合回归的默认输出，关注几个重要的值即可。

拟合的一元线性方程式为“Y = -0.00007684X + 0.5706”，且该方程式是显著的。R2 约为 0.4，这个数值不算低，这表明北大西洋 40m 水深处的 FL 浮游细菌群落存在明显的线性距离衰减趋势。

期望从中提取关键信息，或者进行一些其他统计时，这有一些参考用法。

#例如
coefficients(fit)[1]  #获取截距
获取 site_dis_km 的斜率，需通过 coefficients(fit)[2] 来实现。
 
可以在 names(summary(fit))之后查看主要的内容项，接着从中进行提取，例如
summary(fit) 的 adj.r.squared  #校正后 R2

最后，绘制出散点图，以此来展示北半球大西洋 FL 浮游细菌群落的距离与衰减之间的关系。

#ggplot2 绘制带线性拟合线的散点图
拟合线既考虑了“黑点”，也考虑了“红点”。“红点”代表着与 NADR 地区中采样地点相关的距离关系。
library(ggplot2)
 
p <- ggplot(comm_dis, aes(site_dis_km, comm_sim)) +
使用 aes 函数并将 color 参数设置为 as.character(edge)，然后与 geom_point 进行组合。
设置颜色手动取值，值为 c('red', 'black')，并设置颜色范围限制为 c('0', '1')，然后加上
使用“lm”方法进行平滑处理（geom_smooth）且不显示标准误差（se = FALSE）。
主题（面板网格 = 元素空白()，面板背景 = 元素矩形（颜色 = '黑色'，填充 = '透明'））
legend 的位置为“none”，绘图的标题为具有水平居中文本属性（hjust = 0.5）的元素文本。再加上
设置 y 轴的连续刻度范围为从 0 到 1 。
设置 x 为“Distance (km)”，y 为“Bray-curtis similarity”，标题为“North”的 labs 操作。
 
#提取上述结果中重要的统计值，用于在图中将方程式标注出
label_fit <- data.frame(
formula 等于通过 sprintf 函数生成的字符串，该字符串为“italic(Y) == %.7f*italic(X) + %.3f”，其中 %.7f 对应的是 coefficients(fit)[2]的值，%.3f 对应的是 coefficients(fit)[1]的值
R2 等于 sprintf 函数生成的字符串，该字符串内容为“italic(R^2) == %.3f”，其中的 %.3f 会被 summary(fit)$adj.r.squared 的值所替换，即 R2 等于“italic(R^2) == 某个具体的三位小数数值”，这个具体的数值是 summary(fit)$adj.r.squared 的值。
    p_value = sprintf('italic(P) < %.3f', 0.001))
 
p +
geom_text 的 x 坐标为 3000，y 坐标为 0.9，其标签由公式决定，数据来自 label_fit，并且开启了解析功能。
hjust 的值为 0，颜色为“黑色”，显示图例被设置为 FALSE 。
geom_text 的 x 坐标为 4000，y 坐标为 0.8，其标签由 R2 决定，数据来源于 label_fit，并且开启了 parse 功能。
hjust 等于 0，颜色为“黑色”，显示图例为假。 +
geom_text 中的 x 坐标为 4000，y 坐标为 0.7，其标签由 p_value 决定，数据来源于 label_fit，并且开启了 parse 功能。
    hjust = 0, color = 'black', show.legend = FALSE)

南大西洋FL浮游细菌群落的距离-衰减关系

接下来讨论南半球的情形。

作者发现，来自最南部地区 FKLD 的采样地点的 FL 浮游细菌群落结构与南半球其它采样地的群落结构有明显区别，这种区别在图中以红色的点呈现。所以在拟合方程式时，未考虑与 FKLD 中样本有关的距离关系。

为什么 FKLD 中的样本偏离明显呢？可能是因为它的地理位置特殊以及环境特殊。作者在文中进行了一些讨论，但是由于篇幅原因，只节选了一小部分来帮助理解。

FKLD 地区与南半球大西洋的其它采样地点相比，处于多个洋流带的交界处，它具有明显区别于其他地区的盐度和温度环境，其独特的水文特征或许是塑造更为别致的群落结构的原因。

首先选择剔除数据中与 FKLD 有关的距离关系，接着用剩余数据在 R 函数 lm()中按照 lm(Y~X)的方程式进行拟合。

#读取上述输出文件，南半球的样点距离与群落相似度数据
comm_dis <- read.delim('comm_dis.depth40.South.txt', sep = '\t')
 
#注：site_dis 列的单位是米（m），在这种大尺度范围下，以千米（km）作为度量会更好一些
comm_dis$site_dis_km <- comm_dis$site_dis/1000
 
#lm() 拟合群落组成相似度（comm_sim）地理距离（site_dis_km）的一元线性关系
#拟合时不再考虑与 FKLD 地区中采样地点有关的距离关系
comm_dis_select <- subset(comm_dis, edge == '1')
 
fit <- lm(comm_sim~site_dis_km, data = comm_dis_select)
summary(fit)  #展示拟合模型的简单统计

拟合的一元线性方程式为“Y = -0.0001032X + 0.8098”，且该方程式是显著的。R2 约为 0.67，这个数值很高，表明南大西洋 40m 水深处 FL 浮游细菌群落存在明显的线性距离衰减趋势。

绘制散点图，用来展示南半球大西洋 FL 浮游细菌群落的距离与衰减之间的关系。

#ggplot2 绘制带线性拟合线的散点图
拟合线仅将“黑点”纳入考量，“红点”体现了与 FKLD 地区中采样地点相关的距离关系，然而因其偏离程度过大而被剔除。
library(ggplot2)
 
p <- ggplot() +
使用 comm_dis 作为数据，通过 geom_point 函数进行绘图。在 aes 映射中，x 轴为 site_dis_km，y 轴为 comm_sim，颜色由 as.character(edge)来决定。最后加上相应的绘图操作。
scale_color_manual(values = c('red', 'black'), limit = c('0', '1')) +
geom_smooth 对 comm_dis_select 这个数据进行操作，将其 x 轴设为 site_dis_km，y 轴设为 comm_sim，使用方法为 'lm'，并且不显示标准误差（se = FALSE）。
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'),
    legend.position = 'none', plot.title = element_text(hjust = 0.5)) +
scale_y_continuous(limit = c(0, 1)) +
设置 x 为“Distance (km)”，y 为“Bray-curtis similarity”，标题为“South”的 labs 操作。
 
#提取上述结果中重要的统计值，用于在图中将方程式标注出
label_fit <- data.frame(
formula 等于通过 sprintf 函数生成的字符串，该字符串的内容为“italic(Y) 等于 %.6f 乘以 italic(X) 加上 %.3f”，其中 %.6f 被替换为 coefficients(fit)[2]，%.3f 被替换为 coefficients(fit)[1]
    R2 = sprintf('italic(R^2) == %.3f', summary(fit)$adj.r.squared),
    p_value = sprintf('italic(P) < %.3f', 0.001))
 
p +
geom_text 的 x 坐标为 4000，y 坐标为 0.9，其标签由 formula 决定，数据来自 label_fit，且 parse 为 TRUE。
    hjust = 0, color = 'black', show.legend = FALSE) +
geom_text 的 x 坐标为 5000，y 坐标为 0.8，其标签由 aes 中的 R2 决定，数据来源于 label_fit，并且 parse 为 TRUE。
    hjust = 0, color = 'black', show.legend = FALSE) +
geom_text 的 x 坐标为 5000，y 坐标为 0.7，其标签由 p_value 决定，数据来源于 label_fit，并且需要进行解析（parse = TRUE）。
    hjust = 0, color = 'black', show.legend = FALSE)

一个多项式回归示例

多项式指的是预测变量呈现多次方的形式，像“Y~X2”以及“Y~X3”等这样的情况。并且，多项式等式依然被认为是线性回归模型，原因在于等式仍然是预测变量的加权和形式。

继续通过该文的示例理解它的用法。

综合考虑南北大西洋FL浮游细菌群落的距离-衰减关系

上述分析表明南、北半球大西洋中的群落相似度距离衰减呈递减趋势。然而，当作者综合考虑南北半球时发现，FL 浮游细菌群落的组成相似度并非随着采样地间地理距离的增加而持续降低，在距离达到 6000km 之后，其相似度反而会上升。

此时，二者间的关系难以再用“Y =β1X+β0”这种形式进行描述。作者还发现，直接拟合简单线性回归的 R2 值是很低的。整体呈现先下降后上升的趋势，这暗示可以通过一个弯曲的曲线来提高精度。随后，作者添加了二次项（X2），也就是拟合了“Y =β1X2+β2X+β0”的关系，此时拟合曲线能够较好地呈现群落相似度与距离的关系。

6000km 之后为何开始出现群落相似度随距离增加而上升的趋势呢？10000km 之后这种趋势尤为明显。其原因在于此时环境选择过程成为主要因素，而扩散限制变得不那么重要了。作者在文中对这些进行了详细讨论，由于篇幅原因，仅节选了一小部分来帮助理解。

扩散限制的影响较小，一部分原因可解释为海洋的动态极为剧烈。这导致在海洋微生物群落中很难明确界定真正的地方性物种。例如，有研究显示，每天每平方米有数百万种微生物在移动，并且能够在大约 4 天的时间内传播超过 11000 千米。

南北半球存在“对称性”，即由南到北跨越温带 - 热带 - 温带，这是与仅单独观察南半球或北半球时的一大区别，单独观察南半球或北半球时，由南到北仅表现在温带 - 热带或热带 - 温带过渡，这一点不可忽视。考虑环境选择。例如，对于超过 10000km 的距离，在本研究中，这基本上对应了最南部和最北部采样点间的距离。虽然它们地理距离跨度很大，但当考虑所处的纬度时，可以看到它们位于相似的温带地区。相似的温度、季节变化等因素是部分群落组成相似的原因。

接下来，以这种二次关系作为示例，来展示在 R 语言中拟合多项式回归的方式。

带多次方的线性回归在 R 函数 lm()中通过 I()来指定。其中，I(X^2)表示二次方，I(X^3)表示三次方，依此类推。

读取上述输出文件，然后对南北半球大西洋的样点距离数据和群落相似度数据进行综合考虑。
comm_dis <- read.delim('comm_dis.depth40.Both.txt', sep = '\t')
 
#注：site_dis 列的单位是米（m），在这种大尺度范围下，以千米（km）作为度量会更好一些
comm_dis$site_dis_km <- comm_dis$site_dis/1000
 
#lm() 用于拟合群落组成相似度（comm_sim）与地理距离（site_dis_km）之间的线性关系
#首先不妨先看一下简单一次线性回归的结果
fit1 <- lm(comm_sim~site_dis_km, data = comm_dis)
展示拟合模型的简单统计，发现 R2 的值很低。
 
#然后上述简单线性回归基础上添加个二次项
添加二次项到 I(X^2) ，对于 I(X^3) 能够添加三次项，按照这样的规律依次类推。
fit2 <- lm(comm_sim~site_dis_km+I(site_dis_km^2), data = comm_dis)
展示拟合模型的简单统计，它相比一次线性回归能有效提高回归精度。

根据输出结果，可知拟合的一元二次方程式为：“Y 等于 0.00000000673 乘以 X 的平方减去 0.0000986 乘以 X 再加上 0.6066”。因为两组变量的数量级相差太大了，所以会有这么多小数点。

），并且是显著的。R2约为0.3也还是个可观的数值。

这个关系是同时考虑南北半球大西洋时的。

#ggplot2 绘制带线性拟合线的散点图
library(ggplot2)
 
p <- ggplot(comm_dis, aes(x = site_dis_km, y = comm_sim)) +
geom_point(color = 'black') +
使用方法为“lm”，不显示标准误差，颜色为“blue”，公式为 y~x 的 geom_smooth 加上
使用 method 为 'lm'，se 为 FALSE，颜色为 'blue4，公式为 y~poly(x, 2) 进行 geom_smooth 操作，然后再加上其他相关内容（如果有其他相关内容的话）。
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'),
    legend.position = 'none', plot.title = element_text(hjust = 0.5)) +
scale_y_continuous(limit = c(0, 1)) +
在 labs 中，x 为“Distance (km)”，y 为“Bray-curtis similarity”，标题为“Both”。
 
这次不把具体的拟合式标出来，只是把简单一次和二次回归的 R2 进行比较。
label_fit <- data.frame(
R1_2 等于使用 sprintf 函数生成的字符串，该字符串内容为“italic(R^2) == 保留三位小数的 summary(fit1)$adj.r.squared 的值”
R2_2 等于通过 sprintf 函数生成的字符串，该字符串为“italic(R^2) == %.3f”，其中的 %.3f 表示保留三位小数，而这个值是通过 summary(fit2)$adj.r.squared 计算得出的。
 
p +
geom_text 中 x 的值为 10000，y 的值为 0.9，通过 aes 函数指定标签为 R1_2，数据来自 label_fit，并且设置了 parse 为 TRUE
hjust 等于 0，颜色为“蓝色”，不显示图例。
geom_text 的 x 坐标为 10000，y 坐标为 0.8，其标签由 R2_2 决定，数据来自 label_fit，并且开启了 parse 选项。
hjust 等于 0 ，颜色为“blue4”，不显示图例。