脑中风预测

通过对一组人员的脑中风情况进行分析,探索影响脑中风的因素
用了逻辑回归,决策树,朴素贝叶斯
数据来于和鲸社区,上传数据的人是从kaggle下的

一、数据清洗

A<-read.csv("E:/R语言练习/脑中风预测/healthcare-dataset-stroke-data.csv",head=T);   #读取数据

查看每个字段的数据类型,方便起见在后面写上说明

mode(A[1,1])      # id                      人员编号

[1] “numeric”

mode(A[1,2])		#gender 				    性别

[1] “character”

mode(A[1,3])		#age 					年龄

[1] “numeric”

mode(A[1,4])		#hypertension 			0:无高血压,1:有高血压

[1] “numeric”

mode(A[1,5])		#heart_disease 			0:无心脏病,1:有心脏病

[1] “numeric”

mode(A[1,6])		#ever_married 			是否婚配

[1] “character”

mode(A[1,7])		#work_type 				工作类型

[1] “character”

mode(A[1,8])		#Residence_type 			居住类型

[1] “character”

mode(A[1,9])		#avg_glucose_level		血液中平均葡萄糖水平

[1] “numeric”

mode(A[1,10])		#bmi  					体重指数

[1] “character”

mode(A[1,11])		#smoking_status 			吸烟状况

[1] “character”

mode(A[1,12])		#stroke					0:不中风,1:中风

[1] “numeric”

dim(A) #查看数据量

[1] 5110 12
一共有5110条记录,每条记录有12个字段。

#性别处理:将男性改为1,女性改为0:

table(A$gender)   

Female Male Other
2994 2115 1
#可以看到有第三性别,不过只有一例,这里设为女性,因为女性数量较男性多,这样不至于使误差过大。

A[,2][A[,2]=="Male"]=1;
A[,2][A[,2]=="Female"|A[,2]=="Other"]=0;
A[,2]<-as.numeric(A[,2])

#婚配状况处理:将已婚配改为1,未婚配改为0:

table(A[,6]);

No Yes
1757 3353

A[,6][A[,6]=="Yes"]=1;
A[,6][A[,6]=="No"]=0;
A[,6]<-as.numeric(A[,6]);

#工作类型处理

table(A[,7]);

children Govt_job Never_worked Private Self-employed
687 657 22 2925 819

A[,7][A[,7]=="children"]=1;
A[,7][A[,7]=="Govt_job"]=2;
A[,7][A[,7]=="Never_worked"]=3;
A[,7][A[,7]=="Private"]=4;
A[,7][A[,7]=="Self-employed"]=5;
A[,7]<-as.numeric(A[,7]);

#居住类型处理:城镇改为1,农村改为0

table(A[,8]);

Rural Urban
2514 2596

A[,8][A[,8]=="Rural"]=0;
A[,8][A[,8]=="Urban"]=1;

#体重指数处理:

A[,10]<-as.numeric(A[,10]);

Warning message:
强制改变过程中产生了NA

sum(is.na(A1[,10]))

[1] 201
#可以看到一共有201个空值,201:5110约为1:25,不是个小数目。

sum(A1$gender=="Male")/sum(A1$gender=="Female") #缺失值的男女比例

[1] 1.072165

sum(A$gender=="Male")/sum(A$gender=="Female")  #整体的男女比例

[1] 0.7064128
#两者相差0.3,为了避免删除缺失记录后,样本不能准确反映中风与性别的关系,此处选择填充数据。填充需要对一些变量的数值分类,一会再填充。

#吸烟状况处理

table(A1[,11]);
A1[,11][A1[,11]=="formerly smoked"]=0;
A1[,11][A1[,11]=="never smoked"]=1;
A1[,11][A1[,11]=="smokes"]=2;
A1[,11][A1[,11]=="Unknown"]=3;
A1[,11]<-as.numeric(A1[,11]);
table(A1[,11]);

#填充bmi
#填充过程要先对连续型变量分类,再将对象的各个特征所对应的平均bmi算出,再将这些特征的bmi求一个平均数作为对象的bmi

A2<-A1; #接下来工作都在A2中执行,不破坏A1

#现将年龄分类:

hist(A2$age,breaks=10,labels=T)
 A2$age[A2$age<=10]=1;
 A2$age[A2$age>10&A2$age<=20]=2;
 A2$age[A2$age>20&A2$age<=30]=3;
 A2$age[A2$age>30&A2$age<=40]=4;
 A2$age[A2$age>40&A2$age<=50]=5;
 A2$age[A2$age>50&A2$age<=60]=6;
 A2$age[A2$age>60&A2$age<=70]=7;
 A2$age[A2$age>70]=8;

