第5章
Chapter 5线性回归——使用R
学习成果
通过本章的学习,您将能够: 
 解释回归分析,通常用于根据预测变量预测一个结果(目标或响应)变量的值；
 创建一个简单线性回归模型;
 使用残差与拟合图(residual vs. fitted plot)、正态分位数图(normal QQ plot)、位置尺度图(scale location plot)和残差与杠杆图(residuals vs. leverage plot)验证模型。
5.1概述
回归分析是一种估计变量之间关系的统计过程。当回归分析关注的是一个因变量(也称目标或响应变量)和一个或多个自变量(也称预测变量)之间的关系时,会包括许多建模和分析变量的技术。简单线性回归用于确定一个因变量和单个自变量之间的线性关系的程度。通常,回归分析用于以下三种目的: 对目标变量的预测(预告);对x和y之间关系的建模;对假设的测试。
5.2模型拟合
R语言中的模型是数据点序列的表示,具有看起来嘈杂的点云。模型拟合是指选择最适合的模型描述一组数据,R有不同类型的模型，下面列出这些模型及其对应的命令。
 线性模型(Linear Model,LM): lm()是R中的一个线性模型函数,用于构建一个简单回归模型。
 广义线性模型(Generalised Linear Model,GLM): 通过给出线性预测器的一种符号描述和误差分布的描述进行定义。
 混合效应线性模型(Linear Model for Mixed Effect,LME)。
 非线性最小二乘(Nonlinear Least Square,NLS): 确定一个非线性模型的参数的非线性(加权)最小二乘估计。
 广义加法模型(Generalised Additive Model,GAM): GAM只是一种统计模型,通常,响应变量和预测变量之间的线性关系会被几个用于数据建模和捕获非线性特征的非线性平滑函数代替。
每种模型都有特定的函数,数据点的分布基于描述模型的函数。
5.3线性回归
R中的线性回归由两个主要变量组成,它们通过两个变量的指数(幂)为1的方程相互关联。在数学上,当绘制图形时,线性关系表示一条直线。在存在非线性关系的情况下,变量的指数不等于1,它会在图上创建一条曲线。
一般的线性回归方程为
y=ax+b
其中,y是一个响应变量,x是一个预测变量,a和b是常数,称为系数。
大数据分析——基于R语言第5章线性回归——使用R5.3.1R中的lm()函数lm(formula, data, subset, weights, na.action,
method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,
singular.ok = TRUE, contrasts = NULL, offset, …)其中,
 formula表示x和y之间的关系。
 data包含模型中的变量。
 subset是一个可选的向量,它指定了用于模型拟合的观测值的一个子集。
 weights是一个可选的向量,它指定了模型拟合过程中的权重,它是一个数字向量或NULL。
 na.action是一个可选的函数,它指定了对于包含NA的数据如何反应的动作。
 method是用于拟合的方法。
 model,x,y,qr: 如果对应参数为TRUE,则返回模型矩阵、模型框和QR分解。
 singular.ok: 如果该参数为FALSE,则奇异拟合(singular fit)会出错。
 contrasts: 一个可选的列表偏移量,用于指定包含在线性预测器中的先前已知的成分。
