ggplot2的回归模型拟合以及回归曲线做图方法杂烩

ggplot2的回归模型拟合以及回归曲线做图方法杂烩

本篇对 ggplot2 绘制回归曲线的方法做个汇总。

在先前的不少回归分析相关的推文中，大多都在过程当中说起了一些简单的可视化方法，以帮助理解回归中的变量响应关系。在这些做图方法中，有些是特定功能 R 包自带的可视化方法，有些经过 ggplot2 绘制，还有部分直接使用 base 做图。相对而言， ggplot2 是 R 里面很是专业的做图包，其同时涵括了多种类型的回归模型在内，这使回归曲线的实现更加轻松。

您是否常用 geom_smooth() 拟合线性回归或 LOESS 平滑？除了这两种模型外，您是否还了解过 ggplot2 实现其它回归类型的方法，或者其它方式的回归曲线绘制函数呢？下文就列举一些示例帮助您加深印象吧。尽管如此，并非全部回归均可以绘制图形展现。例如在二维平面或者三维空间中，可表达的信息有限，多元回归（具备多维的变量信息）就难以呈现出来。所以，本篇就以最多见的双变量关系分析（一元回归）为例，展现回归曲线经过 ggplot2 的绘制方法。

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

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

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

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

经过ggplot2函数stat_smooth()、geom_smooth()统计变量关系并拟合回归机器学习

ggplot2内置了回归统计函数，想必同窗们也都很熟悉，就是stat_smooth()以及geom_smooth()。两个函数的使用方法大体一致，可基于给定的双变量关系拟合指定的回归，默认方法包含局部加权回归（loess）、通常线性模型（lm）、广义线性模型（glm）、广义加性模型（gam）以及一些非线性回归（nls，如指数回归）等，并同时将双变量关系的回归线绘制在二维图中表示出来。此外，结合一些拓展包也能够实现更多非线性回归的做图。编辑器

接下来展现几个简单示例，了解不一样模式回归曲线经过ggplot2函数的绘制方法。函数

示例1.1-局部加权回归（loess）的平滑拟合学习

MASS包的内置数据集Boston，记录了波士顿郊区的房屋价值与当地社会环境、人口组成等信息。这里但愿关注当地低收入人口比例和其所居住地房屋价格是否有关，不妨首先拟合一个双变量关系的平滑曲线进行探索。url

平滑回归方法中，局部加权回归（LOESS）是最经常使用的平滑器。在函数stat_smooth()或geom_smooth()中经过method = 'loess'指定局部加权回归的平滑拟合，并经过span参数控制曲线平滑度。spa

#MASS 包的示例数据集，波士顿郊区的房屋价值，详情 ?Boston
library(MASS)
data(Boston)
 
#ggplot2 绘制低收入人口比例和居住地房屋价格的关系
library(ggplot2)
 
p <- ggplot(data = Boston, aes(x = medv, y = lstat)) +
geom_point() +
theme_bw()
 
p
 
##拟合局部加权回归（loess），并展现 95% 置信区间
#stat_smooth() 参数 span 可控制曲线平滑度，值越大拟合曲线越平滑，详情 ?stat_smooth
#注意平滑参数值应当谨慎设置，它强烈影响曲线意义的解读
p1 <- p + stat_smooth(method = 'loess', span = 0.5, se = TRUE, level = 0.95) 
p1
 
#相似地，geom_smooth() 参数中，span 控制曲线平滑度
p2 <- p + geom_smooth(method = 'loess', span = 0.75, se = TRUE)
p2

平滑曲线展示了居住地房屋价格越高，当地低收入人口比例越低。可能与低收入人群难以负担得起较高居住地的成本，而更多选择廉价的房屋租赁或长期居住有关。

示例1.2-通常加性模型的平滑拟合

实现非参数平滑回归的另外一方法是使用加性模型。在函数stat_smooth()或geom_smooth()中经过method = 'gam'指定加性模型的平滑拟合，并经过在formula参数中选择平滑器参数实现。默认状况下，该函数假定响应变量服从正态分布，并执行通常加性模型的形式。

