广义线性回归模型之0,1变量回归（logit/probit回归）—R语言实现
1、广义线性回归
广义线性模型有三个组成部分：
(1) 随机部分, 即变量所属的指数族分布 族成员, 诸如正态分布, 二项分布, Poisson 分布等等. (2) 线性部分, 即 η = x⊤β. (3) 连接函数 g(µ) = η。 R 中的广义线性模型函数glm() 对指数族中某分布的默认连接函数 是其典则连接函数, 下表列出了 R 函数glm() 所用的某些指数族分布的 典则连接函数.
2、0-1因变量的回归模型
对于因变量为0,1变量的问题，可以考虑两种模型来解决
经过Probit变换和Logit变换，两种模型可以写成：
多变量情况：
logit回归
probit回归
3、案例分析
R中的广义线性回归函数为：glm()
语法为：glm(formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, offset, control = glm.control(…), model = TRUE, method = “glm.fit”, x = FALSE, y = TRUE, contrasts = NULL, …)
与线性回归lm不同之处就在于参数family，这个参数的作用在于定义一个族以及连接函数，使用该连接函数将因变量的期望与自变量联系起来，例如：上面讲到的logit模型和probit模型的参数分别为family= binomial(link=logit)和family= binomial(link=probit).
数据集介绍
数据来自新竹市输血服务中心的记录http://archive.ics.uci.edu/ml/datasets/Blood+Transfusion+Service+Center变量有Recency（上次献血距离研究时的月份），Frequency（总献血次数），Time（第一次献血是多少个月之前），Donate（是否将在2007年3月再献血，1为会，0为不会）
分别三个变量与因变量的相关性
通过图可以看出，三个变量与因变量还是有一定的关系的。
完整代码附上，备注的比较详细，也比较简单。
a<-read.csv("Trans.csv",header=T)
a<-a[,-3] #去掉第三列
names(a)<-c("x1","x2","x3","y") #设置变量名
a1<-a[c(1:400),] #取前400行数据，赋值给a1
a2=a[c(401:748),] #取出从401行开始剩下所有行的数据，赋值给a2
a1[c(1:5),] #展示a1前5行数据
#画图
par(mfrow=c(2,2)) #设置画图模式为2*2的格式
boxplot(x1~y,data=a1,main="Recency") #画出Recency与Donated的盒状图
boxplot(x2~y,data=a1,main="Frequency") #画出Frequency与Donated的盒状图
boxplot(x3~y,data=a1,main="Time") #画出Time与Donated的盒状图
par(mfrow=c(1,1)) #设置画图模式为1x1的格式
#拟合回归图
glm0.a=glm(y~1,family=binomial(link=logit),data=a1) #拟合logistic回归，不使用任何变量的空模型
glm1.a=glm(y~x1+x2+x3,family=binomial(link=logit),data=a1) #拟合logit回归,使用所有变量全模型
anova(glm0.a,glm1.a) #计算glm0.a与glm1.a的deviance
1-pchisq(30.13,3) #计算模型显著性检验的P值
glm1.b=glm(y~x1+x2+x3,family=binomial(link=probit),data=a1) #拟合porbit回归，使用所有变量的全模型
#参数估计
library(car)
Anova(glm1.a,type="III") #对模型glm1.a做三型方差分析
summary(glm1.a) #显示模型glm1.a的各方面细节，包括参数估计值、P值等
Anova(glm1.b,type="III") #对模型glm1.b做三型方差分析
summary(glm1.b) #显示模型glm1.b的各方面细节，包括参数估计值、P值等
AIC(glm1.a)
AIC(glm1.b)
#预测与评估
p=predict(glm1.a,a2) #利用模型glm1.a对数据a2进行预测
p=exp(p)/(1+exp(p)) #计算预测得到的概率
p1 = predict(glm1.b,a2)#标准正分布
p1 = pnorm(p1)
a2$y.pred=1*(p>0.5) #以0.5为阈值生成预测值
table(a2[,c(4,5)]) #计算预测值与真实值的2维频数表
a21=a2[a2$y==1,]
nrow(a21)
nrow(a2)
a2$y.pred=1*(p>0.18) #以0.18为阈值生成预测值
table(a2[,c(4,5)]) #计算预测值与真实值的2维频数表
ngrids=100 #设置格点数为100
TPR=rep(0,ngrids) #为TPR(true positive ratio)赋初值
FPR=rep(0,ngrids) #为FPR(false positive ratio)赋初值
for(i in 1:ngrids){
p0=i/ngrids; #选取阈值p0
y.true=a2$y #从a2中取出真实值并赋值给y.true
y.pred=1*(p>p0) #以p0为阈值生成预测值
TPR[i]=sum(y.pred*y.true)/sum(y.true) #计算TPR
FPR[i]=sum(y.pred*(1-y.true))/sum(1-y.true) #计算FPR
}
#画出ROC图
plot(FPR,TPR,type="l",col=2) #画出FPR与TPR的散点图，即ROC曲线
points(c(0,1),c(0,1),type="l",lty=2) #添加对角线
p=matrix(0,length(a2[,1]),3) #生成矩阵，用于存储各模型的预测值
p[,1]=predict(glm0.a,a2) #利用模型glm0.a对数据a2进行预测
p[,2]=predict(glm1.a,a2) #利用模型glm1.a对数据a2进行预测
p[,3]=predict(glm1.b,a2) #利用模型glm1.b对数据a2进行预测
p[,c(1:2)]=exp(p[,c(1:2)])/(1+exp(p[,c(1:2)])) #计算预测得到的概率
p[,3]=pnorm(p[,3]) #计算预测得到的概率
plot(c(0,1),c(0,1),type="l",main="FPR vs. TPR",xlab="FPR",ylab="TPR") #画图，生成基本框架
FPR=rep(0,ngrids) #为FPR赋初值
TPR=rep(0,ngrids) #为TPR赋初值
for(k in 1:3){
prob=p[,k] #ȡ??p?е?K?е?ֵ??????K??ģ?͵?Ԥ??????
for(i in 1:ngrids){
p0=i/ngrids #ѡȡ??ֵ
y.hat=1*(prob>p0) #??????ֵ????Ԥ??ֵ
FPR[i]=sum((1-y.true)*y.hat)/sum(1-y.true) #????FPR
TPR[i]=sum(y.true*y.hat)/sum(y.true) #????TPR
}
points(FPR,TPR,type="b",col=k,lty=k,pch=k) #??ͼ?????ӵ?k??ģ?͵?TPR??FPR??ɢ??ͼ
}
legend(0.6,0.5,c("logit??ģ??","logitȫģ??",
"probitȫģ??"),lty=c(1:3),col=c(1:3),pch=c(1:3))