线性关系r范围_一个简单线性回归和多项式回归在R中的实现示例

一个简单线性回归和多项式回归在R中的实现示例 常见线性回归的原理就不多说了，大家都懂，就是普通最小二乘法( Ordinary Least Square ， OLS )回归法，主要用于简单线性回归、多项式回归和多元线性回归。

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

本篇就先从简单线性回归和多项式回归开始，对如何在R中执行回归分析作简介。

某文献中的回归示例(部分节选)

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

作者在大西洋表层深度(20-200m)采集水样，采样地点横跨南北半球覆盖了约12000 km的距离，并对浮游细菌进行生物地理学相关研究。在正文result的第3部分，作者通过回归方程(包括简单一次和二次线性回归)，描述大西洋表层浮游细菌群落相似度的距离衰减模式。

例如下图节选自于原文图3B，横轴为采样地间的地理距离，纵轴代表了采样地间40m 水深处FL浮游细菌(文中对一类小型浮游细菌类群的统称)群落的组成相似度。South、North、Both分别代表了单独考虑南半球、单独考虑北半球、综合考虑南北半球时的群落相似度距离-衰减趋势。

备注：红色的点代表一些与边界地区(大西洋中靠近最南和最北地带)有关的样本，作者在分析中将它们特殊对待，下文中提到。

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

示例数据和预处理过程

在原文献附表“Table S3”中，作者提供了研究中采样地点的地理位置(经纬度)以及浮游细菌物种组成等数据，这些可以搜索原文找到下载。

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

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

https://pan.baidu.com/s/1UcPqDNOtNGYFDZk-_SzDuA

节选的部分示例数据

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

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

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

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

(1)根据记录的各地点的经纬度坐标，计算出采样地间的地理距离；

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

##根据经纬度计算采样点的地理距离
#读取采样点的地理位置数据
site #site  0), ]  #只选择北半球
#site  
#计算采样点间的地理距离
#geosphere 包 distm() 根据经纬度计算地理距离(默认距离单位，米)
#distm() 要求两列数据，第一列是经度，第二列是纬度
site_dis rownames(site_dis) colnames(site_dis)  
#将采样点地理距离矩阵转换为两两对应数值的数据框结构
site_dis site_dis head(site_dis)
##计算群落间物种组成相似度，根据原文描述，使用 Bray-curtis 相似度
#浮游细菌群落物种组成数据
spe spe spe  
#vegan 包 vegdist() 计算群落间物种组成 Bray-curtis 相异度矩阵
#并通过 1-Bray-curtis 相异度获得 Bray-curtis 相似度
comm_sim  
#将矩阵转换为两两群落对应数值的数据框结构
diag(comm_sim) comm_sim[upper.tri(comm_sim)] comm_sim comm_sim head(comm_sim)
##采样点距离和群落相似度数据合并
comm_dis names(comm_dis)  
#标识与最北 NADR 或最南 FKLD 地区中采样地点有关的距离对
site_edge comm_dis[which(as.character(comm_dis$site1) %in% site_edge | as.character(comm_dis$site2) %in% site_edge),'edge'] comm_dis[which(! comm_dis$edge %in% '0'),'edge'] head(comm_dis)
write.table(comm_dis, 'comm_dis.depth40.Both.txt', sep = '\t', row.names = FALSE, quote = FALSE)
#write.table(comm_dis, 'comm_dis.depth40.North.txt', sep = '\t', row.names = FALSE, quote = FALSE)
#write.table(comm_dis, 'comm_dis.depth40.South.txt', sep = '\t', row.names = FALSE, quote = FALSE)

根据选择南半球或北半球还是南北半球整体考虑，分别输出文件“comm_dis.depth40.South.txt”、“comm_dis.depth40.North.txt”和“comm_dis.depth40.Both.txt”，内容如上。

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

备注：这里使用R包geosphere计算地理距离，所得结果与原文中计算的地理距离似乎有微小区别(数值不完全相同，可能是在考虑地球椭球形状时存在估计偏差)，因此下文线性回归的拟合R²也不完全一样(其实差别很小，整体趋势是一致的)，所以无需见怪……

一个简单线性回归示例

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

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

在文中，作者将南北半球分开讨论，归因于不同的生物地理学的驱动力。

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

首先考虑北半球的情形。