##gam 默认状况下执行通常加性模型的平滑回归，并展现 95% 置信区间
#stat_smooth() 和 geom_smooth() 的用法大体类似
#经过 formula = y~s(x) 指定函数关系，s() 中的参数 k 可调整平滑度，详情 ?s()
#注意平滑参数值应当谨慎设置，它强烈影响曲线意义的解读
p + stat_smooth(method = 'gam', formula = y~s(x, k = 10), se = TRUE, level = 0.95) 
p + geom_smooth(method = 'gam', formula = y~s(x, k = 5), se = TRUE)

和loess平滑显示了一致的趋势，平滑曲线展示了居住地房屋价格越高，当地低收入人口比例越低。

示例1.3-通常线性模型（lm）之简单线性回归

上图的非参数平滑曲线显示了低收入人口的居住比例随当地房屋价格的升高而下降，不妨使用参数回归代替描述这个趋势，例如最常使用的简单线性回归。函数stat_smooth()或geom_smooth()中，经过method = 'lm'指定线性回归拟合变量关系，默认绘制Y=β₁X+β₀形式的简单线性回归形式。

##通常线性模型（lm），并展现 95% 置信区间
#stat_smooth() 拟合简单线性回归，详情 ?stat_smooth
p1 <- p1 + stat_smooth(method = 'lm', se = TRUE, level = 0.95, color = 'red')
p1
 
#geom_smooth() 拟合简单线性回归，详情 ?geom_smooth
p2 <- p2 + geom_smooth(method = 'lm', se = TRUE, color = 'red')
p2

该图显示了波士顿郊区低收入人口的居住比例随当地房屋价格的升高而呈现明显的线性下降趋势。蓝线，上文的非参数loess平滑拟合；红线，普通线性回归拟合。

示例1.4-通常线性模型（lm）之带屡次项的线性回归

不过上述非参数平滑曲线图的走势代表，彷佛还能够用一个弯曲的曲线来提升拟合精度，即拟合含屡次项的线性回归。所以继续考虑添加二次项关系，即Y=β₁X+β₂X²+β₀。相似地，在函数stat_smooth()或geom_smooth()中经过method = 'lm'指定线性回归，并经过formula指定变量间屡次项关系。

#添加低收入人口比例和居住地房屋价格的二次线性回归，并展现 95% 置信区间
#geom_smooth() 默认 formula=y~x 对应一次关系，可更改 formula = y~poly(x, n) 拟合 n 次关系，详情 ?stat_smooth
p1 + stat_smooth(method = 'lm', formula = y~poly(x, 2), se = TRUE, level = 0.95, color = 'green3')
 
#一样地，geom_smooth() 默认 formula=y~x 对应一次关系，可更改 formula = y~poly(x, n) 拟合 n 次关系，详情 ?geom_smooth
p2 + geom_smooth(method = 'lm', formula = y~poly(x, 2), se = TRUE, color = 'green3')

能够看到，考虑了低收入人口比例和居住地房屋价格的二次线性关系后，拟合精度更优，更能描述出两者实际的趋势。蓝线，上文的非参数loess平滑拟合；红线，上文普通线性回归拟合；绿线，二次线性回归拟合。

示例1.5-广义加性模型（gam）之泊松加性模型

网盘示例数据“fish_data.txt”，来自前文泊松响应的广义加性模型的示例，记录了河流中某鱼类个体数量与水体环境的关系。前文已经经过非参数的泊松加性模型探索了影响鱼类物种丰度的环境因素，并描述了响应曲线的特征，在这里将继续经过ggplot2绘制泊松加性模型的平滑拟合曲线。

如上文所述，在函数stat_smooth()或geom_smooth()中经过method = 'gam'指定加性模型的平滑拟合。因为函数默认响应变量服从正态分布并执行通常加性模型的形式，所以在这里须要额外在参数method.args中指定响应变量的泊松分布，以执行泊松加性模型。此外，还需在formula参数中设定平滑器参数。

