c语言logistic函数,转载几个R语言中实现Logistic回归模型的案例

案例一：本文用例来自于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)