Ｒ语言 一元线性回归

一元线性回归分析

首先介绍回归分析中最基础的情况：一元线性回归分析。它规定模型ｆ函数只能是ｙ=k*x+b的形式，即只使用一个变量ｘ（故称为一元）的线性形式来预测目标变量ｙ。

6.1.1引例

利用某网站历次促销活动中促销让利费用和销售金额的数据（单位是十万元），将使用该数据集来说明线性回归分析的应用。

使用如下语句来绘制其散点图：

cost<-c(1.8,1.2,0.4,0.5,2.5,2.5,1.5,1.2,1.6,1.0,1.5,0.7,1.0,0.8)

sales<-c(104,68,39,43,127,134,87,77,102,65,101,46,52,33)

data<-data.frame(cost=cost,sales=sales)

plot(data,pch=16,xlab="cost促销让利费用（十万元）",ylab="sales 促销销量（十万元)")

sol.lm<-lm(sales~cost,data)

abline(sol.lm,col="red")

由上图可以看出，促销让利费用cost变量和促销销量sales变量之间呈直线形式。

６.１.２　一元线性回归分析的原理及Ｒ语言实现

１.　最小二乘法

最小二乘法旨在使获得的模型最大限度地拟合样本数据

２．样本方差和协方差

３.　计算模型的参数ｋ和ｂ

（１）估算参数ｋ和ｂ的典型值

(k<-cov(cost,sales)/cov(cost,cost))#模型参数ｋ

(b<-mean(sales)-k*mean(cost))        #模型参数ｂ

也可使用回归方程建立函数lm来计算。代码如下：

(sol.lm<-lm(sales~cost,data))

Call:

lm(formula= sales ~ cost, data = data)

 

Coefficients:

(Intercept)         cost 

      13.82        48.60 

（２）估算参数ｋ和ｂ的取值范围

在计算得到的回归模型ｙ＝k*x+b中，系数ｂ和ｋ只是一个真实模型的典型估算值，其估算范围是：

[k_i-sd(k_i)t_α/2(n-2),k_i+sd(k_i)t_α/2(n-2)]

其中，k₀表示回归模型ｙ＝ｋ×ｘ＋ｂ中的ｂ，k₁表示ｋ。sd(k_i)是参数的标准差，可以使用summary(sol.lm)$coefficients[,2]直接读取，代码如下：

df<-sol.lm$df.residual

alpha=0.05

left<-summary(sol.lm)$coefficients[,1]-summary(sol.lm)$coefficients[,2]*qt(1-alpha/2,df);left

(Intercept)        cost

 1.667702  40.182861

right<-summary(sol.lm)$coefficients[,1]+summary(sol.lm)$coefficients[,2]*qt(1-alpha/2,df);right

(Intercept)        cost

 25.97978   57.01138

在上述代码中，df为自由度，等于样本数n减去自变量数p,再减１，即为n-2，该自由度可以通过sol.lm$df.residual直接读取，alpha是置信度，典型取值0.05,left是取值范围的最小值，right是取值范围的最大值。

４.　衡量相关程度

（１）相关系数ｒ和判定系数r²

对于引例，可以使用如下代码来计算相关参数。

(r<-cor(cost,sales))  #相关系数r

(r2<-r^2)             #判定系数r2

在ｌｍ函产生的结果中也包含了判定系数ｒ的信息。

summary(sol.lm)

Call:

lm(formula= sales ~ cost, data = data)

 

Residuals:

    Min     1Q  Median      3Q    Max

-19.701  -3.566  1.430   4.873  14.281

 

Coefficients:

            Estimate Std. Error t valuePr(>|t|)   

(Intercept)  13.824      5.579  2.478   0.0291 * 

cost          48.597      3.862 12.584 2.84e-08 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

 

Residualstandard error: 9.106 on 12 degrees of freedom

MultipleR-squared:  0.9296,       Adjusted R-squared:  0.9237

F-statistic:158.4 on 1 and 12 DF,  p-value: 2.843e-08

 

上述的Multiple R-squared就是判定系数ｒ２，该系数也可以使用summary(sol.lm)$r.squared来直接读取。

（２）修正判定系数adjusted.r²

事实上，判定系数ｒ２在用于多元回归分析时有一个缺点，自变量越多，判定系数ｒ２越大。为了消除这个缺点，可以引入修正判定系数adjusted.r²的概念。

lm函数产生的结果中也包含了修正判定系数adjusted.r²的信息

summary(sol.lm)

Call:

lm(formula= sales ~ cost, data = data)

 

Residuals:

    Min     1Q  Median      3Q    Max

-19.701  -3.566  1.430   4.873  14.281

 

Coefficients:

            Estimate Std. Error t valuePr(>|t|)   

