目录
一、数据
脊柱数据(Column_2C.csv、Column_3C.csv)有两个版本,区别在于分为两类还是3类。数据可以从网站http://archive.ics.uci.edu/ml/datasets/Vertebral+Column下载,不过是.dat文件,需要进行相应的转换
参考我这篇博客Python将.dat文件转换成.csv文件。或者直接下载我上传的文件,是已经对格式和数据经过处理,可以直接进行分析的csv文件。
数据具有6个自变量(生物力学特征V1-V6,都是数量),和类别(V7)
两个数据内容分别为
类别 | 数量(人) | 代码 |
C2 | ||
正常 | 100 | NO |
不正常 | 210 | AB |
C3 | ||
正常 | 100 | NO |
椎间盘突出 | 60 | DH |
腰椎滑脱 | 150 | SL |
二、logistic回归
1.拟合
#logistic
#拟合
w2=read.csv("column_2C.csv")
a=glm(V7~.,w2,family="binomial") #默认logit连接函数
b=step(a) #做逐步回归筛选变量
summary(b) #结果汇总
在逐步回归过程中,变量3和4被淘汰,输出的AIC为188.95
2.预测
通过逐步回归或其他方法挑选的模型可能有助于解释模型,但不一定对预测有帮助。因此在我们分析的步骤中,还是利用逐步回归之前的a拟合模型
R自动把p(概率)作为排序靠后的类别的概率,于是运行代码levels(w2$V7)时,得到“AB”“NO”,p指的是NO的概率(R按照字典排序)
按照总误判率来寻找p值得阈值的函数如下:
BI=function(D,w,ff,fm="binomial"){
a=glm(ff,w,family=fm)
z=predict(a,w,type="response")
ee=NULL
for(p in seq(.1,.9,.01)){
u=rep(levels(w[,D])[2],nrow(w))
u[!(z>p)]=levels(w[,D])[1]
e=sum(u!=w[,D])/nrow(w)
#e=sum(u=="NO"&w[,D]=="AB")/sum(w[,D]=="AB")
ee=rbind(ee,c(p,e))
}
I=which(ee[,2]==min(ee[,2]))
P=ee[min(I),1]
return(ee[min(I),])
}
其中D代表因变量是第几个变量,w为使用数据,ff是拟合公式,fm是模型。在此函数中,p值从0.1到0.9按0.01的间隔搜索。调用函数
p.2C=BI(7,w2,V7~.)
得到阈值0.63,误判率0.1322。运行以下代码查看判别结果
#判别结果
D=7;p=p.2C[1];a=glm(V7~.,w2,family = "binomial")
z=(predict(a,w2,type="response")>p)
u=rep(levels(w2[,D])[2],nrow(w2))
u[!(z>p)]=levels(w2[,D])[1]
zz=table(w2[,7],u); #混淆矩阵
sum(diag(zz))/sum(zz) #正确率
混淆矩阵如图,正确率为0.871
还可以修改函数以改进效果,可以改变e的判别方式
三、probit回归
由于该数据的变量V6有出错信息,但我们又不愿意去掉该变量,因此代码中改用probit选项,此时不能计算AIC,也不能用逐步回归选择变量,但还是用BI()函数选择阈值
D=7;w=w2;ff=V7~.;fm=quasibinomial(link="probit")
pp.2C=BI(D,w2,ff,fm)
D=7;p=pp.2C[1];a=glm(V7~.,w2,family=fm)
summary(a)
z=(predict(a,w,type="response")>p)
u=rep(levels(w[,D])[2],nrow(w))
u[!(z>p)]=levels(w[,D])[1]
zz=table(w2[,7],u); #混淆矩阵
sum(diag(zz))/sum(zz) #正确率
得到的probit回归的估计系数及近似正态z检验的结果和分类展示如下,可知训练出来的模型:
四、经典判别分析(线性、混合线性、灵活线性)
三者代码基本相似,只是部分地方有微小差别
#线性判别分析
library(MASS)
w2=read.csv("column_2C.csv")
a=lda(V7~.,data=w2)
b=predict(a,w2)$class
zz=table(w2[,7],b) #混淆矩阵
sum(diag(zz))/sum(zz) #正确率
#混合线性判别分析
library(mda)
w2=read.csv("column_2C.csv")
set.seed(1010)
a=mda(V7~.,data=w2)
b=predict(a,w2)
zz=table(w2[,7],b) #混淆矩阵
sum(diag(zz))/sum(zz) #正确率
#灵活线性判别分析
library(splines)
library(Matrix)
library(fda)
w2=read.csv("column_2C.csv")
set.seed(1010)
a=fda(V7~.,data=w2)
b=predict(a,w2)
zz=table(w2[,7],b) #混淆矩阵
sum(diag(zz))/sum(zz) #正确率
五、交叉验证与比较
把数据随机分成Z份(此处Z=10),使用以下函数划分数据:
#产生数据集
Fold=function(Z=10,w,D,seed=7777){
n=nrow(w);d=1:n;dd=list()
e=levels(w[,D])
t=length(e)
set.seed(seed)
for(i in 1:t){
d0=d[w[,D]==e[i]]
j=length(d0)
ZT=rep(1:Z,ceiling(j/Z))[1:j]
id=cbind(sample(ZT,length(ZT)),d0)
dd[[i]]=id
}
#每个dd[[i]]是随机1:Z及i类的下标集组成的矩阵
mm=list()
for(i in 1:Z){
u=NULL
for(j in 1:t)u=c(u,dd[[j]][dd[[j]][,1]==i,2])
mm[[i]]=u #mm[[i]]为第i个下标集
}
return(mm) #输出Z个下标集
}
运行logistic、lda、mda、fda判别分析的交叉验证,代码如下:
BI=function(D,w,ff,fm="binomial"){
a=glm(ff,w,family=fm)
z=predict(a,w,type="response")
ee=NULL
for(p in seq(.1,.9,.01)){
u=rep(levels(w[,D])[2],nrow(w))
u[!(z>p)]=levels(w[,D])[1]
e=sum(u!=w[,D])/nrow(w)
#e=sum(u=="NO"&w[,D]=="AB")/sum(w[,D]=="AB")
ee=rbind(ee,c(p,e))
}
I=which(ee[,2]==min(ee[,2]))
P=ee[min(I),1]
return(ee[min(I),])
}
#读数据
w=read.csv("column_3C.csv")
D=7;Z=10
mm=Fold(Z,w,D,7777)
ff=paste(names(w)[D],"~.",sep="")
ff=as.formula(ff) #logistic回归用最优阈值
p=BI(D,w,ff)[1]
ERROR=matrix(0,Z,4)
J=1
for(i in 1:Z){
m=mm[[i]]
a=glm(V7~.,w,family="binomial")
z=(predict(a,w[m,],type="response")>p) #测试集预测
u=rep(levels(w[m,D])[2],length(m))
u[!z]=levels(w[m,D])[1]
ERROR[i,J]=sum(w[m,D]!=u)/length(m)
}
J=J+1
for(i in 1:Z){
m=mm[[i]]
a=lda(ff,w[-m,])
ERROR[i,J]=sum(w[m,D]!=predict(a,w[m,])$class)/length(m)
}
J=J+1
for(i in 1:Z){
m=mm[[i]]
set.seed(1010)
a=mda(ff,w[-m,])
ERROR[i,J]=sum(w[m,D]!=predict(a,w[m,]))/length(m)
}
J=J+1
for(i in 1:Z){
m=mm[[i]]
a=fda(ff,w[-m,])
ERROR[i,J]=sum(w[m,D]!=predict(a,w[m,]))/length(m)
}
ERROR=data.frame(ERROR)
names(ERROR)=c("logistic","lda","mda","fda")
#ERROR
ME=apply(ERROR,2,mean)
得到各种方法的误判率和平均误判率
三分类数据同上,只需改变读取的数据集。同时三分类不能使用logistic方法