#现将平均葡萄糖水平分类:

hist(A2[,9],labels=T)
A2[,9][A2[,9]<=60]=1;
A2[,9][A2[,9]>60&A2[,9]<=80]=2;
A2[,9][A2[,9]>80&A2[,9]<=100]=3;
A2[,9][A2[,9]>100&A2[,9]<=120]=4;
A2[,9][A2[,9]>120&A2[,9]<=140]=5;
A2[,9][A2[,9]>140&A2[,9]<=160]=6;
A2[,9][A2[,9]>160&A2[,9]<=180]=7;
A2[,9][A2[,9]>180&A2[,9]<=200]=8;
A2[,9][A2[,9]>200&A2[,9]<=220]=9;
A2[,9][A2[,9]>220&A2[,9]<=240]=10;
A2[,9][A2[,9]>240]=11;

#将空值与非空值分开:

B1<-subset(A2,is.na(A2[,10]))
B2<-subset(A2,!is.na(A2[,10]))

#计算各个特征在各个水平下的平均值

C<-as.numeric(tapply(B2[,10],B2[,2],mean))
D<-as.numeric(tapply(B2[,10],B2[,3],mean))
E<-as.numeric(tapply(B2[,10],B2[,4],mean))
F<-as.numeric(tapply(B2[,10],B2[,5],mean))
G<-as.numeric(tapply(B2[,10],B2[,6],mean))
H<-as.numeric(tapply(B2[,10],B2[,7],mean))
I<-as.numeric(tapply(B2[,10],B2[,8],mean))
J<-as.numeric(tapply(B2[,10],B2[,9],mean))
K<-as.numeric(tapply(B2[,10],B2[,11],mean))
L<-as.numeric(tapply(B2[,10],B2[,12],mean))

#循环填充

for(i in 1:nrow(B1)){
B1[i,10]=(C[B1[i,2]+1]+D[B1[i,3]]+E[B1[i,4]+1]+F[B1[i,5]+1]+G[B1[i,6]+1]
+H[B1[i,7]]+I[B1[i,8]+1]+J[B1[i,9]]+K[B1[i,11]+1]+L[B1[i,12]+1])/10
}   

#计算各个空值并赋值给B1相对应的位置

AA1<-subset(A1,is.na(A1[,10]));
AA2<-subset(A1,!is.na(A1[,10]));
AA1[,10]<-B1[,10];  #对应位置赋值
AA1[(nrow(AA1)+1):nrow(A2),]<-AA2;
A3<-AA1;

#正式填充完毕!
#A3现在是填充过的尚未对连续变量进行分类的数据集

二、数据分析:

描述性分析:

A4<-subset(A3,A3[,12]==1); #将中风的人赋值给A4
A5<-subset(A3,A3[,12]==0); #将不中风的人赋值给A5

先查看总体中风情况:

hist(A3[,12],labels=T);  #直方图
AA3<-table(A3[,12])
per<-paste(names(AA3),"___",round(100*(AA3)/sum(AA3),2),"%")
pie(table(A3[,12]),labels=per);  #饼图

在这里插入图片描述
在这里插入图片描述
从中可以看到中风人数占调查总人数4.87%,共有249人,不中风占95.13%,有4861人。

#性别:

hist(A3[,2],col='red',labels=T);   #所有人按性别分类的直方图
hist(A4[,2],col='green',add=T,labels=T);   #中风者按性别分类的直方图
table(A4[,2])/table(A3[,2]); #男女各自中风率

在这里插入图片描述
在这里插入图片描述
可以看到一共调查了2995名女性 2115名男性,其中女性有141人中风,男性有108人中风,中风率分别为0.047和0.051,相差不大。

#按年龄分析:

hist(A3[,3],col='red',labels=T);   #所有人按年龄分类的直方图
hist(A4[,3],col='green',add=T,labels=T);   #中风者按年龄分类的直方图

在这里插入图片描述
明显看到年龄越大,中风的人数越多,50~80岁是中风高发段,中风率也基本呈增高趋势。
#高血压分析:

hist(A3[,4],col='red',labels=T);   #所有人按是否高血压分类的直方图
hist(A4[,4],col='green',add=T,labels=T);   #中风者按高血压分类的直方图
table(A4[,4])/table(A3[,4]); #各自中风率

在这里插入图片描述
在这里插入图片描述
可以看到高血压者有498人,其中中风的有66人,占比13.3%
血压正常的有4612人,其中中风的有183人,占比3.97%,明显低于高血压者
#心脏病分类:

