拟合线性模型

R 中最简单的模型就是线性模型,即使用一个线性模型描述两个随机变量在一系列假
设下的关系。在以下例子中,我们创建了一个线性函数,将 x 映射到 3+2*x(即 f(x))。
然后,生成一个正态分布的随机数值向量 x,再用 f(x)加上一些独立的噪声生成 y:
f <- function(x) 3 +2 * x
x <- rnorm(100)
y <- f(x) +0.5 * rnorm(100)

假设不知道 y 是如何通过 x 生成的,那么我们还能用线性模型还原它们之间的关系,
即还原线性模型的系数吗?以下代码中,lm( )用线性模型拟合 x 和 y 的关系。需要注意
的是 y~x 表达式告诉 lm( )这个线性回归是关于因变量 y 和单个回归元(regressor)x 的:
model1 <- lm(y ~x)
model1
##
## Call:
## lm(formula = y ~ x)
##
## Coefficients:
## (Intercept) x
## 2.969 1.972
真实模型的系数是 3(截距)和 2(斜率),用样本数据 x 和 y 拟合的模型系数为
2.9692146(截距)和 1.9716588(斜率),这和模型的真实系数非常接近。
我们将模型保存为 model1。要得到模型的系数,可以用以下代码:
coef(model1)
## (Intercept) x
## 2.969215 1.971659
也可以用 model1$查看模型的系数,因为 model1 本质上是一个列表。
此外,还可以调用 summary( )了解更多关于这个线型模型的统计性质:
summary(model1)
##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.96258 -0.31646 -0.04893 0.34962 1.08491
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.96921 0.04782 62.1 <2e-16 ***
## x 1.97166 0.05216 37.8 <2e-16 ***
## ---
## Signif. codes:
## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##
## Residual standard error: 0.476 on 98 degrees of freedom
## Multiple R-squared: 0.9358, Adjusted R-squared: 0.9352
## F-statistic: 1429 on 1 and 98 DF, p-value: < 2.2e-16
要理解以上统计内容,最好回顾一两本统计学教材的线性回归章节。图 7-31 将数据和
拟合模型放到了一起:
plot(x, y, main = "A simple linear regression")
abline(coef(model1), col = "blue")

 


图 7-31
在上述代码中,我们直接为 abline( )提供包含回归系数估计的两元素数值向量,
用来添加拟合的回归线。
接下来,我们可以调用 predict( )对所拟合模型进行预测。在一个标准误差下,预
测 x=-1 和 x=0.5 时的 y,运行以下代码:
predict(model1, list(x = c(-1, 0.5)), se.fit = TRUE)
## $fit
## 1 2
## 0.9975559 3.9550440
##
## $se.fit
## 1 2
## 0.06730363 0.05661319
##
## $df
## [1] 98
##
## $residual.scale
## [1] 0.4759621

预测结果是一个列表,其中包含了 y 的预测值($fit)、拟合值的标准误差($se.fit)、
自由度($df)和残差标准差($residual.scale)。
现在,我们已经基本上学会了如何对给定数据拟合线性模型,接下来就看一些现实生
活中的数据。下面的例子中,我们试着用不同复杂度的线性模型预测航班的飞行时间。其
中,对于预测飞行时间最为显著的变量是飞行距离。
首先,我们加载数据集并绘制飞行距离和飞行时间 air_time 的散点图。因为数据集
的记录条数非常多,我们用 pch="." 使每个数据点变得非常小:
data("flights", package = "nycflights13")
plot(air_time ~distance, data = flights,
pch = ".",
main = "flight speed plot")
图 7-32 清楚地表明飞行距离和飞行时间之间有正向相关关系,因此用线性模型拟合这
两个变量是合理的。

 


图 7-32
在将整个数据集输入线性模型之前,我们将数据集分成两部分:训练集和测试集。划分
测试集的目的是我们不仅想要进行模型的样本内评估,还要进行样本外评估。更具体地说,
我们把 75%的数据放入训练集,把剩下 25%的数据放入测试集。在以下代码中,sample( )
取出原始数据记录 75%的随机样本,并用 setdiff( )得到剩下的数据记录:
rows <- nrow(flights)
rows_id <- 1:rows
sample_id <- sample(rows_id, rows * 0.75, replace = FALSE)
flights_train <- flights[sample_id,]
flights_test <- flights[setdiff (rows_id, sample_id), ]
注意,setdiff(rows_id, sample_id)返回的是 rows_id 的索引而不是 sample_id

的索引。
现在,flights_train 是训练集,flights_test 是测试集。使用划分好的数据集,
模型的拟合和评估流程就很简单了。首先,我们使用训练集拟合模型,然后做一个样本内
预测来看样本内误差的大小:
model2 <- lm(air_time ~distance, data = flights_train)
predict2_train <- predict(model2, flights_train)
error2_train <- flights_train$air_time -predict2_train
为了计算误差的大小,我们定义一个 evaluate_error( )函数计算平均绝对误差和
误差的标准差:
evaluate_error <- function(x) {
c(abs_err = mean(abs(x), na.rm = TRUE),
std_dev = sd(x, na.rm = TRUE))
}
用这个函数,我们可以计算模型 model2 的样本内预测误差:
evaluate_ _error(error2_train)
## abs_err std_dev
## 9.413836 12.763126
平均绝对误差表示平均来看预测值偏离真实值的绝对值约为 9.41 分钟,且误差的标准
差为 12.8。
接下来,我们用模型对测试集进行预测,做一个简单的样本外评估:
predict2_test <- predict (model2, flights_test)
error2_test <- flights_test$air_time -predict2_test
evaluate_ _error(error2_test)
## abs_err std_dev
## 9.482135 12.838225
预测结果被保存在一个含预测值的数值向量中。平均绝对误差和误差的标准差都有轻
微的上升,表明样本外预测的质量没有明显地变差,看来 model2 没有过拟合。
因为 model2 只有一个回归元——路程(distance),很自然地就会想到包含更多回归
元是否会提高模型的预测效果。以下代码拟合一个新的线性模型,它不仅包含路程,还包
含航空公司(carrier)、月份(month)和起飞时间(dep_time):
model3 <- lm(air_time ~ carrier + distance + month + dep_time,
data = flights_train)

predict3_train <- predict(model3, flights_train)
error3_train <- flights_train$air_time - predict3_train
evaluate_ _error(error3_train)
## abs_err std_dev
## 9.312961 12.626790
样本内误差的大小和方差都小幅下降:
predict3_test <- predict(model3, flights_test)
error3_test <- flights_test$air_time - predict3_test
evaluate_ _error(error3_test)
## abs_err std_dev
## 9.38309 12.70168
另外,样本外误差也比 model2 看起来稍好一些。为了比较添加新回归元前后的样本
外误差分布,我们将两条密度曲线相互叠加到一起,如图 7-33 所示。
plot(density(error2_test, na.rm = TRUE),
main = "Empirical distributions of out-of-sample errors")
lines(density(error3_test, na.rm = TRUE), lty = 2)
legend("topright", legend = c("model2", "model3"),
lty = c(1, 2), cex = 0.8,
x.intersp = 0.6, y.intersp = 0.6)

 


图 7-33
由以上密度图可知,model2 到 model3 模型拟合效果的提升小到可以忽略,即没有
得到明显的改善。

posted @ 2019-01-22 14:38  NAVYSUMMER  阅读(1098)  评论(0编辑  收藏  举报
交流群 编程书籍