在前面文章中介绍了无序多分类Logistic回归分析(Multinomial Logistic Regression Analysis)的假设检验理论，本篇文章将实例演示在R软件中实现无序多分类Logistic回归分析的操作步骤。
关键词：R语言; R软件; 无序多分类logistic回归; 无序logistic回归; 无序逻辑回归
一、案例介绍
欲探索性别(1=男性、2=女性)与年龄(1=<40岁、2=40-59岁、3=≥60岁)是否对某中医证型(1=A型、2=B型、3=C型)的分类有影响，从医院数据库中随机选择了200例样本进行分析。部分数据见图1。本文案例可从“附件下载”处下载。
二、问题分析
本案例中，因变量“中医证型”为无序多分类变量，欲探索“性别”与“年龄”是否对“中医证型”分类有影响，可将“中医证型”的某一类别设置为对照组，通过无序多分类logistic回归分析将另外两种不同类别证型的样本分别与对照组进行对比，得到“性别”、“年龄”与这两种类别证型的暴露-风险关系。无序多分类logistic回归需要满足3个条件：
条件1：因变量唯一，且为无序多分类变量，本案例符合。
条件2：存在一个或多个自变量，可为定性与定量变量，本案例符合。
条件3：一般要求例数较少类的样本量为自变量个数的10~15倍(EPV原则)且经验上每组的样本例数最好多于30例，参照水平组不应少于30或50例。
三、软件操作及结果解读
(一) 导入数据
mydata <- read.csv("无序多分类Logistic回归.csv") #导入CSV数据
View(mydata) #查看数据
在数据栏目中可以查看全部数据情况，数据集中共有3个变量和200个观察数据，3个变量分别代表性别(gender)、年龄(age)及中医证型(TCM)。
如果数据集较大也可使用如下命令查看数据框结构：
str(mydata) #查看数据框结构
(二) 适用条件判断
1. 条件3判断(因变量样本例数)
(1)软件操作
##给因变量添加标签##
mydata$TCM <- factor(mydata$TCM, levels = c(1, 2, 3), labels=c("A", "B", "C"))
##因变量计数##
table(mydata$TCM)
(2)结果解读
由图3的TCM(中医证型)”可见，因变量“A”、“B”、“C”三组例数分别为60、74和66，均>30例。根据“例数较少类的因变量例数为自变量个数的10~15倍(EPV原则)”，本案例可纳入4~6个自变量进行多因素无序多分类logistic回归分析。
2. 条件3判断(组别样本例数)
逐一计算分类变量各类别的因变量例数。
(1) 软件操作
##给自变量添加标签label ##
mydata$age <- factor(mydata$age,levels = c(1,2,3),labels=c("<40","40-59","≥60"))
mydata$gender <- factor(mydata$gender, levels = c(1, 2), labels=c("male", "female"))
##组别样本计数##
table(mydata$TCM, mydata$gender)
table(mydata$TCM, mydata$age)
(2) 结果解读
由图4可知，“年龄”水平为“<40”和“≥60”时因变量的例数<30，如果“年龄”在多因素分析过程中进入模型，应注意避免例数较少的水平被选为参照，若实在不能避免，在结果解释时应注意其局限性。
(三) 变量筛选
1. 软件操作
##层次回归分析## library(nnet) #载入包“nnet” fit1.1<-multinom(TCM~age, data=mydata) fit1.2<-multinom(TCM~gender, data=mydata) fit<-multinom(TCM ~ age + gender, data=mydata) ##依次与全变量模型比较## anova(fit1.1,fit,test = "Chisq") anova(fit1.2,fit,test = "Chisq")
2. 结果解读
图5显示了每个模型去除某1个自变量后与包含全部自变量的模型相比差异是否有统计学意义，也即表示去除的某个变量是否有统计学意义。如图5所示，模型“fit1.1”与模型“fit”相比差异尚无统计学意义，说明是否移除“gender”变量无统计学意义，即变量“gender”可以被移除模型。但考虑到本案例只有2个自变量，此次分析过程中仍然考虑保留在多因素分析模型中。
(四) 模型拟合
1. 软件操作
“中医证型”的参照水平先选择为“C”，即分析结果为“A”与“C”对比的情况、“B”与“C”对比的情况。调整自变量参照组，以年龄“≥60”为参照组(尽管“≥60”水平的例数<30，此处为了演示更多有统计学意义的结果，仍选择该水平为参照)。
##设置参照水平并进行模型拟合## mydata$age<-relevel(mydata$age,ref = "≥60") mydata$TCM<-relevel(mydata$TCM,ref = "C") fit<-multinom(TCM ~ age + gender, data=mydata) summary(fit) #查看模型拟合详细结果
##计算模型中变量系数的95%CI ##
confint(fit) #计算系数的95%CI
##计算OR## exp(coef(fit)) #计算OR exp(confint(fit)) #计算OR的95%CI
##计算P值## z <- summary(fit) $ coefficients/summary(fit) $ standard.errors p <- (1-pnorm(abs(z),0,1))*2 p
2. 结果解读
图6展示了模型的参数，列出了截距和自变量的“Coefficients (非标准化系数)”、“Std.Error (标准误)”、“Residual Deviance(残差)”和AIC值。图7则显示了模型拟合的所有变量系数的95%CI。图8计算了OR及其95%CI。图9计算了各系数与参照组比较对应的P值。
由图9可知，在“中医证型”A与C比较中，性别的P=0.040，但这并不能表示性别在模型中具有统计学意义，因为该变量的整体检验并无统计学意义。由图8可知，年龄<40岁的患者出现A型的风险是年龄≥60岁的7.133倍 (95%CI：2.063~24.664；P=0.002)，年龄40-59岁的患者出现A型的风险是年龄≥60岁的4.718倍(95%CI：1.916~11.619；P=0.001)。同理，在“中医证型”B与C的比较中，性别无统计学意义。年龄<40岁的患者出现B型的风险是年龄≥60岁的7.234倍 (95%CI：2.351~22.260；P<0.001)，年龄40-59岁的患者出现B型的风险是年龄≥60岁的2.902倍 (95%CI：1.281~6.573；P=0.011)。
(五) 补充分析
上述分析是以因变量C水平为参照得出的分析结果，可知B与C、A与C相比的情况，但是不知道A与B相比的情况。此处以水平A为参照进行补充分析。
1. 软件操作
分析过程、解读结果与上文类似。
##设置参照水平并进行模型拟合## mydata$TCM<-relevel(mydata$TCM,ref = "A") fit1<-multinom(TCM ~ age + gender, data=mydata) summary(fit1) #查看模型拟合详细结果。
##计算模型中变量系数的95%CI ##
confint(fit1) #计算系数的95%CI
##计算OR## exp(coef(fit1)) #计算OR exp(confint(fit1)) #计算OR的95%CI
##计算P值## z <- summary(fit1) $ coefficients/summary(fit) $ standard.errors p <- (1-pnorm(abs(z),0,1))*2 p
2. 结果解读
图10展示了模型的参数，列出了截距和自变量的“Coefficients (非标准化系数)”、“Std.Error (标准误)”、“Residual Deviance(残差)”和AIC值。图11则显示了模型拟合的所有变量系数的95%CI。图12计算了OR及其95%CI。图13计算了各系数与参照组比较对应的P值。
从图13可以看出，在“中医证型”B与A比较中，性别的P=0.0419，但同样这并不能表示性别在模型中具有统计学意义，因为该变量的整体检验并无统计学意义。年龄<40岁的患者和年龄40-59岁的患者出现B型的风险与年龄≥60岁的患者相比，差异也均无统计学(P=0.980和P=0.244)。C与A比较的结果等价于图6中A与C比较的结果。
四、结论
本研究采用无序多分类Logistic回归探讨性别和年龄是否对某中医证型的分类有影响。因变量例数分布满足样本量需求。
“Likelihood Ratio Tests(似然比)”检验结果表明，自变量“gender”在模型中无统计学意义 (χ²=5.721，P=0.057)，“年龄”在模型中有统计学意义 (χ²=21.875，P<0.001)。在“中医证型”A与C比较中，性别无统计学意义。年龄<40岁的患者出现A型的风险是年龄≥60岁的7.133倍 (95%CI：2.063~24.664；P=0.002)，年龄40-59岁的患者出现A型的风险是年龄≥60岁的4.718倍(95%CI：1.916~11.619；P=0.001)。在“中医证型”B与C的比较中，性别无统计学意义。年龄<40岁的患者出现B型的风险是年龄≥60岁的7.234倍 (95%CI：2.351~22.260；P<0.001)，年龄40-59岁的患者出现B型的风险是年龄≥60岁的2.902倍 (95%CI：1.281~6.573；P=0.011)。
在“中医证型”B与A比较中，性别和年龄均无统计学意义。
五、知识小贴士
- 无序多分类logistic回归的因变量为多分类变量，内部分析过程是选择一个参照组将因变量拆分为多个二分类变量，拟合多个二分类logistic回归。
- 由于无序多分类logistic回归的本质是多个二分类logistic回归，因此其结果解读可参考二分类logistic回归，只是每个变量需要在多个二分类logistic回归中分别解读。
- 在“中医证型”B与C的比较中，可发现自变量“性别”有统计学意义，但由于整体检验该变量并无统计学意义，因此尚不能认为“性别”对中医证型的分类有影响。