#读取示例的鱼类物种丰度和水体环境数据
dat <- read.delim('fish_data.txt', sep = '\t', row.names = 1)
 
#ggplot2 绘制鱼类物种丰度与水域流域面积的关系
library(ggplot2)
 
p <- ggplot(data = dat, aes(x = acre, y = fish)) +
geom_point() +
theme_bw()
 
p
 
##计数型变量，尝试使用泊松加性模型拟合平滑曲线探索变量关系，并展现 95% 置信区间
#stat_smooth() 和 geom_smooth() 的用法大体类似
#经过 method.args 指定响应变量的类型，这里使用泊松加性模型
#经过 formula = y~s(x) 指定函数关系，平滑器 s() 中的参数 k 可调整平滑度，详情 ?s()
#注意平滑参数值应当谨慎设置，它强烈影响曲线意义的解读
p1 <- p + stat_smooth(method = 'gam', formula = y~s(x, k = 5), method.args = list(family = 'poisson'), se = TRUE, level = 0.95)
p1
 
p2 <- p + geom_smooth(method = 'gam', formula = y~s(x, k = 3), method.args = list(family = 'poisson'), se = TRUE)
p2

整体来讲，随水域流域面积（acre）的增长，鱼类物种数量将提高。

示例1.6-广义线性模型（glm）之泊松回归

如上非参数的泊松加性模型显示了鱼类物种数量随水域流域面积的增长而提高，这种提高效应彷佛能够经过线性关系来描述。所以，接下来不妨使用参数化的广义线性模型描述这个趋势，即泊松回归。

不管使用函数stat_smooth()或geom_smooth()，均经过method = 'glm'指定使用广义线性回归拟合变量关系，并需额外在参数method.args中指定响应变量的泊松分布。

##计数型变量，尝试使用广义线性模型的泊松回归描述变量关系，并展现 95% 置信区间
#stat_smooth() 和 geom_smooth() 中，经过 method.args 指定响应变量的类型，这里使用泊松回归
p1 + stat_smooth(method = 'glm', method.args = list(family = 'poisson'), se = TRUE, level = 0.95, color = 'red') 
p2 + geom_smooth(method = 'glm', method.args = list(family = 'poisson'), se = TRUE, color = 'red')

泊松回归有效体现出了鱼类物种数量随水域流域面积的增长而提高的线性响应趋势。蓝线，上文非参数泊松加性模型的平滑拟合；红线，泊松回归的线性拟合。

示例1.7-广义线性模型（glm）之logistic回归型

考虑到广义线性模型的种类很是多，这里再展现一例广义线性模型的其它常见类型，logistic回归在ggplot2中的绘制方法，加深对ggplot2直接绘制广义线性模型拟合曲线方法的理解。

网盘示例数据“breast.csv”，来自前文广义线性模型之logistic回归的示例，记录了乳腺癌恶性肿瘤相关的细胞特征得分。前文已经经过logistic回归肯定了多种细胞特征得分与肿瘤恶性状态（良性或恶性，属于二分响应变量）的关联，并评估了经过这些细胞特征得分区分肿瘤良性或恶性的准确度。在这里将继续经过ggplot2绘制logistic回归的拟合曲线。

相似地，在函数stat_smooth()或geom_smooth()中经过method = 'glm'指定使用广义线性回归拟合变量关系，并额外在参数method.args中指定响应变量的二项分布。

#威斯康星州乳腺癌数据集
breast <- read.csv('breast.csv')
 
#使用 0-1 的二分数值从新定义肿瘤的两种状态（0，良性；1，恶性）
breast[which(breast$class == 2),'response'] <- 0
breast[which(breast$class == 4),'response'] <- 1
 
#例如以肿瘤厚度得分与肿瘤良性或恶性状态的关系为例
library(ggplot2)
 
p <- ggplot(data = breast, aes(x = clumpThickness, y = response)) +
geom_jitter(height = 0.05, width = 0, alpha = 0.5) +  #做图时对点添加了少许上下抖动效果仅为减小点的重叠
theme_bw()
 
p
 
