data <- read.csv("I:/R/R2/kechengsheji.csv")    #读取数据
a1 <- data[1:50,]          #取前50行的数据建立模型
a2 <- data[51:60,]         #取后10行的数据进行预测
par(mfrow = c(2,2))        #设置画图格式为2*2
boxplot(sex~yonxhsx,data = a1,main ="sex")
boxplot(zy~yonxhsx,data = a1,main ="zy")
boxplot(ah~yonxhsx,data = a1,main ="ah")
boxplot(swsj~yonxhsx,data = a1,main ="swsj")
par(mfrow = c(2,2))
boxplot(yxsj~yonxhsx,data = a1,main ="yxsj")
boxplot(zyqk~yonxhsx,data = a1,main ="zyqk")
boxplot(sxcj~yonxhsx,data = a1,main ="sxcj")
glm0.a <- glm(yonxhsx~1,family = binomial(link = logit),data = a1)       #拟合logit回归，不使用任何变量的空模型
glm1.a <- glm(yonxhsx~sex+zy+ah+swsj+yxsj+zyqk+sxcj, family = binomial(link = logit),data = a1)  #logit回归全模型
anova(glm0.a,glm1.a)     #计算glm0.a和glm1.a的dexiance
1-pchisq(33.978,6)       #计算模型的显著性检验的P值
glm0.b <- glm(yonxhsx~1,family = binomial(link = probit),data = a1)      #拟合probit回归，不使用任何变量的空模型
glm1.b <- glm(yonxhsx~sex+zy+ah+swsj+yxsj+zyqk+sxcj, family = binomial(link = probit),data = a1)    #probit回归全模型
anova(glm0.b,glm1.b)    #计算glm0.b和glm1.b的dexiance
1-pchisq(33.834,6)      #计算模型的显著性检验的P值
library(car)            #载入程序包car
anova(glm1.a,type="lll")   #对模型glm1.a作三型方差分析
summary(glm1.a)   #显示模型glm1.a的各方面细节，包括参数估计值、P值等
anova(glm1.b,type="lll")    #对模型glm1.b作三型方差分析
summary(glm1.b)    #显示模型glm1.b的各方面细节，包括参数估计值、P值等
glm2.b <- glm(yonxhsx~zy+ah+zyqk, family = binomial(link = probit),data = a1)     #拟合简化后的probit 回归模型
summary(glm2.b)      #显示模型glm2.b的各方面细节，包括参数估计值、P值等
glm2.a <- glm(yonxhsx~zy+ah+zyqk, family = binomial(link = logit),data = a1)    #拟合简化后的logit 回归模型
summary(glm2.a)    #显示模型glm2.a的各方面细节，包括参数估计值、P值等
deviance(glm0.a)       #计算模型glm0.a的deviance
deviance(glm1.a)       #计算模型glm1.a的deviance
AIC(glm0.a)            #计算模型glm0.a的AIC取值
AIC(glm1.a)            #计算模型glm1.a的AIC取值
AIC(glm0.a,k = log(length(a1[,1])))       #计算模型glm0.a的BIC取值
AIC(glm1.a,k = log(length(a1[,1])))       #计算模型glm1.a的BIC取值
logit.aic <- step(glm1.a,trace=0)    #根据AIC准则选择最优模型，并赋值给logit.aic
summary(logit.aic)    #显示模型logit.aic的各方面细节，包括参数估计值、P值等
n <- length(a1[,1])    #样本大小
logit.bic <- step(glm1.a,k = log(n),trace=0)    #根据BIC准则选择最优模型，并赋值给logit.bic
summary(logit.bic)       #显示模型logit.bic的各方面细节，包括参数估计值、P值等
probit.aic = step(glm1.b,trace = 0)    #根据AIC准则选择最优模型，并赋值给probit.aic
summary(probit.aic)     #显示模型probit.aic的各方面细节，包括参数估计值、P值等
probit.bic = step(glm1.b,k = log(n),trace = 0)   #根据BIC准则选择最优模型，并赋值给probit.bic
summary(probit.bic)     #显示模型probit.bic的各方面细节，包括参数估计值、P值等
summary(glm1.a)     #显示模型glm1.a的各方面细节，包括参数估计值、P值等
p <- predict(glm2.a,a2)    #利用模型glm1.a对数据a2进行预测
p <- exp(p)/(1+exp(p))    #计算预测得到的概率
a2$yonxhsx.pred <- 1*(p>0.5)   #以0.5为阈值生成预测值
table(a2[,c(7,8)])           #计算预测值与真实值的2维数表
a2$yonxhsx.pred <- 1*(p>0)    #以0为阈值生成预测值
table(a2[,c(7,8)])           #计算预测值与真实值的2维数表
a2$yonxhsx.pred <- 1*(p>0.05)   #以0.05为阈值生成预测值
table(a2[,c(7,8)])         #计算预测值与真实值的2维数表
ngrids <- 100           #设置格点数为100
TPR <- rep(0,ngrids)    #为TPR（true positive rate）赋初值
FPR <- rep(0,ngrids)    #为FPR（false positive rate）赋初值
for(i in 1:ngrids){
        p0 = i/ngrids                #选取初值p0
        yonxhsx.true = a2$yonxhsx    #从a2中选出真实值并赋给yonxhsx.true
        yonxhsx.pred = 1*(p > p0)      #以0.05为阈值生成预测值
        FPR[i] = sum(yonxhsx.pred*(1-yonxhsx.true))/sum(1-yonxhsx.true)   #计算FPR
        TPR[i] = sum(yonxhsx.pred*yonxhsx.true)/sum(yonxhsx.true)  #计算TPR
}
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]),6)    #生成矩阵，用于存储各模型的预测值
p[,1] <- predict(glm1.a,a2)      #利用模型glm1.a对数据a2进行预测
p[,2] <- predict(logit.aic,a2)    #利用模型logit.aic对数据a2进行预测
p[,3] <- predict(logit.bic,a2)    #利用模型logit.bic对数据a2进行预测
p[,c(1:3)] <- exp(p[,c(1:3)])/(1+exp(p[,c(1:3)]))  #计算所得到的概率
p[,4] <- predict(glm1.b,a2)   #利用模型glm1.b对数据a2进行预测
p[,5] <- predict(logit.aic,a2)   #利用模型logit.aic对数据a2进行预测
p[,6] <- predict(logit.bic,a2)   #利用模型logit.bic对数据a2进行预测
p[,c(4:6)] <- exp(p[,c(4:6)])/(1+exp(p[,c(4:6)]))    #计算所得到的概率
plot(c(0,1),c(0,1),type = "l",main = "FPR vs TPR",xlab = "FPR",ylab = "TPR")      #画图，生成基本框架
TPR <- rep(0,ngrids)     #为TPR赋初值
FPR <- rep(0,ngrids)     #为FPR赋初值
for(k in 1:6){
        prob = p[,k]  #取出p中的k列的值，即第k个模型的预测概率
        for(i in 1:ngrids){
                p0 = i/ngrids         #选取阈值
                yonxhsx.hat = 1*(prob > p0)    #根据阈值生成预测值
                TPR[i] = sum(yonxhsx.true*yonxhsx.hat)/sum(yonxhsx.true)   #计算TPR
                FPR[i] = sum(yonxhsx.hat*(1-yonxhsx.true))/sum(1-yonxhsx.true)   #计算FPR
        }
        points(FPR,TPR,type = "b",col = k,lty = k,pch = k)   #向图上添加第k个模型的TPR与FPR的散点图
}
legend(0.6,0.5,c("LOGIT FULL MODEL","LOGIT AIC MODEL",
"LOGIT BIC MODEL","PROBIT FULL MODEL","PROBIT AIC MODEL",
"PROBIT BIC MODEL"),lty = c(1:6),col = c(1:6),pch = c(1:6)) #为6个模型添加标识，区分6个模型

好久了。这个问题好像还是没有解决。。。谢谢各位啦。

还没毕业么？

见楼下的回复。。。本楼重复了，请版主删掉？

问题在于仿真出来的y均值太大，导致它变成0-1变量时取0值的样本量太少。这样的y跟x交互之后的列联表长这样：

	y=0	y=1
x=0	54	6628
x=1	2	3316

可想而知，这样的数据分布对计算机的参数估计来说挑战非常大；即使x和y相关，计算机也算不出来。一般来说，每个单元格的占比不应少于总样本的5%。

解决问题的办法也比较简单，就是把y的均值调小，让x和y的分布更均衡些。具体见下面的code：