张金龙的博客分享 http://blog.sciencenet.cn/u/zjlcas 物种适应性、分布与进化

博文

广义线性模型中的逻辑斯蒂回归

已有 8178 次阅读 2017-10-19 21:06 |个人分类:科研笔记|系统分类:科研笔记|关键词:学者

其实题目不太准确,这篇短文介绍的是二项分布线性模型, 也就是当响应变量为01时的线性模型,例如物种是否出现,病人是否需要住院,学生是否成功考入大学,英语四六级成绩是否及格等。逻辑斯蒂回归模型,是二项分布线性模型中的一种,也是最广泛的一种。

事情还要从1935年说起。那一年,英国伦敦大学学院的C.I. Bliss 在《Annals of Applied Biology》上发表论文The Calculation of the Dosage-Mortality Curve这篇论文的表4,就是现在广泛用于逻辑斯蒂回归教学的beetle数据。该数据是C.I. Bliss从Strand在1930发表的论文中摘出来的。

Bliss的这篇论文主要以杂拟谷盗(Tribolium confusum)为例,系统介绍拟合存活比例曲线的方法。杂拟谷盗是一种粮食害虫。当时生物学家将481只杂拟谷盗分成8每一组分别在暴露在不同浓度的二硫化碳5小时,之后记录每个组有多少只甲虫死了由于甲虫对二硫化碳的反应只有两种可能,要么存活,要么死亡,所以这个数据集是典型的二项分布数据。

在这篇论文中,Bliss提出了一种连接函数,也就是事件发生与否,用转换以后的正态分布累计曲线去估算,称为Probit实际上,二项分布线性模型的连接函数有三种,分别为logit, probit以及cloglog。 logit、probit以及cloglog都是统计变换,目的是将二项分布数据转换为一种形式,以便能用线性模型估计直接估计模型的参数。其中最常见的为logit使用logit连接函数的二项分布回归称为逻辑斯蒂回归logit是事件发生与是否发生的比例 (odds)的对数值,形为log(y/(1-y))称为log oddscloglogcomplementary log log的缩写,是将响应变量进行两次log变换, 公式为log(-log(1-y)) cloglog常用来拟合数据出现大量极端数值的数据, 例如,大样地的幼苗的负密度制约研究中,部分胸径等级的树苗只有很少一部分死亡,cloglog能较准确捕捉到与树苗死亡相关性较强的环境因子, 如同种密度制约、光环境等。


下面用R代码演示如何进行二项分布回归,分别用logit, probit以及cloglog三种连接函数

library(broom)
####表格形式的数据
####tabular form
####数据来源www.stat.ualberta.ca/~kcarrier/STAT562/comp_log_log
####beetle
####   LOGDOSE     NUMBER_OF_BEETLES  NUMBER_KILLED  COMP_LOG_LOGPROBITLOGIT
####1   1.691                59             6          5.7   3.4 35.0
####2   1.724                60            13         11.3  10.7  9.8
####3   1.755                62           18         20.9   23.422.4
####4   1.784                56            28         30.3  33.8 33.9
####5   1.811                63            52         47.7  49.6 50.0
####6   1.837                59            53         54.2  53.4 53.3
####7   1.861                62            61         61.1  59.7 59.2
####8   1.884                60            60         59.9  59.2 58.8
####>
####二硫化碳的浓度
LOGDOSE          <-c(1.691, 1.724, 1.755, 1.784, 1.811, 1.837, 1.861, 1.884 )
####各小组杂拟谷盗的头数
NUMBER_OF_BEETLES<-c(59,    60,    62,    56,    63,    59,    62,    60    )    
####杂拟谷盗经过二硫化碳熏蒸5小时后,死亡的个体数
NUMBER_KILLED    <-c(6,     13,    18,    28,    52,    53,    61,    60    )
####合并成一个数据框
beetle<-data.frame(NUMBER_OF_BEETLES,NUMBER_KILLED, LOGDOSE)
####第一种方法,数据表格式
####调用glm函数,注意连接函数的选择
mod_logit  <-glm(cbind(NUMBER_KILLED, NUMBER_OF_BEETLES -NUMBER_KILLED)~LOGDOSE,data=beetle,  family=binomial(link="logit") )
mod_probit <-glm(cbind(NUMBER_KILLED, NUMBER_OF_BEETLES -NUMBER_KILLED)~LOGDOSE,data=beetle,  family=binomial(link="probit"))
mod_cloglog<-glm(cbind(NUMBER_KILLED, NUMBER_OF_BEETLES -NUMBER_KILLED)~LOGDOSE,data=beetle,  family=binomial(link="cloglog"))
####Showing significance of the variables
tidy(mod_logit )

##         term  estimate std.error statistic      p.value
##1 (Intercept) -60.74013  5.181870 -11.72166 9.871561e-32
##2     LOGDOSE  34.28593  2.913213  11.76911 5.631446e-32

tidy(mod_probit)

##         term  estimate std.error statistic      p.value
##1 (Intercept) -34.95614  2.648984 -13.19605 9.245343e-40
##2     LOGDOSE  19.74103  1.488046  13.26640 3.625628e-40