##二分响应变量，使用广义线性模型的二项 logistic 回归描述变量关系，并展现 95% 置信区间
p + stat_smooth(method = 'glm', method.args = list(family = 'binomial'), se = TRUE, level = 0.95, color = 'red') 
p + geom_smooth(method = 'glm', method.args = list(family = 'binomial'), se = TRUE, color = 'red')

二项响应的logistic回归经常使用于描述响应变量被分配为两个类别（0或1）的几率。在这个实例中，0表明肿瘤良性，1表明肿瘤恶性。图中红色曲线为logistic回归得到的几率拟合曲线，形象呈现了肿瘤为恶性的几率随肿瘤厚度得分增长的非线性响应。

示例1.8-非线性参数回归示例之指数回归

再展现一例常见的参数非线性回归方法，指数回归在ggplot2中的绘制。

仍以上文提到的波士顿郊区的房屋价值数据集为例，经过局部加权回归（LOESS）的平滑曲线观察居住地房屋价格和人均城镇犯罪率的关系时，发现两者彷佛存在指数关联。随人均城镇犯罪率的升高，房屋价格开始呈现剧烈的降低，随后达到平缓。所以，考虑在图中拟合指数回归曲线。

经过在stat_smooth()或geom_smooth()中指定method = 'nls'便可拟合参数非线性回归，formula中设置非线性回归式，method.args中设定一些必要的补充信息。

#MASS 包的示例数据集，波士顿郊区的房屋价值，详情 ?Boston
library(MASS)
data(Boston)
 
#ggplot2 绘制居住地房屋价格和人均城镇犯罪率的关系
library(ggplot2)
 
p <- ggplot(data = Boston, aes(x = crim, y = medv)) +
geom_point() +
theme_bw() +
stat_smooth(method = 'loess', span = 0.5, se = TRUE, level = 0.95) 
 
p
 
##参数非线性回归（nls）的拟合，以指数回归为例
#stat_smooth() 和 geom_smooth() 的用法大体类似
#参数中 a 和 b 的初始值手动指定，可参考 https://mp.weixin.qq.com/s/fxkMlEBs8UEwuGlUi50uLA
p + stat_smooth(method = 'nls', formula = y ~ a*x^b, method.args = list(start = list(a = 2, b = 1.5)), se = FALSE, color = 'red')
p + geom_smooth(method = 'nls', formula = y ~ a*x^b, method.args = list(start = list(a = 2, b = 1.5)), se = FALSE, color = 'red')

这个指数回归大体描述出了随人均城镇犯罪率的升高，房屋价格首先呈现剧烈的降低，并在随后达到平缓的趋势。蓝线，上文的非参数loess平滑拟合；红线，指数回归拟合曲线。

经过ggplot2函数stat_function()指定回归函数绘制拟合线

不少回归方法，特别对于大多数非线性回归而言，ggplot2及其拓展包中缺乏做图方案，难以经过stat_smooth()、geom_smooth()直接做图。这时候，能够考虑使用stat_function()根据指定函数绘制拟合线。

若是已经提早计算出了回归式的各参数，则能够直接将已知的回归式指定给ggplot2函数stat_function()。stat_function()可以在做图时将自变量代入至已知的回归式中拟合响应变量的预测值，并使用平滑线链接响应变量的预测值得到回归线。理论上，stat_function()能够实现对任意给定回归式的拟合线绘制。

示例2.1-一个带屡次项的线性回归示例

例如在上文“示例1.4-通常线性模型（lm）之带屡次项的线性回归”中，已经描述了低收入人口比例和居住地房屋价格的二次线性关系。如要经过stat_function()绘制拟合曲线，则操做过程以下。

#MASS 包的示例数据集，波士顿郊区的房屋价值，详情 ?Boston
library(MASS)
data(Boston)
 
#ggplot2 绘制低收入人口比例和居住地房屋价格的关系
library(ggplot2)
 
p <- ggplot(data = Boston, aes(x = medv, y = lstat)) +
geom_point() +
theme_bw()
 
p
 
