R回归分析

qq_57518649

已于 2023-05-01 00:08:27 修改

阅读量4.3k

点赞数 11

文章标签： r语言回归开发语言前端

于 2023-04-25 22:21:30 首次发布

本文链接：https://blog.csdn.net/qq_57518649/article/details/130373844

版权

回归分析是处理变量间相互依存关系的一种有效方法。线性回归是最基本的回归方法。适合探索两个变量或多个变量之间的依存关系，如：一元线性回归和多元线性回归。

# 1.一元线性回归

1.1 模型：

y=a+bx+e，e是误差项。

最常用的估计参数a，b的方法是普通最小二乘法，即通过残差平方和Q最小来确定回归直线。

1.2 lm()函数实现回归：

# 一般形式
lm(formula,data,subset,weights,method="qr")
formula:指定用于拟合的模型形式的公式，一般以“响应变量~解释变量”的形式。
data:用于回归的数据对象
weights：指定用于回归的每个观测值权重，一般用于加权最小加权二乘回归。
method：指定用于回归拟合的方法，qr是拟合。

1.3 实例：iris数据集的一元回归建模

(lm1<-lm(Sepal.Length~Petal.Width,data=iris))
# 结果
Call:
lm(formula = Sepal.Length ~ Petal.Width, data = iris)

Coefficients:
(Intercept)  Petal.Width  
     4.7776       0.8886 
其中，Intercept为a的估计值，0.8886为b的估计值

回归拟合模型为：Sepal.Length=4.7776+0.8886*Petal.Width

1.4 对模型进行统计分析

