其他
R语言多项逻辑回归-因变量是无序多分类
“医学和生信笔记，专注R语言在临床医学中的使用、R语言数据分析和可视化。主要分享R语言做医学统计学、临床研究设计、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
💡专注R语言在🩺生物医学中的使用
R语言二项逻辑回归：R语言logistic回归的细节解读
多项逻辑回归
因变量是无序多分类资料（＞2）时，可使用多分类逻辑回归（multinomial logistic regression）。
使用孙振球医学统计学第四版例16-5的数据，课本电子版及数据已上传到QQ群，自行下载即可。
某研究人员欲了解不同社区和性别之间居民获取健康知识的途径是否相同，对2个社区的314名成人进行了调查，其中X1是社区，社区1用0表示，社区2用1表示；X2是性别，0是男，1是女，Y是获取健康知识途径，1是传统大众传媒，2是网络，3是社区宣传。
df <- read.csv("例16-05.csv",header = T)
psych::headtail(df)
## X1 X2 Y
## 1 0 0 1
## 2 0 0 1
## 3 0 0 1
## 4 0 0 1
## ... ... ... ...
## 311 1 1 3
## 312 1 1 3
## 313 1 1 3
## 314 1 1 3
首先变为因子型，无需多分类的logistic回归需要对因变量设置参考，我们这里直接用factor()函数变为因子，这样在进行无序多分类的logistic时默认是以第一个为参考。也可以使用relevel()重新设置参考。
df$X1 <- factor(df$X1,levels = c(0,1),labels = c("社区1","社区2"))
df$X2 <- factor(df$X2,levels = c(0,1),labels = c("男","女"))
# 因变量设置参考，这里选择第1个（传统大众传媒）为参考
df$Y <- factor(df$Y,levels = c(1,2,3),labels = c("传统大众传媒","网络","社区宣传"))
str(df)
## 'data.frame': 314 obs. of 3 variables:
## $ X1: Factor w/ 2 levels "社区1","社区2": 1 1 1 1 1 1 1 1 1 1 ...
## $ X2: Factor w/ 2 levels "男","女": 1 1 1 1 1 1 1 1 1 1 ...
## $ Y : Factor w/ 3 levels "传统大众传媒",..: 1 1 1 1 1 1 1 1 1 1 ...
使用nnet::multinom进行无序多分类的logistic回归：
library(nnet)
fit <- multinom(Y ~ X1 + X2, data = df, model = T)
## # weights: 12 (6 variable)
## initial value 344.964259
## iter 10 value 316.575399
## iter 10 value 316.575399
## iter 10 value 316.575399
## final value 316.575399
## converged
summary(fit)
## Call:
## multinom(formula = Y ~ X1 + X2, data = df, model = T)
##
## Coefficients:
## (Intercept) X1社区2 X2女
## 网络 0.5484998 -1.3743147 0.4321069
## 社区宣传 0.3940422 -0.9933526 1.2266459
##
## Std. Errors:
## (Intercept) X1社区2 X2女
## 网络 0.2583299 0.3201514 0.3265384
## 社区宣传 0.2574175 0.2952083 0.2991714
##
## Residual Deviance: 633.1508
## AIC: 645.1508
可以看到结果比二项逻辑回归的结果简洁多了，少了很多信息，只给出了Coefficients/Std. Errors/Residual Deviance/AIC。
不过也是两个模型的结果，分别是 社区宣传 和 传统大众传媒 比，网络 和 传统大众传媒 比。
自变量的Z值（wald Z, Z-score）和P值需要手动计算:
z_stats <- summary(fit)$coefficients/summary(fit)$standard.errors
p_values <- (1 - pnorm(abs(z_stats)))*2
p_values
## (Intercept) X1社区2 X2女
## 网络 0.03373263 1.765117e-05 1.857371e-01
## 社区宣传 0.12583082 7.656564e-04 4.128929e-05
但其实可以调包解决：
res <- broom::tidy(fit)
res
## # A tibble: 6 × 6
## y.level term estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 网络 (Intercept) 0.548 0.258 2.12 0.0337
## 2 网络 X1社区2 -1.37 0.320 -4.29 0.0000177
## 3 网络 X2女 0.432 0.327 1.32 0.186
## 4 社区宣传 (Intercept) 0.394 0.257 1.53 0.126
## 5 社区宣传 X1社区2 -0.993 0.295 -3.36 0.000766
## 6 社区宣传 X2女 1.23 0.299 4.10 0.0000413
OR值及其95%的可信区间也没有给出来，需要手动计算OR值和可信区间：
# 计算OR值
OR <- exp(coef(fit))
OR
## (Intercept) X1社区2 X2女
## 网络 1.730655 0.2530129 1.540500
## 社区宣传 1.482963 0.3703330 3.409774
# 计算OR值的95%的可信区间
OR.confi <- exp(confint(fit))
OR.confi
## , , 网络
##
## 2.5 % 97.5 %
## (Intercept) 1.0430848 2.8714501
## X1社区2 0.1350919 0.4738666
## X2女 0.8122910 2.9215386
##
## , , 社区宣传
##
## 2.5 % 97.5 %
## (Intercept) 0.8953982 2.4560912
## X1社区2 0.2076398 0.6605021
## X2女 1.8970133 6.1288742
模型整体的假设检验：
# 先构建一个只有截距的模型
fit0 <- multinom(Y ~ 1, data = df, model = T)
## # weights: 6 (2 variable)
## initial value 344.964259
## final value 338.603448
## converged
# 两个模型比较,Likelihood ratio tests
anova(fit0, fit)
## Likelihood ratio tests of Multinomial Models
##
## Response: Y
## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 626 677.2069
## 2 X1 + X2 622 633.1508 1 vs 2 4 44.0561 6.245931e-09
P<0.001，模型具有统计学意义。
获取模型预测的类别：
pred <- predict(fit, df, type = "class")
head(pred)
## [1] 网络 网络 网络 网络 网络 网络
## Levels: 传统大众传媒 网络 社区宣传
获取模型预测的概率：
prob <- predict(fit, df, type = "probs") # 或者使用 fitted(fit)
head(prob)
## 传统大众传媒 网络 社区宣传
## 1 0.2373257 0.4107289 0.3519453
## 2 0.2373257 0.4107289 0.3519453
## 3 0.2373257 0.4107289 0.3519453
## 4 0.2373257 0.4107289 0.3519453
## 5 0.2373257 0.4107289 0.3519453
## 6 0.2373257 0.4107289 0.3519453
模型拟合优度的检验，这里使用卡方检验：
chisq.test(df$Y, pred)
##
## Pearson's Chi-squared test
##
## data: df$Y and pred
## X-squared = 39.521, df = 4, p-value = 5.436e-08
计算伪R^2：
DescTools::PseudoR2(fit, which = "all")
## McFadden McFaddenAdj CoxSnell Nagelkerke AldrichNelson
## 0.06505559 0.04733575 0.13090778 0.14803636 NA
## VeallZimmermann Efron McKelveyZavoina Tjur AIC
## NA NA NA NA 645.15079819
## BIC logLik logLik0 G2
## 667.64715610 -316.57539909 -338.60344772 44.05609725
不仅给出了伪R^2，还给出了超多的值，每一项的意义可以参考下面这张图：
结果解读可以参考二项逻辑回归。
参考资料
https://blog.csdn.net/weixin_41744624/article/details/105506951 https://zhuanlan.zhihu.com/p/113403422 https://duanku.pai-hang-bang.cn/kuzi_1046977453210716059 https://bookdown.org/chua/ber642_advanced_regression/ https://peopleanalytics-regression-book.org/
咨询交流等，欢迎加入🐧QQ交流群：613637742
往期推荐