线性回归中的lm()函数的简单语法为lm(formula, data),其中删除了可选参数。
确定一个学生数据集(student data set)的预测变量和响应变量之间的关系模型。预测向量存储着学生在学习上投入的时间,响应向量存储着学生的分数。
1. 检查数据集中的数据
考虑表5.1所给出的数据集D:\\student.csv,显示了学生在学习上投入的时间(NoOfHours)和分数(Freshmen_Score)。表5.1student数据集中的数据小时数成绩2552.5623653.570477小时数成绩4.5825755.5836856.5882. 将数据集中的数据读取到数据框中
使用read.table()函数以表格形式读取文件D:\\student.csv,并从中创建一个数据框HS,使用与文件中的行和字段变量相对应的大小写。> HS <- read.table("D:/student.csv", sep=",",header=TRUE)
> HS
NoOfHoursFreshmen_Score12.05522.56233.06543.57054.07764.58275.07585.58396.085106.5883. 查看数据框中数据的结果汇总
使用summary()函数生成结果汇总。这里,对所有数字变量计算出最小值、第一四分位数、中值、均值、第三四分位数和最大值。> summary(HS)
NoOfHoursFreshmen_Score
Min.: 2.000Min.: 55.001st Qu. : 3.1251st Qu. : 66.25Median : 4.250Median : 76.00Mean: 4.250Mean: 74.203rd Qu. : 5.3753rd Qu. : 82.75Max.: 6.500Max.: 88.004. 查看数据框的内部结构
显示R对象HS的内部结构,它表明有2个变量NoOfHours和Freshmen_Score的10组观测值。> str(HS)
'data.frame' : 10 obs. of 2 variables:$ NoOfHours : num 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5$ Freshmen_Score: int 55 62 65 70 77 82 75 83 85 885. 绘制R对象
绘制R对象HS,x轴使用HS$NoOfHours,y轴使用HS$Freshmen_Score,如图5.1所示。> plot (HS$NoOfHours, HS$Freshmen_Score)图5.1预测变量与响应变量的散点图
在图5.2所示的均值处(Freshmen_Score的均值为74.20)绘制一条水平线横穿该图。> abline (h=mean (HS$Freshmen_Score))图5.2在均值处绘制一条直线的预测变量与响应变量的散点图
当使用均值预测Freshmen_Score的分数时,在某些情况下可以观察到实际(观测)值和预测值之间的显著差异。
例如,第一个学生的Freshmen_Score为55,如果使用均值预测分数,则会预测它是74.20。在这里,观测值小于期望值。对于第10个学生,观测值为88,观测值大于预测分数74.20。
这表明要想预测期望的分数,则还应该考虑其他因素。
6. 相关系数r=n∑xy-∑x∑yn∑x2-∑x2|n∑y2-∑y2对于学生数据集,其计算为r=10×3295.5-42.5×742/square root of(［10×201.25 - 1806.25］×［10×56130 - 550564］)= 32955 - 31545/1488.05= 1410/1488.05= 0.947548805= 0.957. 使用R中的cor()函数确定线性关联的度和方向
线性关联的度(degree)和方向(direction)可以使用相关性确定。学习用时数和GPA分数之间关联的Pearson相关系数显示如下: > cor(HS$NoOfHours, HS$Freshmen_Score)
\[1\] 0.9542675这里的相关性值表明学习的用时数和学生的分数之间有一个强关联关系。
但是,对于相关性有以下一些需要注意的问题。
① 对于非线性关系,相关性不是一个度量关联的合适方法。如果要确定两个变量是否线性相关,则应该使用散点图。
② 异常值会影响Pearson相关性,箱线图可用于辨别异常值,异常值对Spearman相关性的影响是很小的。因此,如果没有合适的理由将异常值从分析中删除,则Spearman相关性是首选。
③ 相关性的取值接近于0表明变量不是线性关联的,但是这些变量可能依然是相关的，因此建议绘制这些数据。
④ 相关性并不意味着因果关系,即基于相关性的值并不能判断一个变量可以引起另一个变量。
⑤ 相关性分析仅仅有助于确定关联度。
因为相关性分析在确定因果关系时可能是不合适的,因此需要使用回归技术度量变量之间关系的本质。
使用一个数学模型y=f(X)表示回归模型,其中y是因变量,X是预测变量(x1,x2…xn)的集合。
一般而言,f(X)是线性或非线性的形式。
① 线性形式: f(X)=β0+β1x1+β2x2+…+βnxn+ε。
② 非线性形式: f(X)=β0+β1xp11+β2xp22+…+βnxpnn+ε。
一些常用的线性形式如下。
① 简单线性形式: 有一个预测变量和一个因变量，即f(X)=β0+β1x1+ε。
② 多重线性形式: 有多个预测变量和一个因变量，即f(X)=β0+β1x1+β2x2+…+βnxn+ε一些常用的非线性形式如下。
① 多项式形式: f(X)=β0+β1xp11+β2xp22+…+βnxpnn+ε。
② 二次型: f(X)=β0+β1x21+β2x1+ε。
③ Logistic形式: f(X)=11+e-(β0+β1x1+β2x2+…βnxn)+ε其中，β0,β1,β2…βn是回归系数,ε是预测中的误差。回归系数和预测中的误差都是实数。
当一个回归模型是线性形式时,这样的回归称为线性回归。类似地,当一个回归模型是非线性形式时,这样的回归称为非线性回归。
由于学生投入的学习时间与学生的分数之间的散点图表明二者存在线性关联,因此建立一个线性回归模型以量化这种关系的本质。
注意: 本章主要讨论线性回归。
在这个例子中,将使用学生投入的学习时间预测他/她的学习成绩。因此,Freshmen_Score可被看作因变量,而NoOfHours可被看作预测变量。这是一个简单线性回归的例子,因为其只有一个预测变量和一个因变量。
因此,用来预测学生成绩的回归模型可以被表示为Freshmen_Score=β0+(β1×NoOfHours)+ε8. 使用lm()函数创建线性模型
下面计算系数: 
(a) Intercept(b)HS$NoOfHours> x<-HS$NoOfHours
> y<-HS$Freshmen_Score
> n <- nrow (HS)
> xmean <- mean(HS$NoOfHours)
> ymean <- mean(HS$Freshmen_Score)
> xiyi <- x  y
> numerator <- sum(xiyi) - n  xmean  ymean
> denominator <- sum(x^2) - n  (xmean ^ 2)
> b1 <- numerator / denominator
> b0 <- ymean - b1  xmean
> b1
\[1\] 6.884848
> b0
\[1\] 44.93939使用lm()构建模型。在这里,HS$Freshmen_Score是响应变量或目标变量,HS$NoOfHours是预测变量。图5.3是对模型的可视化表示。> model_HS <- lm(HS$Freshmen_Score ~ HS$NoOfHours)
> model_HS
Call:
lm(formula = HS$Freshmen_Score ~ HS$NoOfHours)
Coefficients:
(Intercept)HS$NoOfHours44.9396.885
> summary(model_HS)
Call:
lm(formula = HS$Freshmen_Score ~ HS$NoOfHours)
Residuals:
Min1QMedian3QMax-4.3636-1.5803-0.37270.77126.0788
Coefficients
EstimateStd. Errort valuePr(>|t|)(Intercept)44.93943.421013.1361.07e-06 HS$NoOfHours6.88480.76269.0281.81e-05 
---
Signif. codes: 0 '" 0.001 '' 0.01 '' '0.05' '.' 0.1 '' 1
Residual standard error: 3.463 on 8 degrees of freedom
Multiple R-squared: 0.9106, Adjusted R-squared: 0.8995
F- statistic: 81.51 on 1 and 8 DF, p-value: 1.811e-05
> plot(HS$NoOfHours, HS$Freshmen_Score, co1="blue", main = 
"LinearRegression",
+ abline(lm(HS$Freshmen_Score ~ HS$NoOfHours)), cex = 1.3, pch = 16, xlab = "No of hours of study",
+ ylab = "Student Score")图5.3线性回归图
9. 对输出的说明
输出中的第一项是formula(lm(formula=HS$Freshmen_Score～HS$NoOfHours)),R使用它拟合数据。lm()是R中的一个线性模型函数,用于构建一个简单回归模型。HS$NoOfHours是预测变量,HS$Freshmen_Score是目标/响应变量。
模型输出中的下一项描述了residuals。实际观测值(例子中的HS$Freshmen_Score)和模型预测的响应值之间的差称为residuals。模型输出的residuals部分将其汇总成5点,即最小值、1Q(第一四分位数)、中值、3Q(第三四分位数)和最大值。在评估模型对数据的拟合程度时,应该在均值为0上寻找这些点之间的对称分布值，如表5.2所示。表5.2模型输出小时数成绩预测值残差值
（实际值-估计值）25558.70909-3.709092.56262.15152-0.1515236565.59394-0.593943.57069.036360.9636447772.478794.521214.58275.921216.07879 (maximum value)57579.36364-4.36364 (minimum value)5.58382.806060.1939468586.24848-1.248486.58889.69091-1.69091要计算5个汇总点,需要将数据以升序排列。
(-4.36364, -3.70909, -1.69091, -1.24848, -0.59394, -0.15152, 0.19394,
0.96364, 4.52121,6.07879)
最小值: -0.436364。
1Q: 在位置3.25处。
获取3.5位置处的值=(-1.69090-1.24848)/2=-1.46969。
获取3.25位置处的值=(-1.69090-1.46969)/2=-1.580295。
中值: (-0.59394 -0.15152)/2=-0.37273(中值在位置5.5处)。
3Q: 在位置7.75处。
获取位置7.5处的值=(0.19394+0.96364)/2=0.57879。
获取位置7.75处的值=(0.57879+0.96364)/2=0.771215。
最大值: 6.07879。
下一部分描述了模型的系数(coefficient)。从理论上来讲,在简单线性回归中,系数是两个未知的常数,表示线性模型中的截距(intercept)和斜率(slope)。
10. 系数: 估计值
估计值(estimate)包含两行。第一行是截距,所有预测变量X=0时是响应变量Y的均值。请注意,只有当模型中的每个X的实际值为0时均值才有用。系数的第二行是斜率或者HS_NoOfHours对Freshmen_Score的影响。模型中的斜率证明了NoOfHours每增加一小时,Freshmen_Score就会提高6.8848分。
11. 系数: 标准误差
标准误差(standard error)度量了从响应变量的实际平均值估计的系数变化的平均量。理想的情况是: 相对于其系数,这应该是一个较小的数。
12. 系数: tvalue
系数tvalue是指系数估计值距离0有多少标准差(standard deviation)的度量，它应该是远离0的,因为可以声明HS_NoOfHours和Freshmen_Score之间的关系是存在的。tvalue是系数除以标准误差((44.9394/3.4210)=13.1363)。一般来说,tvalue也用来计算pvalue。
13. 系数: Pr(>t)
在模型输出中发现的字母缩写pr(>t)是观测到的任何值等于或大于t的概率。偶然情况下的一个小的pvalue表明不太可能观察到预测变量(HS_NoOfHours)和响应变量(Freshmen_Score)之间的关系。典型的情况是: 5%及以下的pvalue是一个很好的临界点。
注意: signif.Codes与每个估计值相关联。
3颗星(或星号)表示一个非常显著的pvalue。
由标明的系数是pvalue<0.001的系数。
由标明的系数是pvalue<0.01的系数。
14. 残差标准误差(Residual Standard Error)
残差标准误差是对线性回归拟合质量的度量。从理论上而言,每个线性模型都假设包含一个误差项E,它可以防止从预测变量中完美地预测出响应变量。下面计算均方根误差(Root Mean Squared Error,RMSE),它是均方残差的平方根。
考虑如表5.3所示的学生数据集。表5.3学生数据集小时数成绩预测值残差值
(实际值-估计值)残差平方和25558.70909-3.7090913.757348632.56262.15152-0.151520.0229583136565.59394-0.593940.3527647243.57069.036360.963640.9286020547772.478794.5212120.441339864.58275.921216.0787936.9516878657579.36364-4.3636419.041354055.58382.806060.193940.03761272468586.24848-1.248481.558702316.58889.69091-1.690912.859176628注意: 以下部分将会展示利用predict()和resid()函数计算预测值和残差值的方法。
残差标准误差=(残差的平方之和/模型自由度)的平方根。
 残差的平方之和=95.95154715。
 模型的自由度=8。
