前言
本次学习旨在为了解决某些复杂散点的拟合问题，尝试找到一种回归模型并进行可视化
学习内容提炼自STHDA(http://www.sthda.com/english/)
背景
回归分析（或者回归模型）由一系列机器学习方法组成，允许用户可以通过一个或者多个变量x来预测连续结果变量y，简要的来说，回归模型的目的通过建立一个数学方程式，将y定义为变量x的函数，随后，该公式会被用来基于变量x的新值来预测结果y。
一般来说，线性回归(linear-regression)是最简单也最流行的一种回归模型，用以预测一个连续变量，它假设预测变量和结果之间存在着线性关系。
在有些情况下，预测变量和结果之间的关系并不是线性的，这个时候你就需要建立一个非线性的回归模型(non-linear regression)，例如多项式回归和样条回归(polynomial and spline regression)
你当然可以将多个模型应用在你的数据上，随后进行比较后选择一个最好的模型去解释你的数据，为了进行比较，我们需要一些统计指标来比较不同模型在解释现有数据和结果预测的表现差异；最好的模型应该体现在在最低的预测误差上，用来评价模型的指标常见的有：
- Root Mean Squared Error，用以衡量模型预测误差。它对应着已知观察值与模型预测值之间的平均差异，RMSE的计算公式为
mean(observeds - predicteds)^2 %>% sqrt()，RMSE的值越低，模型也就越适合。 - Adjusted R-square，说到adjusted R方，咱们就不得不说一下R方(R square)，R方用于单变量拟合，在多变量拟合时，就需要对R方做一个（正向或者负向的）惩罚，那就是Adjusted R-square，简单来说R方代表着模型可以解释的方差占原始数据方差的比例，比例越大也就证明模型的预测性越好，至于计算公式咱们容后再表
- 以上两个指标可以简单地通过
caret包的RMSE()和R2()函数来完成计算
需要注意的是，以上提到的统计指标应该用一批新的数据来进行计算（不是用来建立模型的那一批），常见的做法是把一批数据的80%用来作为训练集(training set)，剩下的20%用来作为测试集(test set)，一种粗暴且流行的模型表现衡量方法叫做k-fold cross-validation，它在小数据集的情况下的表现也是不错的，其原理如下
- 随机将数据分为k群(k-fold)，例如5群
- 保留一个子集作为测试集，并使用剩下的k-1个集合合并后进行模型建立
- 用保留子集对训练得到的模型进行评价
- 重复前三步流程k-1次（也就是每个子集都有一次成为测试集的机会）
- 计算k次计算得到的误差的均值，就被成为cross-validation error，用以评价模型表现
使用线性回归举例
本次演示涉及到两个包的使用caret和tidyverse：
# install.packages("caret")
# install.packages("tidyverse")
library(tidyverse)
library(caret)
theme_set(theme_bw())
示例数据来自于GitHub的datarium包，我们使用caret包的createDataPartition函数来创建训练集和测试集：
# Load the data
data("marketing", package = "datarium")
# Inspect the data
sample_n(marketing, 3)
# youtube facebook newspaper sales
# 1 219.12 55.44 70.44 25.44
# 2 332.28 58.68 50.16 32.40
# 3 176.76 28.68 22.92 17.52
# Split the data into training and test set
set.seed(123)
training.samples <- marketing$sales %>%
createDataPartition(p = 0.8, list = FALSE)
train.data <- marketing[training.samples, ]
test.data <- marketing[-training.samples, ]
随后我们使用R语言自带的lm()函数来对训练集进行拟合，使用predict函数来对测试集做出预测，随后将预测数据和实际数据进行比较得到模型的表现：
# Build the model
model <- lm(sales ~., data = train.data)
# Summarize the model
summary(model)
# Call:
# lm(formula = sales ~ ., data = train.data)
#
# Residuals:
# Min 1Q Median 3Q Max
# -10.7142 -0.9939 0.3684 1.4494 3.3619
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 3.594142 0.420815 8.541 1.05e-14
# youtube 0.044636 0.001552 28.758 < 2e-16
# facebook 0.188823 0.009529 19.816 < 2e-16
# newspaper 0.002840 0.006442 0.441 0.66
#
# Residual standard error: 2.043 on 158 degrees of freedom
# Multiple R-squared: 0.8955, Adjusted R-squared: 0.8935
# F-statistic: 451.2 on 3 and 158 DF, p-value: < 2.2e-16
# Make predictions
predictions <- model %>% predict(test.data)
predictions
# 2 3 6 14 15 24 32 40 41 46 51 58
# 15.036290 15.151964 15.395866 10.563146 22.137863 19.741211 13.715581 24.457935 19.601260 18.178608 15.116375 15.296489
# 65 66 68 81 86 91 93 104 107 116 124 129
# 20.412706 9.400329 14.375762 13.812250 18.335633 11.929636 23.046574 17.616955 7.526894 15.726891 18.069954 26.507800
# 134 137 138 139 148 156 168 178 181 182 183 184
# 23.111689 13.833975 25.006188 11.835838 27.874451 6.461594 15.915354 14.597916 12.599533 16.614622 7.997147 28.986837
# 189 192
# 22.075364 10.105750
我们可以看到，使用summary()函数来查看model的时候，事实上也是有三个评价指标的：
- Residual standard error：可以理解为模型预测结果平均偏出训练集实际结果几个单位（有点像是方差），越低越好
- Adjusted R-squared：在多元线性模型里，使用Adjusted R-squared来评价模型，越高越好
- F-statistic：F检验给出了模型的整体重要性评估，它评估至少一个预测变量是否具有非零系数
得到预测结果后，我们就可以使用RMSE()和R2()对模型在测试集中的表现来进行评估：
# Model performance
# (a) Prediction error, RMSE
RMSE(predictions, test.data$sales)
# [1] 1.965508
# (b) R-square
R2(predictions, test.data$sales)
# [1] 0.9049049
RSEM：类似以上的Residual standard error，是针对测试集得到的结果，越低越好
R2:类似以上的Adjusted R-squared，是针对测试集得到的结果，越高越好
非线性回归
在某些情况下，数据结果和预测变量的关系并不是线性的，那么就有一些方法拓展了线性模型来应对这些情况：
- Polynomial regression：多项式回归，它将多项式（平方，立方）加入到模型里
- Spline regression：样条回归，使用多个多项式线段(polynomial segments)从而使得曲线更加平滑，限定这些segments的值成为Knots
- Generalized additive models (GAM)：广义加性模型，通过自动选择Knots来拟合样条模型
那么我们这里学习就是如何进行非线性回归（不介绍原理），并且比较模型的优劣从而选择一个最适合数据的模型（使用RMSE和R2指标）
读入数据
tidyverse 包用以数据的处理和可视化
caret 包含了一些简单回归函数
示例数据为来自MASS的BOSTON数据，我们将尝试使用1stata变量来预测mdev
library(tidyverse)
library(caret)
theme_set(theme_classic())
# 读入示例数据，来自MASS包的BOSTON数据
# Load the data
data("Boston", package = "MASS")
# Split the data into training and test set
set.seed(111)
training.samples <- Boston$medv %>%
createDataPartition(p = 0.8, list = FALSE)
train.data <- Boston[training.samples, ]
test.data <- Boston[-training.samples, ]
我们首先可以使用ggplot2来画一张scatter plot看一下两者之间的联系：
ggplot(train.data, aes(lstat, medv) ) +
geom_point() +
stat_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
# 我们可以看到这里ggplot2替我们自动选择了一种非线性回归模型~loess regression
线性模型的表现
那么我们首先使用线性模型来测试一下其表现，标准的线性模型可以写成：
medv = b_0 + b_1*lstat
# Build the model
model <- lm(medv ~ lstat, data = train.data)
summary(model)
# Call:
# lm(formula = medv ~ lstat, data = train.data)
#
# Residuals:
# Min 1Q Median 3Q Max
# -15.303 -4.086 -1.573 1.960 24.384
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 34.78563 0.63794 54.53 <2e-16
# lstat -0.96223 0.04401 -21.86 <2e-16
#
# Residual standard error: 6.389 on 405 degrees of freedom
# Multiple R-squared: 0.5413, Adjusted R-squared: 0.5402
# F-statistic: 478 on 1 and 405 DF, p-value: < 2.2e-16
# Make predictions
predictions <- model %>% predict(test.data)
# Model performance
data.frame(
RMSE = RMSE(predictions, test.data$medv),
R2 = R2(predictions, test.data$medv)
)
# RMSE R2
# 1 5.455969 0.5611306
可视化的画，只需要在stat_smooth()函数内指明拟合模型：
ggplot(train.data, aes(lstat, medv) ) +
geom_point() +
stat_smooth(method = lm, formula = y ~ x)
多项式回归模型的表现
多项式回归将多项式添加到了回归方程里：
mdev= b_0 + b_1*1stat+b_2*1stat^2+...+b_n*1stat^n
下面我们先测试一下二项式，可以看到加入二项式后显著性依旧很高，说明添加二项式可以提高模型的准确率：
model <- lm(medv ~ poly(lstat, 2, raw = TRUE), data = train.data)
summary(model)
# Call:
# lm(formula = medv ~ poly(lstat, 2, raw = TRUE), data = train.data)
#
# Residuals:
# Min 1Q Median 3Q Max
# -15.4151 -3.8283 -0.4949 2.3619 25.3261
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 43.211355 0.975879 44.28 <2e-16
# poly(lstat, 2, raw = TRUE)1 -2.361731 0.137782 -17.14 <2e-16
# poly(lstat, 2, raw = TRUE)2 0.043711 0.004127 10.59 <2e-16
#
# Residual standard error: 5.659 on 404 degrees of freedom
# Multiple R-squared: 0.641, Adjusted R-squared: 0.6392
# F-statistic: 360.7 on 2 and 404 DF, p-value: < 2.2e-16
那么我们最多可以加到几项式呢？
model <- lm(medv ~ poly(lstat, 10, raw = TRUE), data = train.data)
summary(model)
# Call:
# lm(formula = medv ~ poly(lstat, 10, raw = TRUE), data = train.data)
#
# Residuals:
# Min 1Q Median 3Q Max
# -14.2679 -3.0880 -0.5218 2.1185 26.2338
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 7.285e+00 2.899e+01 0.251 0.8017
# poly(lstat, 10, raw = TRUE)1 5.325e+01 2.992e+01 1.780 0.0759
# poly(lstat, 10, raw = TRUE)2 -2.601e+01 1.238e+01 -2.101 0.0363
# poly(lstat, 10, raw = TRUE)3 5.865e+00 2.728e+00 2.150 0.0322
# poly(lstat, 10, raw = TRUE)4 -7.584e-01 3.593e-01 -2.111 0.0354
# poly(lstat, 10, raw = TRUE)5 6.108e-02 2.989e-02 2.044 0.0416
# poly(lstat, 10, raw = TRUE)6 -3.168e-03 1.607e-03 -1.971 0.0494
# poly(lstat, 10, raw = TRUE)7 1.058e-04 5.562e-05 1.902 0.0579
# poly(lstat, 10, raw = TRUE)8 -2.197e-06 1.195e-06 -1.838 0.0668
# poly(lstat, 10, raw = TRUE)9 2.578e-08 1.449e-08 1.779 0.0760
# poly(lstat, 10, raw = TRUE)10 -1.305e-10 7.566e-11 -1.725 0.0854
#
# Residual standard error: 5.343 on 396 degrees of freedom
# Multiple R-squared: 0.6863, Adjusted R-squared: 0.6784
# F-statistic: 86.64 on 10 and 396 DF, p-value: < 2.2e-16
这里我们可以看出n=6（但将degree改为6以后发现事实上还是5）：
# Build the model
model <- lm(medv ~ poly(lstat, 5, raw = TRUE), data = train.data)
# Make predictions
predictions <- model %>% predict(test.data)
# Model performance
data.frame(
RMSE = RMSE(predictions, test.data$medv),
R2 = R2(predictions, test.data$medv)
)
# RMSE R2
# 1 4.551925 0.6912013
可以看出使用多项式回归拟合后表现好于简单的线性回归模型，可视化如下：
ggplot(train.data, aes(lstat, medv) ) +
geom_point() +
stat_smooth(method = lm, formula = y ~ poly(x, 5, raw = TRUE))
Log转换
当然面对一些非线性的数据，你也可以使用Log转化进行拟合（尝试）：
# Build the model
model <- lm(medv ~ log(lstat), data = train.data)
summary(model)
# Call:
# lm(formula = medv ~ log(lstat), data = train.data)
#
# Residuals:
# Min 1Q Median 3Q Max
# -14.5104 -3.4671 -0.6271 2.2332 25.9824
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 52.3915 1.0828 48.39 <2e-16
# log(lstat) -12.5858 0.4441 -28.34 <2e-16
#
# Residual standard error: 5.462 on 405 degrees of freedom
# Multiple R-squared: 0.6648, Adjusted R-squared: 0.664
# F-statistic: 803.3 on 1 and 405 DF, p-value: < 2.2e-16
# Make predictions
predictions <- model %>% predict(test.data)
# Model performance
data.frame(
RMSE = RMSE(predictions, test.data$medv),
R2 = R2(predictions, test.data$medv)
)
# RMSE R2
# 1 4.749176 0.6650908
ggplot(train.data, aes(lstat, medv) ) +
geom_point() +
stat_smooth(method = lm, formula = y ~ log(x))
样条回归
多项式回归仅捕获非线性关系中的一定量的曲率。，建模非线性关系的另一种方法（通常是更好的方法）是使用样条回归，简单的解释来说，就是把数据分为多端，每段各自使用多项式拟合，各个多项式回归交接的点便是 Knots；
R包splines里就提供了 bs函数(B-Spline Basis for Polynomial Splines)来创建样条回归，不过你需要指定两个参数：
- 多项式回归的degree，也就是多项式的个数
- knots的位置（一般来说不特定，根据数据的复杂程度来看分成多少段）
library(splines)
# Build the model
knots <- quantile(train.data$lstat, p = c(0.25, 0.5, 0.75))
model <- lm (medv ~ bs(lstat, knots = knots), data = train.data)
summary(model)
# Call:
# lm(formula = medv ~ bs(lstat, knots = knots), data = train.data)
#
# Residuals:
# Min 1Q Median 3Q Max
# -13.8751 -3.1449 -0.7017 2.1611 26.6803
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 50.802 2.719 18.687 < 2e-16
# bs(lstat, knots = knots)1 -14.298 4.191 -3.412 0.000711
# bs(lstat, knots = knots)2 -26.684 2.619 -10.189 < 2e-16
# bs(lstat, knots = knots)3 -28.081 3.104 -9.046 < 2e-16
# bs(lstat, knots = knots)4 -41.369 3.305 -12.519 < 2e-16
# bs(lstat, knots = knots)5 -38.778 4.545 -8.532 3.01e-16
# bs(lstat, knots = knots)6 -39.742 4.393 -9.046 < 2e-16
#
# Residual standard error: 5.354 on 400 degrees of freedom
# Multiple R-squared: 0.6818, Adjusted R-squared: 0.6771
# F-statistic: 142.9 on 6 and 400 DF, p-value: < 2.2e-16
# Make predictions
predictions <- model %>% predict(test.data)
# Model performance
data.frame(
RMSE = RMSE(predictions, test.data$medv),
R2 = R2(predictions, test.data$medv)
)
# RMSE R2
# 1 4.58503 0.6868804
可视化的话就需要调用splines包里的bs函数：
ggplot(train.data, aes(lstat, medv) ) +
geom_point() +
stat_smooth(method = lm, formula = y ~ splines::bs(x, df = 3))
广义加性模型
当你在探索非线性回归的时候，你可能会困扰于多项式回归不够全面，样条回归又需要指定knots，那么GAM(Generalized additive models)就是为了解决这些问题而被开发出来的，A review of spline function procedures in R(https://doi.org/10.1186/s12874-019-0666-3)这篇文章也详解了spline regression在R语言里面的运用，当然也包括了GAM，有兴趣的可以了解一下
mgcv包可以完成GAM的回归拟合，注意以上都是通过lm()函数对数据进行拟合，这里是使用mgcv包的gam函数进行回归拟合，mgcv包的s函数负责选定最优的knots：
library(mgcv)
# Build the model
model <- gam(medv ~ s(lstat), data = train.data)
summary(model)
# Family: gaussian
# Link function: identity
#
# Formula:
# medv ~ s(lstat)
#
# Parametric coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 22.6781 0.2646 85.72 <2e-16
#
# Approximate significance of smooth terms:
# edf Ref.df F p-value
# s(lstat) 7.513 8.448 102 <2e-16
#
# R-sq.(adj) = 0.679 Deviance explained = 68.5%
# GCV = 29.097 Scale est. = 28.488 n = 407
# Make predictions
predictions <- model %>% predict(test.data)
# Model performance
data.frame(
RMSE = RMSE(predictions, test.data$medv),
R2 = R2(predictions, test.data$medv)
)
# RMSE R2
# 1 4.633751 0.6807803
同样，在可视化上，stat_smooth()函数里包括了gam：
ggplot(train.data, aes(lstat, medv) ) +
geom_point() +
stat_smooth(method = gam, formula = y ~ s(x))
模型比较
| 模型 | R2(test.data) | RMSE(test.data) |
|---|---|---|
| 线性回归 | 0.5611306 | 5.455969 |
| 多项式回归 | 0.6912013 | 4.551925 |
| Log转换 | 0.6650908 | 4.749176 |
| 样条回归 | 0.6868804 | 4.58503 |
| 广义加性模型 | 0.6807803 | 4.633751 |
这么看来，对于这批数据，第一次测试出来的结果显多项式回归+样条回归+广义加性模型相对于线性模型和Log转换是有优越性的，并且看起来多项式回归给出的拟合是最好的（思考题：Test error vs. Training error）
推荐阅读
Splines, Additive Models, Model Selection and Validation: https://www.andrew.cmu.edu/user/achoulde/95791/lectures/lecture02/lecture02-95791.pdf