##经过 lm() 计算低收入人口比例和居住地房屋价格的二次线性关系
fit <- lm(lstat~medv+I(medv^2), data = Boston)
summary(fit)
 
#根据 lm() 得到的对回归系数（斜率）和截距的参数估计值
#可知二次回归式书写为 Y = -1.715X + 0.021X^2 + 39.048
#而后将这个回归式书写到 stat_function() 中
p + stat_function(fun = function(x) -1.715*x+0.021*x^2+39.048, color = 'red', size = 1.2)
 
#或者直接在 lm() 结果中提取回归系数（斜率）和截距
#而后将提取出的值代入至 stat_function() 中
b0 <- coefficients(fit)[1]  #截距
b1 <- coefficients(fit)[2]  #一次项的回归系数
b2 <- coefficients(fit)[3]  #二次项的回归系数
 
p + stat_function(fun = function(x) b1*x+b2*x^2+b0, color = 'red', size = 1.2)

如此方法，可将已知的低收入人口比例和居住地房屋价格的二次线性关系添加在图中。

示例2.2-一个指数回归的非线性参数回归示例

相似地，再以上文“示例1.8-非线性参数回归示例之指数回归”的结果为例，展现一个非线性参数回归经过stat_function()绘制拟合曲线的过程。

#ggplot2 绘制居住地房屋价格和人均城镇犯罪率的关系
library(ggplot2)
 
p <- ggplot(data = Boston, aes(x = crim, y = medv)) +
geom_point() +
theme_bw()
 
p
 
##参数非线性回归（nls）的拟合，以指数回归为例
#参数中 a 和 b 的初始值手动指定，可参考 https://mp.weixin.qq.com/s/fxkMlEBs8UEwuGlUi50uLA
mod <- nls(medv ~ a*crim^b, data = Boston, start = list(a = 2, b = 1.5))
summary(mod)
 
#根据 nls() 得到的对 a 和 b 参数估计值
#可知指数回归式书写为 Y = 20.743X^-0.086
#而后将这个回归式书写到 stat_function() 中
p + stat_function(fun = function(x) 20.743*x^-0.086, color = 'red', size = 1.2)
 
#或者直接在 nls() 结果中提取 a 和 b 的值
#而后将提取出的值代入至 stat_function() 中
a <- coefficients(mod)[1]  #a 的值
b <- coefficients(mod)[2]  #b 的值
 
p + stat_function(fun = function(x) a*x^b, color = 'red', size = 1.2)

如此方法，可将已知的房屋价格随人均城镇犯罪率的升高而下降的大体的指数响应趋势添加在图中。

根据已知回归模型预测指望值，并在图中绘制后链接成线

除了经过在stat_function()中手动输入回归关系式外，也能够经过已经在R函数中构建好的回归模型以predict()预测响应变量值，并在图中绘制后以线段或曲线的形式链接起来获得回归线。这种方法更多应用于回归式复杂难以书写的状况，或者应用在其它非参数回归、机器学习方法等难以理解的形式中。

示例3.1-一个非参数的平滑回归示例

例如上文“示例1.1-局部加权回归（loess）的平滑拟合”的平滑回归曲线，是非参数形式的，难以经过具体数学公式来描述。除了直接经过stat_smooth()或geom_smooth()实现外，但若在stat_function()中输入公式则确定行不通。虽然咱们也不必展现其它冗余的方法了，但此处仅仅是想举例告知其它的可选替代形式，所以请容许我继续。

#MASS 包的示例数据集，波士顿郊区的房屋价值，详情 ?Boston
library(MASS)
data(Boston)
 
#ggplot2 绘制低收入人口比例和居住地房屋价格的关系
library(ggplot2)
 
p <- ggplot(data = Boston, aes(x = medv, y = lstat)) +
geom_point() +
theme_bw()
 
p
 
##经过 R 函数 loess() 构建局部加权回归（LOESS）模型
#参数 span 可控制曲线平滑度，值越大拟合曲线越平滑，详情 ?loess
#注意平滑参数值应当谨慎设置，它强烈影响曲线意义的解读
fit <- loess(lstat~medv, data = Boston, span = 0.75)
summary(fit)
 