summary(lm1)
-----结果-----
Call:
lm(formula = Sepal.Length ~ Petal.Width, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.38822 -0.29358 -0.04393  0.26429  1.34521 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       4.77763    0.07293   65.51   <2e-16 ***
iris$Petal.Width  0.88858    0.05137   17.30   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.478 on 148 degrees of freedom
Multiple R-squared:  0.669,	Adjusted R-squared:  0.6668 
F-statistic: 299.2 on 1 and 148 DF,  p-value: < 2.2e-16

1.4.1 参数

Residuals表示残差，是真实值与估计值的差值。DF：自由度。

Coefficients表示系数，表中a为4.777763，b为0.88858。

Std. Error为学生化残差，t/p value主要看p-value，若小于0.05，则表示在显著性水平为0.05的条件下通过显著性检验，即有95%的把握相信参数估计值通过显著性检验，上述a，b均通过检验。

1.4.2 模型优劣判断

利用R-squared判断：

Multiple R-squared（原始值）: 0.669, Adjusted R-squared（调整后）: 0.6668

通常选用Adjusted R-squared，表示模型能够解释原始变量信息的66.7%。在经济领域中，超过60%可以认为模型效果不错。

利用回归诊断图判断：

par(mfrow=c(2,2))
plot(lm1)

结果：

第一个图是残差值与拟合值的散点图，散点均匀分布在x=0的坐标轴两侧表明数据拟合效果较好。

第二个图是qq图，判断原始数据是否符合正态分布。若均匀分布在对角线上，则说明数据较好服从正态分布。

第三个图是标准化残差和拟合值的散点图，和第一个图用途一致。

第四个图是cook距离，用于判断模型中的异常值，图中123和132均为异常值。

~~判断出异常值后，剔除进行回归~~

(lm2<-lm(Sepal.Length~Petal.Width,data=iris[-c(123,132),]))
------结果-------
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.78897    0.06992   68.49   <2e-16 ***
Petal.Width  0.86488    0.04958   17.44   <2e-16 ***
Multiple R-squared:  0.6758,	Adjusted R-squared:  0.6735

相较于0.6668，删除异常值后的调整R-squared更大，对原始数据解释程度更大，拟合效果更好。

# 2.多元线性回归

2.1 模型

实际上，一个响应变量往往伴随多个解释变量。假设响应变量y与n个解释变量x1、...、xn线性相关，多元回归模型为：

$y=\beta_{_0}+\beta_{_1}x_{_1}+\beta_{_2}x_{_2}+\cdots \beta_{_n}x_{_n}+\varepsilon$

第t组数据对应模型：

$y_{_t}=\beta_{_0}+\beta_{_1}x_{_t1}+\beta_{_2}x_{_t2}+\cdots \beta_{_n}x_{_tn}+\varepsilon_{_t}$

其中：

$E(\varepsilon_{_t})=0,Var(\varepsilon)=\delta^{2},Cov(\varepsilon_{_i},\varepsilon_{_j})=0$

矩阵形式：

$Y=C\beta +\varepsilon$

最小二乘估计参数：

$\widehat{\beta}=(\beta_{_0},\beta_{_1},\cdots ,\beta_{_n}) =(C\acute{}C)^{-1}C\acute{}Y$

2.2 实例

2.2.1 lm()函数,根据iris做实例分析。

library(carData)
library(car)
scatterplotMatrix(iris[,1:4],col = "red")

图中明显观察到：Petal.Length与Sepal.Length、Petal.Width呈明显的正相关关系，与sepal.Width呈负相关关系，所有相关关系均为线性关系。

2.2.2 相关系数矩阵

corr=cor(iris[,1:4])
corr

2.2.3 构建模型

将Petal.Length作为响应变量，其他三个变量作为解释变量进行多元线性回归是合理选择。

lm2<-lm(Petal.Length~.,data=iris[,1:4])   
# "."等同于"Sepal.Length+Sepal.Width+Petal.Width"
summary(lm2)
--------结果--------
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.26271    0.29741  -0.883    0.379      
Sepal.Length  0.72914    0.05832  12.502   <2e-16 ***
Sepal.Width  -0.64601    0.06850  -9.431   <2e-16 ***
Petal.Width   1.44679    0.06761  21.399   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.319 on 146 degrees of freedom
Multiple R-squared:  0.968,	Adjusted R-squared:  0.9674 
F-statistic:  1473 on 3 and 146 DF,  p-value: < 2.2e-16

Intercept>0.05，未通过显著性检验，因此需要将回归模型进行修正。

lm3<-lm(Petal.Length~.+0,data=iris[,1:4])
summary(lm3)
--------结果------
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
Sepal.Length  0.69421    0.04284   16.20   <2e-16 ***
Sepal.Width  -0.67286    0.06134  -10.97   <2e-16 ***
Petal.Width   1.46802    0.06315   23.25   <2e-16 ***
Multiple R-squared:  0.9942,	Adjusted R-squared:  0.9941

去掉截距项后，各参数通过显著性检验，且回归的Adjusted R-squared为0.9941，说明回归效果不错。回归方程为：

P.L=0.69421*S.L-0.67286*S.W+1.46802*P.W

由第一个子图可看出拟合值的残差不是完全服从正态分布，而是较为明显的分成两部分散点，每部分都分布在x=0的坐标轴两侧。

2.2.4 模型诊断

为检验拟合值的同方差性，下面利用ncvTest()进行诊断。

par(mfrow=c(2,2))
------结果-------
Non-constant Variance Score Test   非恒定方差分数检验
Variance formula: ~ fitted.values 
Chisquare = 7.831356, Df = 1, p = 0.0051348

由于检验的P值小于0.05，说明回归拟合值不存在同方差性。

此外，还要考虑解释变量之间的交互作用。

lm4<-lm(Petal.Length~Sepal.Length*Petal.Width+Sepal.Width,data=iris[,1:4])
summary(lm4)
------结果-------
Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)              -1.14477    0.54103  -2.116   0.0361 *  
Sepal.Length              0.87530    0.09482   9.231 3.00e-16 ***
Petal.Width               2.03589    0.31038   6.559 8.91e-10 ***
Sepal.Width              -0.61154    0.07013  -8.719 5.93e-15 ***
Sepal.Length:Petal.Width -0.10424    0.05362  -1.944   0.0539 .  
Multiple R-squared:  0.9688,	Adjusted R-squared:  0.968

