案例一：本文用例来自于John Maindonald所著的《Data Analysis and Graphics Using R》一书，其中所用的数据集是anesthetic，数据集来自于一组医学数据，其中变量conc表示麻醉剂的用量，move则表示手术病人是否有所移动，而我们用nomove做为因变量，因为研究的重点在于conc的增加是否会使nomove的概率增加。

首先载入数据集并读取部分文件，为了观察两个变量之间关系，我们可以利cdplot函数来绘制条件密度图

install.packages('DAAG')

library(lattice)

library(DAAG)

head(anesthetic)

move conc   logconc nomove

1    0  1.0 0.0000000      1

2    1  1.2 0.1823216      0

3    0  1.4 0.3364722      1

4    1  1.4 0.3364722      0

5    1  1.2 0.1823216      0

6    0  2.5 0.9162907      1

cdplot(factor(nomove)~conc,data=anesthetic,main='条件密度图',ylab='病人移动',xlab='麻醉剂量')

从图中可见，随着麻醉剂量加大，手术病人倾向于静止。下面利用logistic回归进行建模，得到intercept和conc的系数为-6.47和5.57，由此可见麻醉剂量超过1.16(6.47/5.57)时，病人静止概率超过50%。

anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic)

summary(anes1)

结果显示：

Call:

glm(formula = nomove ~ conc, family = binomial(link = 'logit'),

data = anesthetic)

Deviance Residuals:

Min        1Q    Median        3Q       Max

-1.76666  -0.74407   0.03413   0.68666   2.06900

Coefficients:

Estimate Std. Error    z value        Pr(>|z|)

(Intercept)   -6.469      2.418       -2.675         0.00748 **

conc             5.567      2.044        2.724          0.00645 **

---

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

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 41.455  on 29  degrees of freedom

Residual deviance: 27.754  on 28  degrees of freedom

AIC: 31.754

Number of Fisher Scoring iterations: 5

下面做出模型的ROC曲线

anes1=glm(nomove~conc,family=binomial(link='logit'),data=anesthetic)

对模型做出预测结果

pre=predict(anes1,type='response')

将预测概率pre和实际结果放在一个数据框中

data=data.frame(prob=pre,obs=anesthetic$nomove)

将预测概率按照从低到高排序

data=data[order(data$prob),]

n=nrow(data)

tpr=fpr=rep(0,n)

根据不同的临界值threshold来计算TPR和FPR，之后绘制成图