根据原文中的描述，作者观察到来自最北部地区NADR的采样地的群落结构不存在明显特殊性，它们与其它地点的相似度也不存在明显偏离趋势(图中红色的点)，所以使用文件中所有数据直接拟合。

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

#读取上述输出文件，首先以北半球的为例
comm_dis  
#注：site_dis 列的单位是米(m)，在这种大尺度范围下，以千米(km)作为度量会更好一些
comm_dis$site_dis_km  
#lm() 拟合群落组成相似度(comm_sim)地理距离(site_dis_km)的一元线性关系
fit summary(fit)  #展示拟合模型的简单统计

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

根据输出结果，可知拟合的一元线性方程式为“Y = -0.00007684X + 0.5706”，并且是显著的。R²约为0.4，不算很低，表明北大西洋40m水深处FL浮游细菌群落存在明显的线性距离衰减趋势。

如果期望从中提取关键的信息，或者作一些其它统计，这是一些参考用法。

#例如
coefficients(fit)[1]  #获取截距
coefficients(fit)[2]  #获取 site_dis_km 的斜率
#也可以 names(summary(fit)) 后查看主要的内容项，然后从中提取，例如
summary(fit)$adj.r.squared  #校正后 R2

最后，绘制散点图展示北半球大西洋FL浮游细菌群落的距离-衰减关系。