若各参数通过显著性检验且R-squared变大，则该交互项就可以引入模型中。从结果中观察，lm4模型的R-squared比lm3模型小，认为lm3是最适合的模型，可以进行后续数据探索：预测等。

2.2.5 模型预测

pred<-predict(lm3,data=iris[,1:4])
----对比预测值与实际值-------
par(mfrow=c(1,2))
plot(iris$Petal.Length)
plot(pred)

两图对比，lm3的预测效果非常好。绘制到一张图里观察差值。

2.3 变量的选择

在实际数据分析中，经常会遇到变量非常多的情况，即数据的维度很高，也称为“维数灾难”。如生物统计领域，对于回归建模，少而精的变量建模极为重要，如何选择变量子集是解决问题的关键。下面给出几个变量选择方法，如：逐步回归、岭回归和lasso回归法。

2.3.1 逐步回归

逐步选择，反复添加删除模型中的变量，达到优化模型的目的，需要确定阈值即算法停止的标准。该算法利用AIC准则（赤池信息准则）来度量添加变量的效果，AIC准则的表达式为：

AIC=-2*log(L)+k*edf

其中，L代表似然值，edf代表等效自由度。R中逐步回归函数为step()。

step(object，scope，scale=0，direction=c("both","backward","forward")，trace，steps，k)

object：指定模型的对象，如模型lm。scope：指定变量选择的上下界。

scale：回归模型和方差分析模型中定义的AIC所需要的值。

direction：变量回归方法，向后、向前、同时进行。

trace：是否输出变量选择过程，默认为１，输出。

steps：最大迭代次数。k：惩罚计算中自由度的倍数，默认值为２。

实例

　　利用swiss数据集逐步回归，swiss数据集是瑞士生育率和社会指标，47行观测值，６个变量。

data(swiss)
swiss
summary(swiss)　　＃　统计描述
# 绘制相关矩阵图
corr<-cor(swiss)
library(corrplot)
corrplot(corr)

　　根据相关图，选择Fertility作为响应变量，其余作为解释变量进行回归分析。Examination和Education之间的相关性较强，即解释变量之间存在多重共线性，后续回归分析会存在一些问题。

lm5<-lm(Fertility~.,data=swiss)　　＃　全模型
summary(lm5)
－－－－结果－－－－
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
Examination      -0.25801    0.25388  -1.016  0.31546    
Education        -0.87094    0.18303  -4.758 2.43e-05 ***
Catholic          0.10412    0.03526   2.953  0.00519 ** 
Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
Multiple R-squared:  0.7067,	Adjusted R-squared:  0.671

结果显示，Examination的系数没有通过显著性检验，而与其相关性较强的变量Eduction的p值非常小，为解决此问题，接下来用函数step()进行逐步回归。

tstep<-step(lm5)
-----------结果-----------
Start:  AIC=190.69
Fertility ~ Agriculture + Examination + Education + Catholic + 
    Infant.Mortality
                   Df Sum of Sq    RSS    AIC
- Examination       1     53.03 2158.1 189.86
<none>                          2105.0 190.69
- Agriculture       1    307.72 2412.8 195.10
- Infant.Mortality  1    408.75 2513.8 197.03
- Catholic          1    447.71 2552.8 197.75
- Education         1   1162.56 3267.6 209.36

Step:  AIC=189.86
Fertility ~ Agriculture + Education + Catholic + Infant.Mortality
                   Df Sum of Sq    RSS    AIC
<none>                          2158.1 189.86
- Agriculture       1    264.18 2422.2 193.29
- Infant.Mortality  1    409.81 2567.9 196.03
- Catholic          1    956.57 3114.6 205.10
- Education         1   2249.97 4408.0 221.43