#经过 predict() 根据已构建好的回归模型预测响应变量的指望值
#即根据居住地房屋价格得到低收入人口比例的预测值
lstat_pred <- predict(fit, Boston, type = 'response')
head(lstat_pred)
 
#将低收入人口比例的预测值添加在已有的关系图中，并用折线图链接成“拟合线”
Boston$lstat_pred <- lstat_pred
 
p + geom_line(data = Boston, aes(x = medv, y = lstat_pred), color = 'red') 
 
##和 ggplot2 的默认方法 stat_smooth() 的比较
p + stat_smooth(method = 'loess', span = 0.75, se = TRUE, level = 0.95)

局部加权回归（LOESS）平滑曲线展示了居住地房屋价格越高，当地低收入人口比例越低的趋势。红线，首先构建LOESS模型，并经过已知回归模型预测lstat的值后，绘制在图中链接成线；蓝线，ggplot2默认的stat_smooth()方法所得LOESS曲线；两者是一致的。

示例3.2-一个logistic回归的广义线性模型示例

与通常线性模型相比，广义线性模型涉及到链接函数的转换，使其相对复杂不易理解。所以，某些广义线性模型的函数式可能书写麻烦，不易于经过stat_function()描绘出来。除了直接经过stat_smooth()或geom_smooth()实现做图外，也能够用这种替代形式表现出。以上文“示例1.7-广义线性模型（glm）之logistic回归”的数据为例，给个参考过程。

#威斯康星州乳腺癌数据集
breast <- read.csv('breast.csv')
 
#使用 0-1 的二分数值从新定义肿瘤的两种状态（0，良性；1，恶性）
breast[which(breast$class == 2),'response'] <- 0
breast[which(breast$class == 4),'response'] <- 1
 
#例如以肿瘤厚度得分与肿瘤良性或恶性状态的关系为例
library(ggplot2)
 
p <- ggplot(data = breast, aes(x = clumpThickness, y = response)) +
geom_jitter(height = 0.05, width = 0, alpha = 0.5) +  #做图时对点添加了少许上下抖动效果仅为减小点的重叠
theme_bw()
 
p
 
##二分响应变量，使用广义线性模型的二项 logistic 回归描述变量关系
fit <- glm(response~clumpThickness, data = breast, family = binomial())
summary.glm(fit)
 
#经过 predict() 根据已构建好的回归模型预测响应变量的指望值
#即根据肿瘤厚度得分预测肿瘤为恶性的几率
response_pred <- predict(fit, breast, type = 'response')
head(response_pred)
 
#将预测的肿瘤为恶性的几率添加在已有的关系图中，并用折线图链接成“拟合线”
breast$response_pred <- response_pred
 
p + geom_line(data = breast, aes(x = clumpThickness, y = response_pred), color = 'red') 
 
##和 ggplot2 的默认方法 stat_smooth() 的比较
p + stat_smooth(method = 'glm', method.args = list(family = 'binomial'), se = TRUE, level = 0.95)

在这个实例中，0表明肿瘤良性，1表明肿瘤恶性。logistic回归得到的几率拟合曲线，形象呈现了肿瘤为恶性的几率随肿瘤厚度得分增长的非线性响应。红线，首先构建logistic模型，并经过已知回归模型预测肿瘤为恶性的几率后，绘制在图中链接成线；蓝线，ggplot2默认的stat_smooth()方法所得logistic曲线；两者是一致的。

其它拓展做图

上文主要在几种不一样类型的拟合曲线的绘制过程上做了演示，包含局部加权回归（loess）、通常线性模型（lm）、广义线性模型（glm）、广义加性模型（gam）以及一些非线性回归（nls，如指数回归）等，帮助了解不一样模式回归曲线怎样经过ggplot2函数的实现方法。

而对于ggplot2主题的调整，怎样让图形更美观，就没做啥说起，你们本身看着来吧。最后再展现几种实用的拓展用法。

