【统计模型】大学生恋爱数据分析报告

 a=read.csv("D:/个人成长/学业/课程/大三下课程/统计模型/作业/第三次作业/大学生恋爱数据.csv",header=T)##读入文件

    
 a[c(1:5),]
    
 attach(a)
    
 b=a[,22:32]
    
 JB1=as.factor(JB1)
    
 JB3=as.factor(JB3)
    
  
    
 #统计恋爱状况
    
 x=c(82,211)
    
 color=c('red','orange')
    
 piepercent1=round(100*x/sum(x),1)
    
 pie(x,labels=piepercent1,main="是否恋爱",col=color)
    
 legend("topright",c("否","是"),cex=0.8,fill=color)
    
  
    
 #学生基本情况统计——年级、性别、家乡、身高、体重
    
 x1=c(42,51,57,143)#年级
    
 piepercent2=round(100*x1/sum(x1),1)
    
 x2=c(144,149)#性别
    
 piepercent3=round(100*x2/sum(x2),1)
    
 x3=c(34,113,82,42,22)#家乡
    
 piepercent4=round(100*x3/sum(x3),1)
    
 par(mfrow=c(1,3))
    
 pie(x1,labels=piepercent2,main="年级",col=rainbow(length(x1)))
    
 legend("topright",c("大一","大二","大三","大四"),cex=1.5,fill=rainbow(length(x1)))
    
 pie(x2,labels=piepercent3,main="性别",col=color)
    
 legend("topright",c("男","女"),cex=1.5,fill=color)
    
 pie(x3,labels=piepercent4,main="家乡",col=rainbow(length(x3)))
    
 legend("topright",c("一线城市","二线城市","三线城市","县级市","农村"),cex=1.5,fill=rainbow(length(x3)))
    
 #学生身高体重统计
    
 par(mfrow=c(1,2))
    
 hist(a$身高,main="身高",xlab="组别" ,ylab = "频数")#直方图
    
 hist(a$体重,main="体重",xlab="组别" ,ylab = "频数")#直方图
    
  
    
 #学生恋爱史统计
    
 y1=c(134,159)#是否追求过别人
    
 piepercent11=round(100*y1/sum(y1),1)
    
 y2=c(75,218)#是否被别人追求过
    
 piepercent12=round(100*y2/sum(y2),1)
    
 y3=c(205,88)#寝室同学是否谈过恋爱
    
 piepercent13=round(100*y3/sum(y3),1)
    
 par(mfrow=c(1,3))
    
 pie(y1,labels=piepercent11,main="是否追求过别人",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(y2,labels=piepercent12,main="是否被别人追求过",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(y3,labels=piepercent13,main="寝室同学是否谈过恋爱",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
  
    
 #学生职务担任情况
    
 z1=c(101,192)#班干部
    
 piepercent21=round(100*z1/sum(z1),1)
    
 z2=c(239,54)#党员
    
 piepercent22=round(100*z2/sum(z2),1)
    
 par(mfrow=c(1,2))
    
 pie(z1,labels=piepercent21,main="班干部",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(z2,labels=piepercent22,main="党员",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
  
    
 #学生运动情况
    
 w1=c(257,36)#足球
    
 piepercent31=round(100*w1/sum(w1),1)
    
 w2=c(227,66)#篮球
    
 piepercent32=round(100*w2/sum(w2),1)
    
 w3=c(229,64)#乒乓球
    
 piepercent33=round(100*w3/sum(w3),1)
    
 w4=c(168,125)#羽毛球
    
 piepercent34=round(100*w4/sum(w4),1)
    
 w5=c(144,149)#跑步
    
 piepercent35=round(100*w5/sum(w5),1)
    
 w6=c(220,73)#台球
    
 piepercent36=round(100*w6/sum(w6),1)
    
 par(mfrow=c(2,3))
    
 pie(w1,labels=piepercent31,main="足球",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(w2,labels=piepercent32,main="篮球",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(w3,labels=piepercent33,main="乒乓球",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(w4,labels=piepercent34,main="羽毛球",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(w5,labels=piepercent35,main="跑步",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(w6,labels=piepercent36,main="台球",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
  
    
 #学生才艺情况
    
 v1=c(195,98)#唱歌
    
 piepercent41=round(100*v1/sum(v1),1)
    
 v2=c(270,23)#主持
    
 piepercent42=round(100*v2/sum(v2),1)
    
 v3=c(268,25)#舞蹈
    
 piepercent43=round(100*v3/sum(v3),1)
    
 v4=c(180,110)#乐器
    
 piepercent44=round(100*v4/sum(v4),1)
    
 v5=c(133,160)#其他才艺
    
 piepercent45=round(100*v5/sum(v5),1)
    
 par(mfrow=c(2,3))
    
 pie(v1,labels=piepercent41,main="唱歌",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(v2,labels=piepercent42,main="主持",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(v3,labels=piepercent43,main="舞蹈",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(v4,labels=piepercent44,main="乐器",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 pie(v5,labels=piepercent45,main="其他才艺",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
  
    
 #学生颜值情况
    
 t1=c(90,203)#是否戴眼镜
    
 piepercent51=round(100*t1/sum(t1),1)
    
 par(mfrow=c(1,2))
    
 pie(t1,labels=piepercent51,main="是否戴眼镜",col=color)
    
 legend("topright",c("否","是"),cex=1.5,fill=color)
    
 hist(a$颜值,main="颜值",xlab="组别" ,ylab = "频数")#直方图
    
  
    
 #学生生活规划情况
    
 par(mfrow=c(2,4))
    
 hist(a$每周自习时间,main="每周自习时间",xlab="组别" ,ylab = "频数")
    
 hist(a$每周娱乐时间,main="每周娱乐时间",xlab="组别" ,ylab = "频数")
    
 hist(a$每周睡觉时间,main="每周睡觉时间",xlab="组别" ,ylab = "频数")
    
 hist(a$每周运动时间,main="每周运动时间",xlab="组别" ,ylab = "频数")
    
 hist(a$每月话费,main="每月话费",xlab="组别" ,ylab = "频数")
    
 hist(a$学生组织个数,main="学生组织个数",xlab="组别" ,ylab = "频数")
    
 hist(a$成绩水平,main="成绩水平",xlab="组别" ,ylab = "频数")
    
 hist(a$生活费_百元,main="生活费_百元",xlab="组别" ,ylab = "频数")
    
  
    
 #连续变量之间相关性热力图
    
 library(corrplot)
    
 k=cor(b,use='everything',method='pearson')
    
 par(mfrow=c(1,1))
    
 corrplot(k,addCoef.col = "black")
    
  
    
 #绘制箱线图
    
 par(mfrow=c(2,4))
    
 boxplot(a$每周自习时间~a$是否恋爱,ylab="每周自习时间",xlab="是否恋爱",data=a,names=c("否","是"))
    
 boxplot(a$每周娱乐时间~a$是否恋爱,ylab="每周娱乐时间",xlab="是否恋爱",data=a,names=c("否","是"))
    
 boxplot(a$每周睡觉时间~a$是否恋爱,ylab="每周睡觉时间",xlab="是否恋爱",data=a,names=c("否","是"))
    
 boxplot(a$每周运动时间~a$是否恋爱,ylab="每周运动时间",xlab="是否恋爱",data=a,names=c("否","是"))
    
 boxplot(a$每月话费~a$是否恋爱,ylab="每月话费",xlab="是否恋爱",data=a,names=c("否","是"))
    
 boxplot(a$学生组织个数~a$是否恋爱,ylab="学生组织个数",xlab="是否恋爱",data=a,names=c("否","是"))
    
 boxplot(a$成绩水平~a$是否恋爱,ylab="成绩水平",xlab="是否恋爱",data=a,names=c("否","是"))
    
 boxplot(a$生活费_百元~a$是否恋爱,ylab="生活费_百元",xlab="是否恋爱",data=a,names=c("否","是"))
    
  
    
 #全模型
    
 model.full=glm(Y~as.factor(JB1)+JB2+as.factor(JB3)+JB4+JB5+LA1+LA2+LA3+ZW1+ZW2+YD1+YD2+YD3+YD4+YD5+YD6+CY1+
    
              CY2+CY3+CY4+CY5+YZ1+YZ2+GH1+GH2+GH3+GH4+GH5+GH6+GH7+GH8,family=binomial(link=logit),data=a)
    
 #模型结果，不显著的变量也要解读，加上不具有统计学意义
    
 summary(model.full)
    
 #似然比卡方检验模型整体效果
    
 1-pchisq(30.56,df=7)
    
  
    
 pred=predict(model.full,a)
    
  
    
 #基于AIC准则下变量的选择
    
 c(AIC(model.full),BIC(model.full))
    
 model.aic=step(model.full,trace = F)
    
 summary(model.aic)
    
 ss=length(a[,1])#样本量
    
 #基于BIC准则下变量的选择
    
 model.bic=step(model.full,trace = F,k=log(ss))
    
 summary(model.bic)
    
 #只留了特别显著的变量
    
  
    
 library(pROC)
    
 #多次模拟，去除随机误差的影响（了解即可）
    
 nsimu=100#进行100次模拟
    
 p=0.7#用作训练集的样本概率
    
 ss0=round(ss*p)#训练集样本量
    
 AUC=as.data.frame(matrix(0,nsimu,3))#100行，3列的零数据框
    
 names(AUC)=c("全模型","AIC模型","BIC模型")
    
 #开始模拟
    
 for(i in 1:nsimu){
    
   #打乱a样本顺序，随即编号并排序
    
   aa=a[order(runif(ss)),]
    
   #数据集aa的前70%作为训练集
    
   A0=aa[c(1:ss0),]
    
   #数据集aa的后30%作为测试集
    
   A1=aa[-c(1:ss0),]
    
   
    
   model.1=glm(Y~as.factor(JB1)+JB2+as.factor(JB3)+JB4+JB5+LA1+LA2+LA3+ZW1+ZW2+YD1+YD2+YD3+YD4+YD5+YD6+CY1+
    
             CY2+CY3+CY4+CY5+YZ1+YZ2+GH1+GH2+GH3+GH4+GH5+GH6+GH7+GH8,family=binomial(link=logit),data=A0)
    
   model.2=glm(Y~LA1+LA2+LA3+YD3+YD6+GH5+GH7,family=binomial(link=logit),data=A0)
    
   model.3=glm(Y~LA1+LA2+LA3+GH5,family=binomial(link=logit),data=A0)
    
   
    
   #测试集检验模型效果，计算预测值
    
   pred.1=predict(model.1,A1)
    
   pred.2=predict(model.2,A1)
    
   pred.3=predict(model.3,A1)
    
   
    
   #计算AUC值
    
   y=A1$Y
    
   auc.1=roc(y,pred.1)$auc
    
   auc.2=roc(y,pred.2)$auc
    
   auc.3=roc(y,pred.3)$auc
    
   
    
   #将各个AUC值填充到零矩阵中
    
   AUC[i,]=c(auc.1,auc.2,auc.3)
    
 }
    
 #绘制箱线图看三个模型的AUC分布情况，AUC越大，模型分类效果越好
    
 par(mfrow=c(1,1))
    
 boxplot(AUC,main="外样本AUC对比")
    
 #利用最后一次模拟数据绘制三个模型的ROC曲线
    
 #计算混淆矩阵
    
 roc.1=roc(y,pred.1)
    
 roc.2=roc(y,pred.2)
    
 roc.3=roc(y,pred.3)
    
 #绘制三条ROC曲线，比较效果
    
 par(mfrow=c(1,3))
    
 plot(roc.1,main="全模型")
    
 plot(roc.2,main="AIC模型")
    
 plot(roc.3,main="BIC模型")