逐步回归法选择的基准是AIC值最小，第一步回归的AIC是指剔除该变量后的AIC值，剔除变量Examination的AIC最小，所以第二步剔除了Examination变量，重新计算剔除各变量的AIC值，发现没有比AIC=189.86更小的，故迭代结束，最后只剔除Examination。

利用summary函数展示逐步回归的具体结果，发现全部参数通过显著性检验，且Adjusted R-squared为0.6707，说明模型是有效的。

summary(tstep)
--------------结果-------------
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      62.10131    9.60489   6.466 8.49e-08 ***
Agriculture      -0.15462    0.06819  -2.267  0.02857 *  
Education        -0.98026    0.14814  -6.617 5.14e-08 ***
Catholic          0.12467    0.02889   4.315 9.50e-05 ***
Infant.Mortality  1.07844    0.38187   2.824  0.00722 ** 
---
Multiple R-squared:  0.6993,	Adjusted R-squared:  0.6707

需要注意的是，逐步回归根据AIC最小选择变量，可能会出现挑选变量的参数通不过显著性检验的情况，需要用函数drop1()进行后续处理。

若剔除Examination后各变量系数通不过显著性检验，就人为剔除AIC第二小的Agriculture变量。剔除后各变量系数通过显著性检验且Adjusted R-squared为0.639，说明模型有效。

lm5<-lm(Fertility~.,data=swiss[,-c(2:3)])
summary(lm5)
tstep<-step(lm5)
summary(tstep)
---------结果-------------
Residuals:
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      48.67707    7.91908   6.147 2.24e-07 ***
Education        -0.75925    0.11680  -6.501 6.83e-08 ***
Catholic          0.09607    0.02722   3.530  0.00101 ** 
Infant.Mortality  1.29615    0.38699   3.349  0.00169 ** 
Multiple R-squared:  0.6625,	Adjusted R-squared:  0.639 

drop1(lm5)
---------结果-----
Model:
Fertility ~ Education + Catholic + Infant.Mortality
                 Df Sum of Sq    RSS    AIC
<none>                        2422.2 193.29   # 最小，说明已经是最优模型
Education         1   2380.38 4802.6 223.46
Catholic          1    701.74 3124.0 203.25
Infant.Mortality  1    631.92 3054.2 202.18

若drop1之后系数仍通不过显著性检验，则需考虑原始数据可能不适用于线性模型。

2.3.2 岭回归

岭回归是对每个参数添加一个惩罚项，基于最小化残差平方和与系数的惩罚项总和，系数的惩罚项总和是系数平方和的倍数。岭回归目的是寻找能使RSS最小的参数估计。

$RSS=\sum_{i=1}^{n}(y_{i}-\widehat{y_{i}})^{2}+\lambda \sum_{i=1}^{m}c_{i}^{2}$

library(MASS)   # 声明包：进行各种统计分析
lm.ridge(formula,data,lambda=0,model=FALSE)
formula：用于拟合的模型公式
data：用于岭回归的数据对象
lambda：默认为0

实例：基于R内置数据集lonely

某国16年的宏观经济数据，包括七个变量：物价指数平减后的国民生产总值(GNP.deflator)、国民生产总值(GNP)、失业人数(Unemployed)、军队人数(Armed.Forces)、人口总量(Population)、年份(Year)、就业人数(Employed)。将Employed设置为响应变量，其余变量作为解释变量进行多元回归。

lm6<-lm(Employed~.,data=longley)
summary(lm6)
----------结果---------
Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -3.482e+03  8.904e+02  -3.911 0.003560 ** 
GNP.deflator  1.506e-02  8.492e-02   0.177 0.863141    
GNP          -3.582e-02  3.349e-02  -1.070 0.312681    
Unemployed   -2.020e-02  4.884e-03  -4.136 0.002535 ** 
Armed.Forces -1.033e-02  2.143e-03  -4.822 0.000944 ***
Population   -5.110e-02  2.261e-01  -0.226 0.826212    
Year          1.829e+00  4.555e-01   4.016 0.003037 ** 
Multiple R-squared:  0.9955,	Adjusted R-squared:  0.9925