hist(A3[,5],col='red',labels=T);   #所有人按是否心脏病分类的直方图
hist(A4[,5],col='green',add=T,labels=T);   #中风者按心脏病分类的直方图
table(A4[,5])/table(A3[,5]); #各自中风率

在这里插入图片描述
在这里插入图片描述
患有心脏病的有276人,其中有47人中风,占比17%
心脏正常的有4834人,其中有202人中风,占比4.2%,明显低于心脏病患者。
#是否结婚:

hist(A3[,6],col='red',labels=T);   #所有人按是否结婚分类的直方图
hist(A4[,6],col='green',add=T,labels=T);   #中风者按是否结婚分类的直方图
table(A4[,6])/table(A3[,6]); #各自中风率

在这里插入图片描述
在这里插入图片描述
可以看到有1757人没结婚,其中有29人中风,占比1.7%
有3353人结婚了,其中有220人中风,占比6.6%,明显高于未结婚的。
#工作类型:

hist(A3[,7],col='red',labels=T);   #所有人按工作类型分类的直方图
hist(A4[,7],col='green',add=T,labels=T);   #中风者按工作类型分类的直方图
table(A4[,7])/table(A3[,7]); #各自中风率

在这里插入图片描述
可以看到各个工作类型分别有687,657,22,2925,819人,中风者有2,33,0,149,65人。求中风率时报错
在这里插入图片描述
调整后:
在这里插入图片描述
调整后得到中风率分别如上图所示,可以看到第一类(小孩)与第三类工作(从不工作)明显低于另外三种,第五类(自由职业者)中风率最高接近8%,行政者和民营企业都是5%
#居住类型:

hist(A3[,8],col='red',labels=T);   #所有人按居住类型分类的直方图
hist(A4[,8],col='green',add=T,labels=T);   #中风者按居住类型分类的直方图
table(A4[,8])/table(A3[,8]); #各自中风率

在这里插入图片描述
在这里插入图片描述
可以看到各占一半,中风率也差不多,所以应该不是主要因素。
#平均葡萄糖水平

hist(A3[,9],col='red',labels=T);   #所有人按平均葡萄糖水平分类的直方图
hist(A4[,9],col='green',add=T,labels=T);   #中风者按平均葡萄糖水平分类的直方图

在这里插入图片描述
可以看到大部分人葡萄糖都在50150之间。葡萄糖为170250时的中风率比较高。
#体重指数

hist(A3[,10],col='red',labels=T);   #所有人按体重指数分类的直方图
hist(A4[,10],col='green',add=T,labels=T);   #中风者按体重指数分类的直方图

在这里插入图片描述
各自中风比例都在3%~5%左右,没有明显差异
#是否吸烟
在这里插入图片描述
在这里插入图片描述
最后一类是不知道吸烟的情况,所以不参与讨论。前面的三类中风率基本在5~8,变化不够明显。

逻辑回归:
先看看自变量间的线性关系:
在这里插入图片描述
可以看到是否结婚和年龄之间相关系数最高,为0.68,还没达到0.8,不强,但相关系数只能检测线性关系,如果不是线性的就不行了。很多自变量只取两个值,用相关系数来看相关性是不太合适的,但以我目前的水平还不会检测这个类型自变量和连续型的自变量的关系0.0
查看一下连续型变量之间的关系,看看能不能转化成线性的:

plot(A3[,3],A3[,9],main="年龄和平均葡萄糖水平");
plot(A3[,3],A3[,10],main="年龄和体重指数");
plot(A3[,10],A3[,9],main="体重指数和平均葡萄糖水平");

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
从第一张图可以看出年龄越大,葡萄糖水平比较高的人数越多,但总体还是葡萄糖低的人比较多,变化不明显,看不出有什么其他关系,两个变量间相关系数是0.24,基本可以认为没有线性关系。
从第二张图可以看出0~20岁时,两者逐渐增高,20岁之后就稳定了,20岁后基本没有关系,相关系数为0.33,20岁之前的数据贡献比较大,呈弱相关。
第三张图可以看出不管体重系数是多少,葡萄糖水平高的人都是少数,相关系数是0.17,基本没线性关系。
综上所述没有必要在这三个变量间做合并。
从描述性分析看,工作类型有五类,1,3类是一个水平,2,4是一个水平,5是一个水平,所以把1,3归为一类,24归为1类,5归为1类

A3[,7][A3[,7]==3]=1
A3[,7][A3[,7]==4]=2
A3[,7][A3[,7]==5]=3
table(A3[,7])

在这里插入图片描述

