R语言进行无序多分类Logistic回归

在临床研究中，接触最多的是二分类数据，如淋巴癌是否转移，是否死亡，这些因变量最后都可以转换成二分类0与1的问题。然后建立二元logistic回归方程，可以得到影响因素的OR值。但有时我们也会接触到多分类结局数据，今天咱们来演示一下怎么使用R语言进行多分类结局逻辑回归分析。

咱们先导入数据和R包

library(here)
library(tidyverse)
library(haven)
library(gtsummary)
library(VGAM)
library(kableExtra)
library(HSAUR)
library(reshape2)
setwd("E:/公众号文章2024年/R无序多分类Logistic回归")
load("hsbdemo.rda")
ms<-hsbdemo

这个数据是高中生毕业后的一个就业计划数据，Prog是结局变量，是个三分类变量，预测变量是ses社会经济地位，其他的是一些协变量。

变量比较多，咱们选出需要的变量，并且把字符变量转成因子

ms <- ms %>% 
  select(ses, prog, female , write ) %>% 
  mutate(across(where(is.labelled), as_factor)) 

咱们以prog为分类变量，绘制个基线表，了解相关数据分布

ms %>%
  tbl_summary(by = prog,
              statistic = list(all_continuous() ~ "{mean} ({sd})", 
                               all_categorical() ~ "{n} ({p}%)"),
              type = list(where(is.logical) ~ "categorical")) %>% 
  modify_caption("**Table 1. Survey Participant Characteristic**")  %>%
  modify_header(label ~ "**Variable**") %>%
  modify_spanning_header(c("stat_1", "stat_2", "stat_3") ~ "**Baseline tables for three types of projects**") %>%
  modify_footnote(all_stat_cols() ~ "Mean (SD) or Frequency (%)") %>%
  bold_labels() %>%
  as_gt()

我们以academic为参考，建立新的变量

ms <- ms %>% 
  mutate(prog2 = fct_relevel(prog, 
                               c("academic", 'general', 'vocation')))
levels(ms$prog2)

更改了参考类别后，咱们使用prog2为结局变量建立无序多分类逻辑回归模型，很多包可以建立这个模型，我这里用vglm包来建，我觉得比较简单点

fit <- vglm(prog2~ ses + write, 
                 multinomial, data = ms)
summary(fit)

这样结果就出来啦，seslow:1和seslow:2的这两个系数是针对academic这个结局的
还可以做交互效应的，把交互效应打上去就可以啦，我这里就不弄了。

查看系数和可惜区间

b_mlog <- coef(fit )
ci_mlog <- confint(fit)
b_ci_mlog <- data.frame(b_mlog,ci_mlog) %>%
  rename("log odds" = b_mlog, "Lower CI" = X2.5.., "Upper CI" = X97.5..)
b_ci_mlog %>% 
  kbl(digits = 2, booktabs = T, caption = "Log odds from multinomial logistic regression") %>%
  kable_styling(position = "center")

查看各个结局概率的数据

predict.vgam(fit, type = 'response') %>% 
  head(20)

构造一个新数据

ms2 <- data.frame(ses = rep(c("low", "middle", "high"), each = 41), write = rep(c(30:70), 
                                                                                   3))
ms2<-predict(fit, newdata = ms2,type = 'response') %>% 
  cbind(ms2)

转化数据

ms3 <- melt(ms2, id.vars = c("ses", "write"), value.name = "probability")

最后绘图

ggplot(ms3, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~ 
                                                                                        ., scales = "free")

表明随着write增加，不同的经济基础转向各个计划的概率不同。

参考文献

1. https://xianxiongma.github.io/Clinical-model/chapter2/chapter2.html
2. https://bookdown.org/drki_musa/dataanalysis/multinomial-logistic-regression.html