(Intercept)  13.824      5.579  2.478   0.0291 * 

cost          48.597      3.862 12.584 2.84e-08 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

 

Residualstandard error: 9.106 on 12 degrees of freedom

MultipleR-squared:  0.9296,       AdjustedR-squared:  0.9237

F-statistic:158.4 on 1 and 12 DF,  p-value: 2.843e-08

上述的Adjusted R-squared就是修正判定系数，该系数也可以使用summary(sol.lm)$adj.r.squared来直接读取。

５.　回归系数的显著性检验

（１）Ｔ检验

Ｔ检验可以检验各个模型参数是否等于０，并计算其等于０时的概率。一般来说，当p.value小于０.０５，时，可以认定ｋ不会等于０，其模型结果可用并通过了检验。

summary(sol.lm)

Call:

lm(formula= sales ~ cost, data = data)

 

Residuals:

    Min     1Q  Median      3Q    Max

-19.701  -3.566  1.430   4.873  14.281

 

Coefficients:

            Estimate Std. Error t valuePr(>|t|)   

(Intercept)  13.824      5.579  2.478   0.0291 *  

cost          48.597      3.862 12.584 2.84e-08 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

 

Residualstandard error: 9.106 on 12 degrees of freedom

MultipleR-squared:  0.9296,       Adjusted R-squared:  0.9237

F-statistic:158.4 on 1 and 12 DF,  p-value: 2.843e-08

也可以使用如下代码直接读取结果:

summary(sol.lm)$coefficients[,4]

(Intercept)         cost

2.907896e-022.843305e-08

(2)　Ｆ检验

与Ｔ检验不同，Ｆ检验用于在整体上检验模型参数是否为０，并计算等于０的概率。

summary(sol.lm)

Call:

lm(formula= sales ~ cost, data = data)

 

Residuals:

    Min     1Q  Median      3Q    Max

-19.701  -3.566  1.430   4.873  14.281

 

Coefficients:

            Estimate Std. Error t valuePr(>|t|)   

(Intercept)  13.824      5.579  2.478   0.0291 * 

cost          48.597      3.862 12.584 2.84e-08 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

 

Residualstandard error: 9.106 on 12 degrees of freedom

MultipleR-squared:  0.9296,       Adjusted R-squared:  0.9237

F-statistic:158.4 on 1 and 12 DF,  p-value: 2.843e-08

其值为２．８４３ｅ－０８，远小于０．０５，表示通过了Ｆ检验。

在summary(sol.lm)中只给出了样本自由度＼自变量自由度以及Ｆ值，可以使用如下代码来直接读取Ｆ检验的ｐ-value值。

(f<-summary(sol.lm)$fstatistic[1])

(df1<-summary(sol.lm)$fstatistic[2])

(df2<-summary(sol.lm)$fstatistic[3])

pf(f,df1,df2,lower.tail=F)

       value

2.843305e-08

或者

1-pf(f,df1,df2)

       value

2.843305e-08

６.　模型误差（残差）

回归模型sol.lm的误差（残差residuals）,可用于体现样本点模型预测值与实际数据的差异程度，可通过sol.lm$residuals直接读取误差数据。对已一个正确的回归模型，其误差要服从正态分布。

shapiro.test(sol.lm$residuals)

Shapiro-Wilknormality test

 

data:  sol.lm$residuals

W= 0.9669, p-value = 0.8329

说明残差服从正态分布。

残差的标准误（Residual standard error）可以从整体上体现一个模型的误差情况，它可以用于不同模型间性能的对比。

summary(sol.lm)

Call:

lm(formula= sales ~ cost, data = data)

 

Residuals:

    Min     1Q  Median      3Q    Max

-19.701  -3.566  1.430   4.873  14.281

 

Coefficients:

            Estimate Std. Error t valuePr(>|t|)   

(Intercept)  13.824      5.579  2.478   0.0291 * 

cost          48.597      3.862 12.584 2.84e-08 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

 

Residual standarderror: 9.106 on 12 degrees of freedom

MultipleR-squared:  0.9296,       Adjusted R-squared:  0.9237

F-statistic:158.4 on 1 and 12 DF,  p-value: 2.843e-08

也可以使用summary(sol.lm)$sigma直接读取。

７.　预测

在Ｒ语言中，使用predict函数可以方便地计算预测值和取值范围。

下面的代码返回原有１４个样本的预测值数据。

predict(sol.lm)

     1        2         3         4         5        6         7         8         9        10        11        12        13

101.29856  72.14029 33.26259  38.12230 135.31655135.31655  86.71942  72.14029 91.57914  62.42086  86.71942 47.84173  62.42086

       14

 52.70144