glm.sol<-glm(A3[,12]~A3[,2]+A3[,3]+A3[,4]+A3[,5]+A3[,6]+A3[,7]+
A3[,8]+A3[,9]+A3[,10]+A3[,11], family=binomial, data=A3)
summary(glm.sol)

逻辑回归结果:
在这里插入图片描述
性别的系数没通过检验,即是否中风对性别不明,符合描述性统计
年龄的系数通过了检验,而且大于0,即年龄越大,中风率越高,符合描述性统计
高血压的系数通过了检验,大于0,即血压越高,中风率越高,符合描述性统计
心脏病的系数未通过检验,不符合描述性统计
是否结婚的系数也没有通过检验,不符合描述性统计
工作类型的系数通过检验,符合描述性统计
居住类型的系数没通过检验,符合描述性统计
葡萄糖水平的系数通过检验,符合描述性统计
体重系数的系数没通过检验,符合描述性统计
是否吸烟的系数没通过检验,符合描述性统计
以上结果中有些量不符合描述性统计,先对其进行逐步回归:

step(glm.sol)

在这里插入图片描述
逐步回归结果建议选取年龄,血压,心脏病,工作类型,葡萄糖为自变量。这与我们描述性分析结果相同

glm.sol<-glm(A3[,12]~A3[,3]+A3[,4]+A3[,5]+A3[,7]+
A3[,9], family=binomial, data=A3)
summary(glm.sol)
kappa(glm.sol)

在这里插入图片描述
可以看到心脏病的系数还是不够好,可能是因为样本中患有心脏病的人比例太少了,只有276人,而且是两点分布,而总人数有5110人。其他的指标,年龄是每个水平都有很多人,高血压一共有接近500人,工作类型水平最少的也有819人,葡萄糖含量水平较高的也有500多人
算一下kappa值查看共线性情况
在这里插入图片描述
不算严重共线性
试试其他的模型,因为因变量是离散的,以我目前的水平,没法通过画图来看,有可能用别的模型拟合效果更好。
分别对每个变量试了指数函数,对数函数,二次多项式,发现有些虽然AIC降了一点点,但是没有一个新的变量的系数没有通过检验,所以不采纳。

结论:
令a=exp(0.07age+0.4hypertension+0.32heart_disease-0.36work_type+0.001*avg_glucose_level-6.89)
P=a/(1+a), 其中P为得脑中风的几率
决策树:

library(rpart)
BB1<-rbind(A3[1:200,1:11],A3[4201:5110,1:11]);
BB<-A3[201:4200,1:12];

#以前201至4200行数据作为训练数据
#其他行作为测试数据,因为A3的顺序和源数据是不一样的,所以这么分类使得两边中风比例都在0.2左右

fit<-rpart(stroke~gender+age+hypertension+heart_disease+ever_married+work_type+Residence_type+avg_glucose_level+bmi+smoking_status,method="class",cp=0.0001,minsplit=40,data=BB)
plot(fit,uniform=TRUE,main="Classification Tree for 脑中风");
text(fit,use.n=TRUE,all=TRUE);
result <- predict(fit,BB1,type="class")
source("D:/R软件1/R-4.0.4/bin/i386/数数.R")
count_result(result,A3[101:2000,12])#直接运行,准确率为

在这里插入图片描述
剔除变量:

fit<-rpart(stroke~age+hypertension+heart_disease+work_type+avg_glucose_level,method="class",cp=0.001,minsplit=70,data=BB)
plot(fit,uniform=TRUE,main="Classification Tree for 脑中风");
text(fit,use.n=TRUE,all=TRUE);
result <- predict(fit,BB1,type="class")
source("D:/R软件1/R-4.0.4/bin/i386/数数.R")
count_result(result,A3[101:2000,12])

准确率为
在这里插入图片描述
增大了0.003。

朴素贝叶斯:
用朴素贝叶斯是因为变量弱相关,它的表现也可以,这里来试试效果

library(e1071)
abc<-naiveBayes(stroke~gender+age+hypertension+heart_disease+ever_married+work_type+Residence_type+avg_glucose_level+bmi+smoking_status,data=BB)
result <- predict(abc,BB1,type="class");
count_result(result,A3[101:2000,12]);

这是没有剔除变量的
在这里插入图片描述

abc<-naiveBayes(stroke~age+hypertension+heart_disease+work_type+avg_glucose_level,method="class",data=BB)
result <- predict(abc,BB1,type="class");
count_result(result,A3[101:2000,12]);

在这里插入图片描述
剔除变量后准确率上升了0.004左右

  • 2
    点赞
  • 47
    收藏
    觉得还不错? 一键收藏
  • 7
    评论
评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值