由上述输出结果可知，参数估计不全通过显著性检验，下面利用程辑包car中的函数vif()查看各自变量间的共线情况。

library(car)  # 声明包
vif(lm6)
------结果------
  GNP.deflator     GNP     Unemployed  Armed.Forces  Population 
   135.53244   1788.51348   33.61889      3.58893     399.15102 
     Year 
   758.98060

　　统计学中认为VIF值超过10表明存在多重共线性，上述变量除了Armed.Forces外其他变量均存在多重共线性，且GNP、Population和Year的VIF值超过200，说明存在严重多重共线性。

　　为解决多重共线性，利用函数lm.ridge（）做岭回归进行变量选择。

＃　绘制岭迹图
plot(lm.ridge(Employed~.,longley,lambda=seq(0,0.1,0.001)))
＃　返回三个不同估计值
select(lm.ridge(GNP.deflator~.,longley,lambda=seq(0,0.1,0.001)))
－－－－－－结果－－－－－－
modified HKB estimator is 0.006836982 
modified L-W estimator is 0.05267247 
smallest value of GCV  at 0.006 　　＃　最优值

变量选择准则：随着ｘ轴lambda的增大，回归系数即ｙ轴数据不稳定，振动趋于０的解释变量或回归系数稳定且绝对值很小的解释变量需要剔除。

根据变量选择原则，首先剔除稳定在0左右的黑色线所代表的GNP.deflator。

plot(lm.ridge(Employed~.-GNP.deflator,longley,lambda=seq(0,0.1,0.001)))
select(lm.ridge(GNP.deflator~.-GNP.deflator,longley,lambda=seq(0,0.1,0.001)))
----------结果-----------
modified HKB estimator is 0.006836982 
modified L-W estimator is 0.05267247 
smallest value of GCV  at 0.006    # lambda的最优值

再绘制岭迹图，发现各变量的岭迹仍不稳定，考虑剔除第二个岭迹图中黑色线代表的GNP变量。

plot(lm.ridge(Employed~.-GNP.deflator-GNP,longley,lambda=seq(0,0.1,0.001)))
select(lm.ridge(GNP.deflator~.-GNP.deflator-GNP,longley,lambda=seq(0,0.1,0.001)))
----------结果---------
modified HKB estimator is 0.01413139 
modified L-W estimator is 0.05929863 
smallest value of GCV  at 0.009

上图各变量的岭迹处于较为稳定的状态，此时模型中保留的变量为： Employed、Unemployed、 Armed.Forces和Population Year。对剩下变量进行回归。

lm7<-lm(Employed~.-GNP.deflator-GNP,longley)
summary(lm7)
----------结果---------
Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -2.446e+03  3.400e+02  -7.195 1.76e-05 ***
Unemployed   -1.500e-02  1.515e-03  -9.904 8.14e-07 ***
Armed.Forces -9.344e-03  1.855e-03  -5.038 0.000379 ***
Population   -2.287e-01  1.178e-01  -1.941 0.078295 .    
Year          1.302e+00  1.811e-01   7.191 1.77e-05 ***
***表示最显著,**表示非常显著,*表示比较显著,"."表示显著,而" "表示不显著。
Multiple R-squared:  0.9947,	Adjusted R-squared:  0.9927

各变量系数均通过显著性检验（显著性水平为0.05），且Adjusted R-squared为0.9927，说明该模型是有效的。

2.3.3 lasso回归

与岭回归类似，lasso也添加了针对回归参数的惩罚项，不同之处在于lasso选择惩罚的方式是：用系数的绝对值之和与残差平方和代替RSS，即：