示例4.1-一图展现多组的回归

有时涉及了多分组的状况，自变量的水平一致但各组响应变量的值不一样。为了方便多组之间可以直观的比较，指望在同一张图中展现多组的回归曲线，并经过斜率的变化描述趋势的升高或下降幅度，这也是平时很是常见的情形。

网盘示例数据“plant.txt”，测量了生长于两个不一样温度（低、高）环境中的某种植物的高度，以及土壤中的氮浓度，指望了解土壤氮浓度及温度对植物生长的影响。考虑到线性关系最为直观且易于理解，这里但愿经过创建简单线性回归，探讨植物生长高度与其所处环境中土壤氮浓度及温度的关系。

#植物生长与环境条件的数据集
dat <- read.delim('plant.txt', sep = '\t')
 
#横坐标土壤氮浓度，纵坐标植物生长高度
#按高温和低维条件分别标记分组后，拟合线性回归
library(ggplot2)
 
p <- ggplot(data = dat, aes(x = N, y = plant)) +
geom_point(aes(color = temperature)) +
scale_color_manual(values = c('red', 'blue'), limits = c('High', 'Low')) +
geom_smooth(aes(color = temperature), method = 'lm', se = TRUE, show.legend = FALSE) +
theme_bw() +
theme(legend.key = element_blank()) +
labs(x = 'N', y = 'plant height')
 
p

根据做图结果得知，土壤氮含量高时更有利于该植物生长，而且该植物为喜温植物，在温度容许的范围内，温暖更刺激其生长。

示例4.2-在图中添加文字描述

例如，将上述两组的线性回归的公式、R²和p值额外添加在图中，以告知读者回归的精度和显著性水平。此时，您须要额外经过回归函数获取参数估计，并将它们提取出来后，组合成回归公式添加到图中的合适位置。

#经过 lm() 得到线性回归的参数估计
fit_height <- lm(plant~N, data = subset(dat, temperature == 'High'))
summary(fit_height)
 
fit_low <- lm(plant~N, data = subset(dat, temperature == 'Low'))
summary(fit_low)
 
#提取两个回归的回归系数（斜率）、截距、R二、p 值等
b0_height <- coefficients(fit_height)[1]  #截距
b1_height <- coefficients(fit_height)[2]  #斜率
R2_height <- summary(fit_height)$adj.r.squared  #推荐使用校订后的 R2
p_height <- summary(fit_height)$coefficients[2,4]  #p 值
 
b0_low <- coefficients(fit_low)[1]  #截距
b1_lowt <- coefficients(fit_low)[2]  #斜率
R2_low <- summary(fit_low)$adj.r.squared  #推荐使用校订后的 R2
p_low <- summary(fit_low)$coefficients[2,4]  #p 值
 
#将提取出的结果组合成新的数据集，以便添加在原图中
label_height <- data.frame(
    formula_height = sprintf('italic(Y) == %.3f*italic(X) + %.3f', b1_height, b0_height),
    R2_height = sprintf('italic(R^2) == %.3f', R2_height),
    p_height = ifelse(p_height < 0.001, sprintf('italic(P) < 0.001'), ifelse(p_height < 0.01, sprintf('italic(P) < 0.01'), ifelse(p_height < 0.05, sprintf('italic(P) < 0.05'), sprintf('italic(P) >= 0.05')))),
    formula_low = sprintf('italic(Y) == %.3f*italic(X) + %.3f', b1_lowt, b0_low),
    R2_low = sprintf('italic(R^2) == %.3f', R2_low),
    p_low = ifelse(p_low < 0.001, sprintf('italic(P) < 0.001'), ifelse(p_low < 0.01, sprintf('italic(P) < 0.01'), ifelse(p_low < 0.05, sprintf('italic(P) < 0.05'), sprintf('italic(P) >= 0.05')))))
 
