代码无所谓高端和低端，能解决问题就行。
有序多分类Logistic回归介绍
二分类结局变量，如事件是否发生、是否死亡等，我们可以用二分类logistics回归分析结局变量和自变量因素之间的关系，结果也非常容易解释。但对于有序多分类结局指标，如病情的严重程度、对服务的评价好坏、患者满意度、药物的主观疗效好坏等，这样的资料往往写成如轻、中、重；差、一般、好等，我们就要另辟道路完成变量间的分析。就会用到有序多分类Logistic回归（Ordinal Logistic Regression）。
目录
第四步，批量建立公式，为了使用anova函数求P，所以这里建立空模型
网上的许多关于有序多分类Logistic回归写的都非常好,但是,他们都没能回答一个实际操作过程中不可避免的问题：缺失值。我们来看一组典型数据：
我们要回答的问题是各个基因分别和临床因素有无关系？显著吗？相关性如何？
那么这里就有：
黑色字体区间内是基因，也就是X变量（自变量）
红色字体区间内是临床因素，也就是Y变量（应变量）
我们要回答的是X和Y两个集合内的因素交叉配对建模，同时解决缺失值问题。
这里我们使用函数polr 来自{MASS}包
polr，全称：Ordered Logistic or Probit Regression
Fits a logistic or probit regression model to an ordered factor response. The default logistic case is proportional odds logistic regression, after which the function is named.
Usage
polr(formula, data, weights, start, ..., subset, na.action,
contrasts = NULL, Hess = FALSE, model = TRUE,
method = c("logistic", "probit", "loglog", "cloglog", "cauchit"))z
这里我们写作做简单的polr(Y~ +X, data=data)但是我发现，很多数据中，不少变量是带有缺失值的，那么在建立模型的时候，如果默认删除缺失值，就会有变量数量不一样的问题而报错。
接下来我们主要分几步完成这个任务。
上图只是举个例子，所以有部分变量删除了，我们的数据其实是这样子的
他有420个样本，9个基因，6个结局变量，且变了带有缺失值。
第一步，观察数据，做到心中有数
library(tidyverse)
library(MASS)
setwd("你的路径")
data<- read_csv("data_mullogistic.csv")
names(data)
#这里x变量为：
"hsa.miR.17.3p"
[3] "hsa.miR.17.5p" "hsa.miR.18a.3p"
[5] "hsa.miR.18a.5p" "hsa.miR.19a.3p"
[7] "hsa.miR.19a.5p" "hsa.miR.20a.3p"
[9] "hsa.miR.20a.5p" "MIR17HG"
#一共9个
#这里Y变量为：
[11] "fibrosis_ishak_score" "adjacent_hepatic_tissue_inflammation_extent_type"
[13] "Virus_infection" "neoplasm_histologic_grade"
[15] "pathologic_stage" "vascular_tumor_cell_type"
#一共6个心中知道一共有54个模型需要建立。
第二步，观察6个分类变量到底有哪些
分类变量的字符串往往不规则，且尚未排序，所以我们还要先观察6个分类变量到底有哪些，该如何排序。
#显示变量的种类，准备排序
#Y变量1："fibrosis_ishak_score"
table(data$fibrosis_ishak_score)
na.omit(unique(data$fibrosis_ishak_score))
#Y变量2："adjacent_hepatic_tissue_inflammation_extent_type"
table(data$adjacent_hepatic_tissue_inflammation_extent_type)
na.omit(unique(data$adjacent_hepatic_tissue_inflammation_extent_type))
#Y变量3："Virus_infection"
na.omit(unique(data$Virus_infection))
#Y变量4："neoplasm_histologic_grade"
na.omit(unique(data$neoplasm_histologic_grade))
#Y变量5："pathologic_stage"
na.omit(unique(data$pathologic_stage))
#Y变量6："vascular_tumor_cell_type"
na.omit(unique(data$vascular_tumor_cell_type ))第三步，手动使用factor将这些变量标签排序
这里最重要的是耐心和仔细。
#将定型资料排序
data$fibrosis_ishak_score<-factor( data$fibrosis_ishak_score,
levels = c("0_NoFibrosis","1,2_PortalFibrosis","3,4_FibrousSpeta","5_NodularFormationandIncompleteCirrhosis", "6_EstablishedCirrhosis" ),
labels = c("I","II","III","IV","V"))
data$adjacent_hepatic_tissue_inflammation_extent_type <-factor(data$adjacent_hepatic_tissue_inflammation_extent_type,
levels = c("1None","1Mild","3Severe" ),
labels = c("None","Mild","Severe"))
data$Virus_infection<-factor(data$Virus_infection,
levels = c("1No","2Hepatitis B or C","3Hepatitis B and C"),
labels = c("No","BorC","BandC"))
data$neoplasm_histologic_grade <-factor(data$neoplasm_histologic_grade,
levels = c("G1","G2" ,"G3" ,"G4" ),
labels = c("G1","G2" ,"G3" ,"G4" ))
data$pathologic_stage<-factor(data$pathologic_stage,
levels = c( "Stage I","Stage II", "Stage III" ,"Stage IV" ),
labels = c("I","II","III","IV"))
data$vascular_tumor_cell_type <-factor(data$vascular_tumor_cell_type,
levels = c("1None","2Micro","3Macro" ),
labels = c("None","Micro","Macro" ))做完了上面这一步就确定一下X和Y变量，赋值。
p<-names(data)[2:10]#确定自变量名称
q<-names(data)[11:16]#确定响应变量名称
t<-rep(1,9)#空变量模型中的自变量处填写1，我需要9个1第四步，批量建立公式，为了使用anova函数求P，所以这里建立空模型
#批量写公式 模型
l<-list()
for (y in q){
f<-function(x)as.formula(paste0(y,'~',x))
l[[y]]<-sapply(p, f)
}
ll<-unlist(l)
#批量写公式 空模型
l0<-list()
for (y in q){
f0<-function(x)as.formula(paste0(y,'~',x))
l0[[y]]<-sapply(t, f0)
}
ll0<-unlist(l0)第五步，批量建模
我这里分两种情况建模。因为我第一次做的时候使用了没有缺失值的数据，后来才发现有和没有比复杂很多，我这么写比较符合初学者的探索思路。我比较笨，没啥花花肠子炫技，我的思路其实相当质朴，就是把数据data拆分成了54个数据，一一对应建立模型，但是一样有效。看代码：
#针对没有缺失值的问题 批量建模
univ_models<-lapply(ll,function(x){polr(x,data=data,Hess = TRUE,method = "logistic")})
univ_models0<-lapply(ll0,function(x){polr(x,data=data,Hess = TRUE,method = "logistic")})
#针对有缺失值的问题 批量建模
x<-apply(data, MARGIN = 1, is.na)
x<-!x
XX<-rbind(x,x,x,x,x,x)#6为响应变量的个数
da<-list()
#54为响应变量的个数乘以自变量个数
for (i in 1:54) {
da[[i]]<-data[XX[i,],]
}
univ_models<-map2(ll,da,function(x,y){polr(x,data=y,Hess = TRUE,method = "logistic")})
univ_models0<-map2(ll0,da,function(x,y){polr(x,data=y,Hess = TRUE,method = "logistic")})第六步，对比模型和空模型，取得p值
#计算p
p<-map2(models0,models,anova)
#储存P
aa<-list()
for (i in 1: 54) {
aa[[i]]<-p[[i]][2,7]
}
pvalue<-unlist(aa)
#查看
print（pvalue）
第七步，批量输出结果
univ_results <- lapply(models,function(m){
OR<-round(exp(m$coefficients),7)
CI <- round(exp(confint(m)), 7)
CI<-data.frame(CI[1],CI[2])
colnames(CI) <- c("Lower", "Higher")
pv<- (pnorm(abs( coef(summary(m))[,"t value"]),lower.tail = FALSE)*2)[1]
results_list<-c(pv,CI, OR)
names(results_list)<-c("p.value", "Lower","Higher", "OR")
return(results_list)
})
r<- t(as.data.frame(univ_results, check.names = FALSE))
model_sum1dematrix <- data.frame(matrix(unlist(univ_results), nrow = 54,byrow=T),stringsAsFactors=FALSE)
names(model_sum1dematrix)<-c("p.value", "Lower","Higher", "OR")
rownames(model_sum1dematrix )<-names(ll)
#写出csv
write.csv(model_sum1dematrix,file = "multilogistc.csv")
写完这里我打算理个发，封面使用临高启民元老院flag，下次再见。
