R语言多项逻辑回归-因变量是无序多分类

“医学和生信笔记，专注R语言在临床医学中的使用、R语言数据分析和可视化。主要分享R语言做医学统计学、临床研究设计、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。

R语言二项逻辑回归：R语言logistic回归的细节解读

多项逻辑回归

因变量是无序多分类资料（＞2）时，可使用多分类逻辑回归（multinomial logistic regression）。

使用孙振球医学统计学第四版例16-5的数据，课本电子版及数据已上传到QQ群，自行下载即可。

某研究人员欲了解不同社区和性别之间居民获取健康知识的途径是否相同，对2个社区的314名成人进行了调查，其中X1是社区，社区1用0表示，社区2用1表示；X2是性别，0是男，1是女，Y是获取健康知识途径，1是传统大众传媒，2是网络，3是社区宣传。

代码语言：javascript复制

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()重新设置参考。

代码语言：javascript复制

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回归：

代码语言：javascript复制

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值需要手动计算:

代码语言：javascript复制

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

但其实可以调包解决：

代码语言：javascript复制

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值和可信区间：

代码语言：javascript复制

# 计算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

模型整体的假设检验：

代码语言：javascript复制

# 先构建一个只有截距的模型
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，模型具有统计学意义。

获取模型预测的类别：

代码语言：javascript复制

pred <- predict(fit, df, type = "class")
head(pred)
## [1] 网络 网络 网络 网络 网络 网络
## Levels: 传统大众传媒 网络 社区宣传

获取模型预测的概率：

代码语言：javascript复制

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

模型拟合优度的检验，这里使用卡方检验：

代码语言：javascript复制

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：

代码语言：javascript复制

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/

https 网络安全 r语言

0 人点赞