无序多分类Logistic回归分析(Multinomial Logistic Regression Analysis)——R软件实现

发布于 2022年2月21日 星期一 18:26:57 浏览:6396
原创不易,转载请注明来源,感谢!
附件下载:
无序多分类logistic回归.zip 请勿重复点击,如无响应请耐心等待或稍后再试。

在前面文章中介绍了无序多分类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)  #查看数据
图1

在数据栏目中可以查看全部数据情况,数据集中共有3个变量和200个观察数据,3个变量分别代表性别(gender)、年龄(age)及中医证型(TCM)。

如果数据集较大也可使用如下命令查看数据框结构:

str(mydata)  #查看数据框结构

图2

(二) 适用条件判断

1. 条件3判断(因变量样本例数)

(1)软件操作
##给因变量添加标签##
mydata$TCM <- factor(mydata$TCM, levels = c(1, 2, 3), labels=c("A", "B", "C"))
##因变量计数##
table(mydata$TCM)
图3
(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)
图4
(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")
图5

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) #查看模型拟合详细结果
图6

##计算模型中变量系数的95%CI ##

confint(fit) #计算系数的95%CI

图7
##计算OR##
exp(coef(fit)) #计算OR 
exp(confint(fit)) #计算OR的95%CI
图8
##计算P值##
z <- summary(fit) $ coefficients/summary(fit) $ standard.errors
p <- (1-pnorm(abs(z),0,1))*2
p 
图9

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) #查看模型拟合详细结果。
图10

##计算模型中变量系数的95%CI ##

confint(fit1) #计算系数的95%CI

图11
##计算OR##
exp(coef(fit1)) #计算OR 
exp(confint(fit1)) #计算OR的95%CI
图12
##计算P值##
z <- summary(fit1) $ coefficients/summary(fit) $ standard.errors
p <- (1-pnorm(abs(z),0,1))*2
p 
图13

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的比较中,可发现自变量“性别”有统计学意义,但由于整体检验该变量并无统计学意义,因此尚不能认为“性别”对中医证型的分类有影响。
End
文章目录 沉浸式阅读