#ggplot2 绘制带线性拟合线的散点图
#拟合线同时考虑了“黑点”和“红点”，其中“红点”代表了与 NADR 地区中采样地点有关的距离关系
library(ggplot2)
p geom_point(aes(color = as.character(edge))) +
scale_color_manual(values = c('red', 'black'), limit = c('0', '1')) +
geom_smooth(method = '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)) +
labs(x = 'Distance (km)', y = 'Bray-curtis similarity', title = 'North')
#提取上述结果中重要的统计值，用于在图中将方程式标注出
label_fit     formula = sprintf('italic(Y) == %.7f*italic(X) + %.3f', coefficients(fit)[2], 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 = 3000, y = 0.9, aes(label = formula), data = label_fit, parse = TRUE,
    hjust = 0, color = 'black', show.legend = FALSE) +
geom_text(x = 4000, y = 0.8, aes(label = R2), data = label_fit, parse = TRUE,
    hjust = 0, color = 'black', show.legend = FALSE) +
geom_text(x = 4000, y = 0.7, aes(label = p_value), data = label_fit, parse = TRUE,
    hjust = 0, color = 'black', show.legend = FALSE)

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

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

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

至于为什么FKLD中的样本偏离如此明显，可能归因于其地理位置和环境特殊。作者在文中做了些讨论，篇幅原因仅节选一小部分帮助理解。

例如，相比南半球大西洋的其它采样地点，FKLD地区处于多个洋流带的交界处，并具有区别明显的盐度(Salinity)和温度(Temperature)环境，独特的水文特征可能是塑造更别致群落结构的原因。

最后选择将数据中与FKLD有关的距离关系剔除，然后用剩余数据在R函数lm()中以lm(Y~X)的方程式拟合。

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

根据输出结果，可知拟合的一元线性方程式为“Y = -0.0001032X + 0.8098”，并且是显著的。R²约为0.67，很高的数值了，表明南大西洋40m水深处FL浮游细菌群落存在明显的线性距离衰减趋势。

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

#ggplot2 绘制带线性拟合线的散点图
#拟合线只考虑“黑点”，“红点”代表了与 FKLD 地区中采样地点有关的距离关系，因偏离太大被剔除
library(ggplot2)
p geom_point(data = comm_dis, aes(x = site_dis_km, y = comm_sim, color = as.character(edge))) +
scale_color_manual(values = c('red', 'black'), limit = c('0', '1')) +
geom_smooth(data = comm_dis_select, aes(x = site_dis_km, y = comm_sim), method = '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)) +
labs(x = 'Distance (km)', y = 'Bray-curtis similarity', title = 'South')
#提取上述结果中重要的统计值，用于在图中将方程式标注出
label_fit     formula = sprintf('italic(Y) == %.6f*italic(X) + %.3f', coefficients(fit)[2], 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, aes(label = formula), data = label_fit, parse = TRUE,
    hjust = 0, color = 'black', show.legend = FALSE) +
geom_text(x = 5000, y = 0.8, aes(label = R2), data = label_fit, parse = TRUE,
    hjust = 0, color = 'black', show.legend = FALSE) +
geom_text(x = 5000, y = 0.7, aes(label = p_value), data = label_fit, parse = TRUE,
    hjust = 0, color = 'black', show.legend = FALSE)

一个多项式回归示例

所谓多项式，就是对于预测变量是多次方形式，例如“Y~X²”、“Y~X³”等。多项式等式仍认为是线性回归模型，因为等式仍是预测变量的加权和形式。

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

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

尽管上述分析显示南、北半球大西洋中的群落相似度距离衰减都是递减趋势，但当作者综合考虑南北半球时，发现FL浮游细菌群落的组成相似度并非随着采样地间地理距离的增加而一直降低，距离在6000km之后反而会出现上升。

此时，二者间的关系很难再通过“Y = β₁X+β₀”这种形式描述出来，作者也发现了直接拟合简单线性回归的R²是很低的。考虑到整体是先下降后上升的趋势，暗示可以替换为一个弯曲的曲线提高精度，随后作者添加了二次项(X²)，即拟合了“Y = β₁X² + β₂X+β₀”关系，此时拟合曲线较好地呈现群落相似度和距离的关系。

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

关于扩散限制的影响较小，有一部分可解释于海洋的动态非常剧烈，使得海洋微生物群落中很难定义真正的地方性物种。例如有研究表明，每天每m²存在数百万种微生物的移动，并且可以在大约4天的时间内传播超过11000km。

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

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

对于带多次方的线性回归，R函数lm()中通过I()指定，I(X^2)就是二次方，I(X^3)就是三次方，以此类推。

#读取上述输出文件，综合考虑南北半球大西洋的样点距离与群落相似度数据
comm_dis  
#注：site_dis 列的单位是米(m)，在这种大尺度范围下，以千米(km)作为度量会更好一些
comm_dis$site_dis_km  
#lm() 拟合群落组成相似度(comm_sim)地理距离(site_dis_km)的线性关系
#首先不妨先看一下简单一次线性回归的结果
fit1 summary(fit1)  #展示拟合模型的简单统计，发现 R2 是很低的
#然后上述简单线性回归基础上添加个二次项
#I(X^2) 添加二次项，I(X^3) 可添加三次项，以此类推
fit2 summary(fit2)  #展示拟合模型的简单统计，相比一次线性回归有效提高了回归精度

根据输出结果，可知拟合的一元二次方程式为“Y = 0.00000000673X² - 0.0000986X +0.6066”(这么多小数点是因为两组变量的数量级相差太大了)，并且是显著的。R²约为0.3也还是个可观的数值。

绘制散点图展示同时考虑南北半球大西洋时FL浮游细菌群落的距离-衰减关系。

#ggplot2 绘制带线性拟合线的散点图
library(ggplot2)
p geom_point(color = 'black') +
geom_smooth(method = 'lm', se = FALSE, color = 'blue', formula = y~x) +
geom_smooth(method = 'lm', se = FALSE, color = 'blue4', formula = y~poly(x, 2)) +
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', title = 'Both')
#这次就不标出具体的拟合式了，只标出简单一次和二次回归的 R2 比较
label_fit     R1_2 = sprintf('italic(R^2) == %.3f', summary(fit1)$adj.r.squared),
    R2_2 = sprintf('italic(R^2) == %.3f', summary(fit2)$adj.r.squared))
p +
geom_text(x = 10000, y = 0.9, aes(label = R1_2), data = label_fit, parse = TRUE,
    hjust = 0, color = 'blue', show.legend = FALSE) +
geom_text(x = 10000, y = 0.8, aes(label = R2_2), data = label_fit, parse = TRUE,
    hjust = 0, color = 'blue4', show.legend = FALSE)

友情链接

最大信息系数(MIC)稳健量化相关强度及非线性关联网络构建

几种常见的判别分析分类方法

随机森林分类模型以及对重要变量的选择

分段结构方程建模

层次聚合分类分析

基于距离的差异检验方法

LEfSe分析识别biomarker的一般过程

RDA、db-RDA(CAP)及CCA的变差分解(VPA)

基于距离的冗余分析(db-RDA)

主坐标分析(PCoA)和非度量多维标度(NMDS)分析

常见的群落相似性或距离测度的计算

常见Alpha多样性指数计算