$RSS=\sum_{i=1}^{n}(y_{i}-\widehat{y_{i}})^{2}+\lambda \sum_{i=1}^{m}\left | c_{i} \right |$

目的：找到使RSS最下时的参数估计，R中使用程辑包lars中的函数lars函数。

lars(x,y,type=c("lasso","lar","forward.stagewise","stepwise"))

参数：x为预测变量矩阵，y为响应变量向量，type是拟合模型的类型，lar表示最小交角回归,forward.stagewise极小向前逐段回归,stepwise逐步回归，默认值为lasso回归。

实例：lonely数据集的lasso回归建模

library(lars)
lar.1<-lars(as.matrix(longley[,1:6]),as.matrix(longley[,7]))
lar.1
-------结果--------
Call:
lars(x = as.matrix(longley[, 1:6]), y = as.matrix(longley[, 7]))
R-squared: 0.995 
Sequence of LASSO moves:
     GNP Unemployed Armed.Forces Year GNP Population GNP.deflator
Var    2          3            4    6  -2          5            1
Step   1          2            3    4   5          6            7
     GNP GNP.deflator GNP.deflator
Var    2           -1            1
Step   8            9           10

上述表明共迭代了10次以及每次选择的变量，而根据下面summary给出的cp值最小，选择第六次的迭代结果，即： GNP、Unemployed、Armed.Forces、Year和GNP Population五个变量。

summary(lar.1)
---------结果----------
LARS/LASSO
Call: lars(x = as.matrix(longley[, 1:6]), y = as.matrix(longley[, 7]))
   Df     Rss        Cp
0   1 185.009 1976.7120
1   2   6.642   59.4712
2   3   3.883   31.7832
3   4   3.468   29.3165
4   5   1.563   10.8183
5   4   1.339    6.4068
6   5   1.024    5.0186
7   6   0.998    6.7388
8   7   0.907    7.7615
9   6   0.847    5.1128
10  7   0.836    7.0000

绘制lasso回归图：plot（lar.1）

根据岭回归我们知道五个中有系数不通过显著性检验，因此，需要结合交叉验证来选择变量。

# cv.lars中K为10，做10折交叉验证
cva<-cv.lars(as.matrix(longley[,1:6]),as.matrix(longley[,7]),K=10,plot.it=TRUE)
best<-cva$index[which.min(cva$cv)]   # 最小cv值
coef<-coef.lars(lar.1,mode="fraction",s=best)   # 最小cv值对应系数
coef[coef!=0]   # 通过cv选择的变量
---------结果----------
cva[$index,$cv,$cv.error]
$mode
[1] "fraction"
> best
[1] 0.5454545
> coef[coef!=0]
  Unemployed Armed.Forces   Population         Year 
-0.014641359 -0.008436846 -0.135066279  1.155584381

由上述结果，交叉验证法选择的变量与岭回归的变量选择结果一致，故利用挑选的变量进行多元线性回归得到的回归模型是有效的。

2.3.4 Logistic回归

对于分类数据对象，普通的线性回归方法不再适用，引入改进的线性回归法即广义线性回归法用于构建logistic回归模型。logistic回归中响应变量是二分类变量，如男或女。Y：0/1。

函数glm()：glm(formula，family，data，weights)

参数介绍：

formula-模型；data：回归的数据对象；weights：用于指定每个观测值的权重。

family-指定干扰项的概率分布和模型的连接函数，默认family=gaussian，若需指定logistic回归，需设置binomial(link="logit")。

综合案例：iris数据集的逻辑回归建模

由于iris数据集中Species有三个种类，在建模之前先对数据集进行处理，将属于Setosa的数据剔除，然后用剩余数据建模分析。