tidy(mod_cloglog)

##         term  estimate std.error statistic      p.value
##1 (Intercept) -39.52224  3.235649 -12.21463 2.596958e-34
##2     LOGDOSE  22.01474  1.796996  12.25085 1.662370e-34

####Showing AIC et al.
glance(mod_logit )

##  null.deviance df.null    logLik      AIC     BIC deviancedf.residual
##1      284.2024       7 -18.65681 41.31361 41.4725 11.11558          6

glance(mod_probit)

##  null.deviance df.null   logLik      AIC      BIC deviancedf.residual
##1      284.2024       7 -18.0925 40.18499 40.34388 9.986957          6

glance(mod_cloglog)

##  null.deviance df.null    logLik      AIC      BIC deviancedf.residual
##1      284.2024       7 -14.85619 33.71237 33.87126 3.514334          6

####绘图Visualisingthe results
mortality_rate<-beetle$NUMBER_KILLED/beetle$NUMBER_OF_BEETLES
dose<-beetle$LOGDOSE
plot(mortality_rate~dose,xlab="Dose of Carbon Disulfide",ylab="MortalityRate of Beetles",main="BinomialRegression")
LOGDOSE_new=seq(min(dose),max(dose),length.out=200)
mod_logit_line   <-predict(mod_logit,newdata=data.frame(LOGDOSE=LOGDOSE_new), type="response")
lines(mod_logit_line~LOGDOSE_new,col="blue")
mod_probit_line  <-predict(mod_probit,newdata=data.frame(LOGDOSE=LOGDOSE_new), type="response")
lines(mod_probit_line~LOGDOSE_new,col="red")
mod_cloglog_line <-predict(mod_cloglog,newdata=data.frame(LOGDOSE=LOGDOSE_new), type="response")
lines(mod_cloglog_line~LOGDOSE_new,col="green")
legend("bottomright",title="Linkfunctions",legend=c("logit","probit","cloglog"),lty=c(1,1,1),col=c("blue","red","green"))



#####第二种方法,每个个体生存状态的长数据格式。glm是可以接受这种格式的。
#####按照个体来计
dose_surived<-
c(rep(1.691,59-6),        
rep(1.724,60-13),        
rep(1.755,62-18),        
rep(1.784,56-28),        
rep(1.811,63-52),        
rep(1.837,59-53),        
rep(1.861,62-61),        
rep(1.884,60-60))
dose_dead<-c(
rep(1.691, 6),
rep(1.724,13),
rep(1.755,18),
rep(1.784,28),
rep(1.811,52),
rep(1.837,53),
rep(1.861,61),
rep(1.884,60))
status_survived<-c(
rep("alive",59-6),
rep("alive",60-13),
rep("alive",62-18),
rep("alive",56-28),
rep("alive",63-52),
rep("alive",59-53),
rep("alive",62-61),
rep("alive",60-60))
status_dead<-c(
rep("dead", 6),
rep("dead",13),
rep("dead",18),
rep("dead",28),
rep("dead",52),
rep("dead",53),
rep("dead",61),
rep("dead",60)
)

dose_long<-c(dose_dead,dose_surived)
status   <-as.factor(c(status_dead,status_survived))

incidence_mod_logit  <-glm(status~dose_long,family=binomial(link="logit"))
incidence_mod_probit <-glm(status~dose_long,family=binomial(link="probit"))
incidence_mod_cloglog<-glm(status~dose_long,family=binomial(link="cloglog"))

####Showing significance of the variables
tidy(incidence_mod_logit )

##         term  estimate std.error statistic      p.value
##1 (Intercept) -60.74013  5.181834 -11.72174 9.861985e-32
##2   dose_long  34.28593  2.913192  11.76920 5.625873e-32

tidy(incidence_mod_probit)

##         term  estimate std.error statistic      p.value
##1 (Intercept) -34.95613  2.649010 -13.19592 9.261453e-40
##2   dose_long  19.74102  1.488061  13.26627 3.631968e-40

tidy(incidence_mod_cloglog)

##         term  estimate std.error statistic      p.value
##1 (Intercept) -39.52224  3.235632 -12.21469 2.594819e-34
##2   dose_long  22.01474  1.796986  12.25092 1.660982e-34

####Showing AIC et al.
glance(incidence_mod_logit )

##  null.deviance df.null    logLik      AIC      BIC deviancedf.residual
##1       645.441     480 -186.1771 376.3542 384.7059 372.3542        479

glance(incidence_mod_probit)

##  null.deviance df.null    logLik      AIC      BIC deviancedf.residual
##1       645.441     480 -185.6128 375.2255 383.5773 371.2255        479

glance(incidence_mod_cloglog)

##  null.deviance df.null    logLik      AIC      BIC deviancedf.residual
##1       645.441     480 -182.3765 368.7529 377.1046 364.7529        479




https://m.sciencenet.cn/blog-255662-1081614.html

上一篇:大型森林样地样方的命名与坐标转换
下一篇:繁花满树 (十六首)

0

该博文允许注册用户评论 请点击登录 评论 (3 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2022-10-7 01:59

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部