#添加文字时手动调整下位置
p +
geom_text(x = 1, y = 42, aes(label = formula_height), data = label_height, parse = TRUE, size = 2.5, hjust = 0, color = 'red', show.legend = FALSE) +
geom_text(x = 1, y = 38, aes(label = R2_height), data = label_height, parse = TRUE, size = 2.5, hjust = 0, color = 'red', show.legend = FALSE) +
geom_text(x = 1, y = 34, aes(label = p_height), data = label_height, parse = TRUE, size = 2.5, hjust = 0, color = 'red', show.legend = FALSE) +
geom_text(x = 1, y = 28, aes(label = formula_low), data = label_height, parse = TRUE, size = 2.5, hjust = 0, color = 'blue', show.legend = FALSE) +
geom_text(x = 1, y = 24, aes(label = R2_low), data = label_height, parse = TRUE, size = 2.5, hjust = 0, color = 'blue', show.legend = FALSE) +
geom_text(x = 1, y = 20, aes(label = p_low), data = label_height, parse = TRUE, size = 2.5, hjust = 0, color = 'blue', show.legend = FALSE)

此时将两组线性回归的公式、R²和p值额外添加在图中，形象地告知了回归统计的可靠性等信息。

示例4.3-分面回归图

经常使用于多个自变量，或者多个响应变量的状况下，双变量回归分析的展现。

例如，这里尝试经过多个简单线性回归，分别描述网盘示例数据“fish_data.txt”中全部6种环境因素对鱼类物种丰度的影响。

#读取鱼类物种丰度和水体环境数据
dat <- read.delim('fish_data.txt', sep = '\t', row.names = 1)
 
#绘制二维散点图观测各环境变量与鱼类物种丰度的关系
#并同时添加简单线性回归拟合线
library(ggplot2)
 
dat_plot <- reshape2::melt(dat, id = 'fish')
 
p <- ggplot(dat_plot, aes(value, fish)) +
geom_point() +
facet_wrap(~variable, ncol = 3, scale = 'free') +
geom_smooth(method = 'lm') +
labs(x = '', y = 'fish')
 
p
 
#如下在散点图中添加各环境变量与鱼类物种丰度的一元回归的 R2 和 p 值做为标识
#注：该 R2 和 p 值仅为多个一元回归的 R2 和 p 值，不要将其视为一个多元回归
 
#拟合各环境变量与鱼类物种丰度的一元回归，并提取各自 R2 和 p 值
env <- c('acre', 'do2', 'depth', 'no3', 'so4', 'temp')
R2_adj <- c()
p_value <- c()
 
for (i in env) {
    fit_stat <- summary(lm(dat[['fish']]~dat[[i]]))  #一元线性回归
    R2_adj <- c(R2_adj, fit_stat$adj.r.squared)  #提取校订后 R2
    p_value <- c(p_value, fit_stat$coefficients[2,4])  #提取显著性 p 值
}
 
env_stat <- data.frame(row.names = env, R2_adj, p_value)
env_stat  #数据框中存储了各环境变量与鱼类物种丰度的一元回归的 R2 和 p 值
 
#自动定义文字在图中的展现坐标位置
env_stat$fish <- (max(dat$fish)-min(dat$fish))*0.85 + min(dat$fish)
for (i in env) env_stat[i,'value'] <- (max(dat[[i]])-min(dat[[i]]))*0.7 + min(dat[[i]])
 
#生成文字标签数据集
env_stat$variable <- rownames(env_stat)
env_stat$R2_adj <- paste(sprintf('italic(R^2) == %.3f', env_stat$R2_adj))
env_stat$p_value <- paste(sprintf('italic(P) == %.3f', env_stat$p_value))
 
#和先前的数据框合并，便于做图
env_stat <- env_stat[c('fish', 'variable', 'value', 'R2_adj', 'p_value')]
dat_plot$R2_adj <- NA
dat_plot$p_value <- NA
dat_plot <- rbind(dat_plot, env_stat)
 
p +
geom_text(data = dat_plot, aes(label = R2_adj), parse = TRUE, size = 3, hjust = 0) +
geom_text(data = dat_plot, aes(label = p_value), parse = TRUE, size = 3, hjust = 0, vjust = 3)