用 ggplot2 拟合回归模型和绘制回归曲线的多种方法
下文中所使用的示例数据和R代码的百度盘链接（提取码，egpa）：
https://pan.baidu.com/s/1GwW8EYgFE7EPrW6o9u5ScA
若百度盘失效，也可在GitHub的备份中获取：
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，记录了波士顿郊区的房屋价值与当地社会环境、人口组成等信息。这里希望关注当地低收入人口比例和其所居住地房屋价格是否有关，不妨首先拟合一个双变量关系的平滑曲线进行探索。
平滑回归方法中，局部加权回归（LOESS）是最常用的平滑器。在函数stat_smooth()或geom_smooth()中通过method = 'loess'指定局部加权回归的平滑拟合，并通过span参数控制曲线平滑度。
#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=β1X+β0形式的简单线性回归形式。
##一般线性模型（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=β1X+β2X2+β0。类似地，在函数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