在前面文章中介绍了无序多分类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) #查看数据
![](https://site.cdn.mengte.online/official/2022/02/20220221102619870.png)
在数据栏目中可以查看全部数据情况,数据集中共有3个变量和200个观察数据,3个变量分别代表性别(gender)、年龄(age)及中医证型(TCM)。
如果数据集较大也可使用如下命令查看数据框结构:
str(mydata) #查看数据框结构
![](https://site.cdn.mengte.online/official/2022/02/20220221102604548.png)
(二) 适用条件判断
1. 条件3判断(因变量样本例数)
(1)软件操作
##给因变量添加标签## mydata$TCM <- factor(mydata$TCM, levels = c(1, 2, 3), labels=c("A", "B", "C")) ##因变量计数## table(mydata$TCM)
![](https://site.cdn.mengte.online/official/2022/02/20220221102520716.png)
(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)
![](https://site.cdn.mengte.online/official/2022/02/20220221102435771.png)
(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")
![](https://site.cdn.mengte.online/official/2022/02/20220221102344893.png)
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) #查看模型拟合详细结果
![](https://site.cdn.mengte.online/official/2022/02/20220221102256342.png)
##计算模型中变量系数的95%CI ##
confint(fit) #计算系数的95%CI
![](https://site.cdn.mengte.online/official/2022/02/20220221102234149.png)
##计算OR## exp(coef(fit)) #计算OR exp(confint(fit)) #计算OR的95%CI
![](https://site.cdn.mengte.online/official/2022/02/20220221102208282.png)
##计算P值## z <- summary(fit) $ coefficients/summary(fit) $ standard.errors p <- (1-pnorm(abs(z),0,1))*2 p
![](https://site.cdn.mengte.online/official/2022/02/20220221102135513.png)
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) #查看模型拟合详细结果。
![](https://site.cdn.mengte.online/official/2022/02/20220221102037573.png)
##计算模型中变量系数的95%CI ##
confint(fit1) #计算系数的95%CI
![](https://site.cdn.mengte.online/official/2022/02/20220221102018974.png)
##计算OR## exp(coef(fit1)) #计算OR exp(confint(fit1)) #计算OR的95%CI
![](https://site.cdn.mengte.online/official/2022/02/20220221101946374.png)
##计算P值## z <- summary(fit1) $ coefficients/summary(fit) $ standard.errors p <- (1-pnorm(abs(z),0,1))*2 p
![](https://site.cdn.mengte.online/official/2022/02/20220221101901718.png)
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的比较中,可发现自变量“性别”有统计学意义,但由于整体检验该变量并无统计学意义,因此尚不能认为“性别”对中医证型的分类有影响。