for (i in 1:n){

threshold=data$prob[i]

tp=sum(data$prob>threshold&data$obs==1)

fp=sum(data$prob>threshold&data$obs==0)

tn=sum(data$prob

fn=sum(data$prob

tpr[i]=tp/(tp+fn)  #真正率

fpr[i]=fp/(tn+fp)  #假正率

}

plot(fpr,tpr,type='l')

abline(a=0,b=1)

R中也有专门绘制ROC曲线的包，如常见的ROCR包，它不仅可以用来画图，还能计算ROC曲线下面面积AUC,以评价分类器的综合性能，该数值取0-1之间，越大越好。

library(ROCR)

pred=prediction(pre,anesthetic$nomove)

performance(pred,'auc')@y.values

perf=performance(pred,'tpr','fpr')

plot(perf)

还可以使用更加强大的pROC包，它可以方便的比较两个分类器，并且能自动标出最优临界点，图形看起来比较漂亮：

install.packages('pROC')

library(pROC)

modelroc=roc(anesthetic$nomove,pre)

plot(modelroc,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c('green','red'),max.auc.polygon=TRUE,auc.polygon.col='blue',print.thres=TRUE)

上面的方法是使用原始的0-1数据进行建模,即每一行数据均表示一个个体，另一种是使用汇总数据进行建模，先将原始数据按下面步骤进行汇总

anestot=aggregate(anesthetic[,c('move','nomove')],by=list(conc=anesthetic$conc),FUN=sum)

结果如下：

conc move nomove

1       0.8     6     1

2       1.0     4     1

3       1.2     2     4

4      1.4      2     4

5      1.6      0     4

6      2.5     0      2

anestot$conc=as.numeric(as.character(anestot$conc))

anestot$total=apply(anestot[,c('move','nomove')],1,sum)

anestot$total

[1] 7 5 6 6 4 2

anestot$prop=anestot$nomove/anestot$total

anestot$prop

[1] 0.1428571 0.2000000 0.6666667 0.6666667 1.0000000 1.0000000

对于汇总数据，有两种方法可以得到同样的结果，一种是将两种结果的向量合并做为因变量，如anes2模型。另一种是将比率做为因变量，总量做为权重进行建模，如anes3模型。这两种建模结果是一样的。

anes2=glm(cbind(nomove,move)~conc,family=binomial(link='logit'),data=anestot)

summary(anes2)

结果显示如下：

Call:

glm(formula = cbind(nomove, move) ~ conc, family = binomial(link = 'logit'),

data = anestot)

Deviance Residuals:

1 2 3 4 5 6

0.20147 -0.45367 0.56890 -0.70000 0.81838 0.04826

Coefficients:

Estimate  Std. Error  z value    Pr(>z)

(Intercept) -6.469       2.419    -2.675     0.00748 **

conc          5.567        2.044      2.724      0.00645 **

---

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

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 15.4334 on 5 degrees of freedom

Residual deviance: 1.7321 on 4 degrees of freedom

AIC: 13.811

Number of Fisher Scoring iterations: 5

anes3=glm(prop~conc,family=binomial(link='logit'),weights=total,data=anestot)

结果和上面的一样。

根据logistic模型，我们可以使用predict函数来预测结果，下面根据上述模型来绘图

x=seq(from=0,to=3,length.out=30)

y=predict(anes1,data.frame(conc=x),type='response')

plot(prop~conc,pch=16,col='red',data=anestot,xlim=c(0.5,3),main='Logistic回归曲线图',ylab='病人静止概率',xlab='麻醉剂量')

lines(y~x,lty=2,col='blue')

案例二：利用iris数据集，进行逻辑回归二分类测试，该数据集是R语言自带得数据集，包括四个属性，和三个分类。逻辑回归我们用glm函数实现，该函数提供了各种类型的回归，如：提供正态、指数、gamma、逆高斯、Poisson、二项。我们用的logistic回归使用的是二项分布族binomial。

index ->

将种类为setosa的数据排除出我们需要的数据集

ir ->

levels(ir$Species)[1] ->

生成训练集

split ->

ir_train ->

生成测试集

ir_test ->

通过训练集建立模型

model ->

summary(model)

模型运行结果：

Call:

glm(formula = Species ~ ., family = binomial(link = 'logit'),data = ir_train)

Deviance Residuals:

Min 1Q Median 3Q Max

-1.339e-04 -2.100e-08 2.100e-08 2.100e-08 1.059e-04

Coefficients:

Estimate Std. Error z value Pr(>z)

(Intercept) -1502.72 363247.01 -0.004 0.997

Sepal.Length 12.45 66482.13 0.000 1.000

Sepal.Width -285.61 95437.92 -0.003 0.998

Petal.Length 154.76 115968.97 0.001 0.999

Petal.Width 869.60 204513.80 0.004 0.997

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 9.0949e+01 on 65 degrees of freedom

Residual deviance: 4.0575e-08 on 61 degrees of freedom

AIC: 10

Number of Fisher Scoring iterations: 25

通过anova()函数 对模型进行方差分析

anova(model, test='Chisq')

方差分析如下：

Analysis of Deviance Table

Model: binomial, link: logit

Response: Species

Terms added sequentially (first to last)

Df Deviance Resid. Df Resid. Dev Pr(>Chi)

NULL 65 90.949

Sepal.Length 1 18.934 64 72.015 1.353e-05 ***

Sepal.Width 1 0.131 63 71.884 0.7176

Petal.Length 1 51.960 62 19.924 5.665e-13 ***

Petal.Width 1 19.924 61 0.000 8.058e-06 ***

---

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

下面通过McFadden R2指标进一步对模型进行分析

install.packages('pscl')

library(pscl)

pR2(model)

llh                         llhNull                G2 McFadden      r2ML          r2CU

-2.028752e-08 -4.547461e+01 9.094922e+01 1.000000e+00 7.479224e-01 1.000000e+00

为了得到分类结果，需要设定一个阈值p0——当p大于p0时，认为该样本的响应变量为1，否则为0。阈值大小对模型的预测效果有较大影响，需要进一步考虑。首先必须明确模型预测效果的评价指标。

求解训练模型的最佳阀值

对模型做出预测结果

model ->

pre1=predict(model,type='response')

将预测概率pre1和实际结果放在一个数据框中

data1=data.frame(prob=pre1,obs=ifelse(ir_train$Species=='virginica',1,0))

将预测概率按照从低到高排序

data1=data1[order(data1$prob),]

n=nrow(data1)

tpr=fpr=rep(0,n)

根据不同的临界值threshold来计算TPR和FPR，之后绘制成图

for (i in 1:n){

threshold=data1$prob[i]

tp=sum(data1$prob>threshold&data1$obs==1)

fp=sum(data1$prob>threshold&data1$obs==0)

tn=sum(data$prob

fn=sum(data$prob

tpr[i]=tp/(tp+fn) #真正率

fpr[i]=fp/(tn+fp) #假正率

}

plot(fpr,tpr,type='l')

abline(a=0,b=1)

下面通过pROC包自动标出最优临界点(0.506)

install.packages('pROC')

library(pROC)

modelroc1=roc(ifelse(ir_train$Species=='virginica',1,0),pre1)

plot(modelroc1,print.auc=TRUE,auc.polygon=TRUE,grid=c(0.1,0.2),grid.col=c('green','red'),max.auc.polygon=TRUE,auc.polygon.col='skyblue',print.thres=TRUE)

评估模型的预测效果

predict ->

predict.results 0.506,'virginica','versicolor')

misClasificError ->

print(paste('Accuracy',1-misClasificError))

[1] 'Accuracy 1'

最后一步，我们将通过画ROC曲线，并计算其AUC面积，作为评估二类分类效果的一个典型测量

install.packages('ROCR')

library(gplots)

library(ROCR)

p ->

p.results 0.5,1,0)

pr ->

prf ->

plot(prf)

auc ->

auc <>

0.9285714

auc

real ->

data.frame(real,predict)

res 0.5,'virginca','versicorlor'))

