R语言数据分析与挖掘(第四章):回归分析(4)——logistic回归

DoubleHelix

发布于 2019-12-13 12:14:39

11.9K0

发布于 2019-12-13 12:14:39

前面我们介绍的回归方法，一般适用于数值型数据对象，对于分类数据类型就不再适用。对于分类数据对象，我们需要引入广义线性回归方法，比如logistic回归和poisson回归模型。这里我们介绍logistic回归。

logistic回归又称logistic回归分析，是一种广义的线性回归分析模型，常用于数据挖掘，疾病自动诊断，经济预测等领域。例如，探讨引发疾病的危险因素，并根据危险因素预测疾病发生的概率等。以胃癌病情分析为例，选择两组人群，一组是胃癌组，一组是非胃癌组，两组人群必定具有不同的体征与生活方式等。因此因变量就为是否胃癌，值为“是”或“否”，自变量就可以包括很多了，如年龄、性别、饮食习惯、幽门螺杆菌感染等。自变量既可以是连续的，也可以是分类的。然后通过logistic回归分析，可以得到自变量的权重，从而可以大致了解到底哪些因素是胃癌的危险因素。同时根据该权值可以根据危险因素预测一个人患癌症的可能性。

R语言中用于实现logistic回归的函数是glm()，其基本书写格式为:

glm(formula, family = gaussian, data, weights, subset,
    na.action, start = NULL, etastart, mustart, offset,
    control = list(...), model = TRUE, method = "glm.fit",
    x = FALSE, y = TRUE, singular.ok = TRUE, contrasts = NULL, ...)

参数介绍:

Formula:指定用于拟合的模型公式，类似于Im中的用法:

Family: 指定描述干扰项的概率分布和模型的连接函数，默认值为gaussian, 若需进行logistic同归，则需设置为binomial(link = "logit");

Data:指定用于回归的数据对象，可以是数据框、列表或能被强制转换为数据框的数据对象:

Weights:一个向量，用于指定每个观测值的权重:

Subset:一个向量，指定数据中需要包含在模型中的观测值;

Na.ction:一个函数，指定当数据中存在缺失值时的处理办法，用法与Im中的一致;

Start:一个数值型向量，用于指定现行预测器中参数的初始值;

Etastart:一个数值型向量，用于指定现行预测器的初始值;

Mustart:一个数值型向量，用于指定均值向量的初始值:

Offset:指定用于添加到线性项中的一组系数恒为1的项:

Contol:指定控制拟合过程的参数列表，其中epsilon 表示收敛的容忍度，maxit表示迭代的最大次数，trace 表示每次迭代是否打印具体信息；

Model: 逻辑值，指定是否返回“模型框架”，默认值为TRUE:

Method;指定用于拟合的方法，“glm.ft”表示用于拟合，“model.frame"表示可以返回模型框架；

X:逻辑值，指定是否返回“横型矩阵”，默认值为FALSE:

Y:逻辑值，制度是否能够返回响应变量，默认值为TRUE;

Contrasts:模型中因子对照的列表。

　　下面利用iris 数据集进行操作演练，由于iris数据集中的分类变量Specics中有三种元素:setosa、versicolor 和virginica，即鸢尾花的有三个不同的种类，在建模之前，先对数据集进行处理，将数据集中Species属于setosa类的数据剔除，然后利用剩余的数据进行建模分析，具体操作如下:

> iris<-iris[51:150,]
> iris$Species<-ifelse(iris$Species=="virginica",1,0)
> log1<-glm(Species~.,family=binomial(link='logit'),data=iris)
> summary(log1)

Call:
glm(formula = Species ~ ., family = binomial(link = "logit"), 
    data = iris)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.01105  -0.00541  -0.00001   0.00677   1.78065

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 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(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

上述代码表示:选择iris数据集中第51行到150行的数据，将该数据集中变量

Species列中记录为virginica 的替换为1,否则替换为0，然后利用清洗好的数据进行logistic回归；模型的输出结果显示:解释变量Sepal.Length和Sepal.Width没能通过显著性水平为0.05的检验。

下面基于前面介绍的AIC准则（R语言数据分析与挖掘(第四章):回归分析(3)——变量的选择）进行逐步回归：

> log2<-step(log1)
Start:  AIC=21.9
Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width

               Df Deviance    AIC
- Sepal.Length  1   13.266 21.266
<none>              11.899 21.899
- Sepal.Width   1   15.492 23.492
- Petal.Width   1   23.772 31.772
- Petal.Length  1   25.902 33.902

Step:  AIC=21.27
Species ~ Sepal.Width + Petal.Length + Petal.Width

               Df Deviance    AIC
<none>              13.266 21.266
- Sepal.Width   1   20.564 26.564
- Petal.Length  1   27.399 33.399
- Petal.Width   1   31.512 37.512
Warning messages:
1: glm.fit:拟合機率算出来是数值零或一 
2: glm.fit:拟合機率算出来是数值零或一 
> summary(log2)

Call:
glm(formula = Species ~ Sepal.Width + Petal.Length + Petal.Width, 
    family = binomial(link = "logit"), data = iris)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.75795  -0.00412   0.00000   0.00290   1.92193

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 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 138.629  on 99  degrees of freedom
Residual deviance:  13.266  on 96  degrees of freedom
AIC: 21.266

Number of Fisher Scoring iterations: 10

不难发现逐步回归剔除了变量：Sepal.Length,对逐步回归的结果进行详细展示，可以看到剩余变量的参数估计均通过了显著性水平的0.05的检验，说明构建模型得到了数据的支持。

下面根据WMCR对模型的预测能力进行评估：

> pred<-predict(log2)
> prob<-exp(pred)/(1+exp(pred))
> yhat<-1*(prob>0.5)
> table(iris$Species,yhat)
   yhat
     0  1
  0 48  2
  1  1 49

上述代码表示:首先将模型的预测结果存储到变量pred中，再根据前面介绍的模型进行logit变换的逆变换，输出结果存储到变量prob,此时该变量中的值即为响应变量取值为1的概率值，即变量Species=virginica的概率值，然后分别计算变量prob中大于0.5和小于等于0.5的记录总数，其中0.5即为阈值，由于原始的鸢尾花数据中,种类为versicolor和virginica的记录各有50条，故阈值a取值为0.5。最后利用函数table( )统计原始数据中的记录和预测结果的记录情况(“0”表示versicolor,“1”表示virginica), 不难发现，输出的表格中，数字“48”和“49”均表示预测正确的总数，数字“2”表示真实种类为versicolor, 而预测结果为virginica 的记录总数，

类似地，数字“1”表示真实种类为virginica,而预测结果为versicolor 的记录总数。除此之外，还可以利用图形展示模型的预测效果，业界一般采用ROC曲线对logistic 回归模型的效果进行刻画，R语言的RORC包中有专门的函数用于刻画ROC曲线，具体操作如下:

> library(ROCR)
> pred2<-prediction(pred,iris$Species)
> performance(pred2,'auc')@y.values
[[1]]
[1] 0.9972
>  perf<-performance(pred2,'tpr','fpr')
> plot(perf,col=2,type="l",lwd=2)
> f=function(x){ y=x;return(y)}
> curve(f(x),0,1,col=4,lwd=2,lty=2,add=T)