R语言（五）——横截面数据分类：经典方法（logistic、probit、判别分析）

目录

一、数据

二、logistic回归

1.拟合

2.预测

三、probit回归

四、经典判别分析（线性、混合线性、灵活线性）

五、交叉验证与比较

---

一、数据

脊柱数据（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方法