查看模型效果

plot(res)

基于R的案例

以下这个例子主要参考《Introduction to statistical learning with R》，强烈推荐这本书。尤其是回归部分，讲解的很透彻。

数据以Smarket为例，该数据包含2001-2005年标准普尔500指数1250个观测值，9个变量。9个变量分别为：年份，前5个交易日的收益率，前一个交易日和交易额(单位：billions)，今天的回报率和走势(上升或下降)。

df=read.csv('house.csv',header=TRUE)

str(df)

data.frame': 1250 obs. of  9 variables:

$ Year     : num  2001 2001 2001 2001 2001 ...

$ Lag1     : num  0.381 0.959 1.032 -0.623 0.614 ...

$ Lag2     : num  -0.192 0.381 0.959 1.032 -0.623 ...

$ Lag3     : num  -2.624 -0.192 0.381 0.959 1.032 ...

$ Lag4     : num  -1.055 -2.624 -0.192 0.381 0.959 ...

$ Lag5     : num  5.01 -1.055 -2.624 -0.192 0.381 ...

$ Volume   : num  1.19 1.3 1.41 1.28 1.21 ...

$ Today    : num  0.959 1.032 -0.623 0.614 0.213 ...

$ Direction: Factor w/ 2 levels 'Down','Up': 2 2 1 2 2 2 1 2 2 2 ...

#以前5个交易日的收益率和前一个工作日的交易额为自变量，今天的走势做因变量做Logistic回归；

glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,

data=Smarket,family='binomial')

summary(glm.fit)

Call:

glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +

Volume, family = 'binomial', data = Smarket)

Deviance Residuals:

Min      1Q  Median      3Q     Max

-1.446  -1.203   1.065   1.145   1.326

Coefficients:

Estimate Std. Error z value Pr(>|z|)

(Intercept) -0.126000   0.240736  -0.523    0.601

Lag1        -0.073074   0.050167  -1.457    0.145

Lag2        -0.042301   0.050086  -0.845    0.398

Lag3         0.011085   0.049939   0.222    0.824

Lag4         0.009359   0.049974   0.187    0.851

Lag5         0.010313   0.049511   0.208    0.835

Volume       0.135441   0.158360   0.855    0.392

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 1731.2  on 1249  degrees of freedom

Residual deviance: 1727.6  on 1243  degrees of freedom

AIC: 1741.6

Number of Fisher Scoring iterations: 3

有时候，预测股市就是这么不靠谱。很多自变量没通过验证，接下来我们基于AIC准则用逐步回归做一下变量筛选；

注：AIC一般公式为 AIC=2k-ln(L),其中k为参数的个数，L为估计模型最大化的似然函数值；

step(glm.fit)

Start:  AIC=1741.58

Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume

Df Deviance    AIC

- Lag4    1   1727.6 1739.6

- Lag5    1   1727.6 1739.6

- Lag3    1   1727.6 1739.6

- Lag2    1   1728.3 1740.3

- Volume  1   1728.3 1740.3

1727.6 1741.6

- Lag1    1   1729.7 1741.7

#此处略去一百字；

Direction ~ 1

Call:  glm(formula = Direction ~ 1, family = 'binomial', data = Smarket)

Coefficients:

(Intercept)

0.07363

Degrees of Freedom: 1249 Total (i.e. Null);  1249 Residual

Null Deviance:    1731

Residual Deviance: 1731AIC: 1733

这个结果足以让你崩溃，扔硬币都比这靠谱。原来不用任何变量预测是最靠谱的。

#模型的评估

glm.probs =predict(glm.fit,type ='response')

glm.probs[1:10]

#生成哑变量

contrasts(Smarket$Direction)

glm.pred=rep ('Down ' ,1250)

glm.pred[glm .probs >.5]=' Up'

table(glm .pred ,Direction )

mean(glm.pred== Direction )

[1] 0.5216

通过混淆矩阵计算分类的准确性仅有52%，很不乐观；

挖掘本质上是挖掘有趣的的模式，以上数据说明我们白忙活了一场，但是我们通过实例学习了Logistic模型。

当然我们还可以通过数据变换或构造新的变量来提升拟合降低AIC，最终，模型讲究的还是泛化能力；

有时候：拟合很丰满，泛化很骨感。->->->