iris<-iris[51:150,]   # 剔除Setosa分类数据
iris$Species<-ifelse(iris$Species=="virginica",1,0)   # 剩余响应变量0-1分类
log1<-glm(Species~.,family=binomial(link='logit'),data=iris)  # logistic
summary(log1) 
--------结果---------
Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept)   -42.638     25.707  -1.659   0.0972 .
Sepal.Length   -2.465      2.394  -1.030   0.3032  
Sepal.Width    -6.681      4.480  -1.491   0.1359  
Petal.Length    9.429      4.737   1.991   0.0465 *
Petal.Width    18.286      9.743   1.877   0.0605 .
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 138.629  on 99  degrees of freedom
Residual deviance:  11.899  on 95  degrees of freedom
AIC: 21.899
Number of Fisher Scoring iterations: 10

由结果，Sepal.Length和Petal.Length未通过显著性水平为0.05的检验。下面基于AIC准则进行逐步回归。

log2<-step(log1)
summary(log2)
－－－－结果－－－－
Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept)   -50.527     23.995  -2.106   0.0352 *
Sepal.Width    -8.376      4.761  -1.759   0.0785 .
Petal.Length    7.875      3.841   2.050   0.0403 *
Petal.Width    21.430     10.707   2.001   0.0453 *
AIC: 21.266

剔除掉Sepal.Length变量后，各解释变量系数均通过显著性水平为0.05的检验，说明构建模型有效。模型为：

Y=-50.527-8.376Sepal.Width+7.875Petal.Length+21.430Petal.Width

下面根据WMCR对模型的预测能力进行评估。原始数据中两种类比例为1:1，因此阈值设置为0.5。

pred<-predict(log2)
prob<-exp(pred)/(1+exp(pred))
yhat<-1*(prob>0.5)    # 0或1
table(iris$Species,yhat)　　＃　计数
－－－－结果－－－－－
   yhat
     0  1
  0 48  2　　　错判２
  1  1 49　　　错判１

利用图形展示模型的预测结果，一般采用ROC曲线来展示。AUC（Area Under Curve）被定义为ROC曲线下与坐标轴围成的面积，显然这个面积的数值不会大于1。由于ROC曲线一般处于y=x这条直线的上方，取值范围在0.5和1之间。AUC越接近1.0，检测方法真实性越高。等于0.5时，则真实性最低，无应用价值。本实例的AUC为0.9972，表明该预测模型价值较高。

library(ROCR)  # 画ROC曲线包
pred2<-prediction(pred,iris$Species)
performance(pred2,'auc')@y.values   # "Area under the ROC curve"
----结果------
[[1]]
[1] 0.9972

ROCR包中主要是两个class：prediction和performance。前者是将预测结果和真实标签组合在一起，生成一个prediction对象，然后在用performance函数，按照给定的评价方法，生成一个performance对象，最后直接对performance用plot函数就能绘制出相应的ROC曲线。

perf<-performance(pred2,'tpr','fpr')
plot(perf,col=2,type="l",lwd=2)   # (perf,col=2)：perf@y.values
f=function(x)
  {
  y=x;return(y)
  }
curve(f(x),0,1,col=4,lwd=2,lty=2,add=T)

结果如图，其中实曲线部分代表模型的预测效果，虚直线部分直线代表阈值变化，将实线与虚线进行比较，若实线在虚线上方说明回归效果好，本实例中实曲线位于虚直线上方且曲线坡度变化很快，说明模型效果很好。

qq_57518649

关注

11
点赞
踩
81

收藏

觉得还不错? 一键收藏
0
评论
R回归分析

回归分析是处理变量间相互依存关系的一种有效方法。线性回归是最基本的回归方法。适合探索两个变量或多个变量之间的依存关系，如：一元线性回归和多元线性回归。回归拟合模型为：Sepal.Length=4.7776+0.8886*Petal.Width。最常用的估计参数a，b的方法是普通最小二乘法，即通过残差平方和Q最小来确定回归直线。1.1 模型：y=a+bx+e，e是误差项。1.3 实例：iris数据集的一元回归建模。# 1.一元线性回归。
复制链接

扫一扫