自由度是由数据集的行数-列数或变量数指定的,student数据集中有10行和2列HS$NoOfHours和HS$Freshmen_Score,即10-2=8。
可以使用df.residual()函数在R中计算自由度。> df.residual (model_HS)
\[1\] 8残差标准误差=(95.95154715/8)的平方根。
残差标准误差=(11.99394339)的平方根。
残差标准误差=3.463227309。
15. 拟合优度(Multiple Rsquared)和修正的拟合优度(Adjusted Rsquared)
拟合优度 Rsquared(R2R2)统计提供了模型和实际数据拟合程度的一种度量,它的表示形式为方差比例(proportion of variance)。R2R2是一种用于预测变量(HS_NoOfHours)和响应/目标变量(Freshmen_Score)之间的线性关系的度量,它总是介于0和1之间(也就是说,一个接近于0的数字代表一种不能很好地解释响应变量中方差的回归,而一个接近1的数字则解释了响应变量中观测到的方差)。
拟合优度又称确定系数，它给出了有多少数据点落在回归方程形成的直线内的观点。该系数越大,当绘制数据点和直线时,直线会穿过更高比例的数据点。如果系数是0.80,那么80%的点应该落在回归线内。值1和0表明回归线分别表示全部数据或没有数据。一个较大的系数表明对观测值更好的拟合。