“医学和生信笔记，专注R语言在临床医学中的使用、R语言数据分析和可视化。主要分享R语言做医学统计学、临床研究设计、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。
R语言二项逻辑回归： R语言logistic回归的细节解读
R语言多项逻辑回归：R语言多项逻辑回归-因变量是无序多分类
有序逻辑回归
ordinal logistic regression适用于因变量为等级资料。使用孙振球版医学统计学例16-4的数据。
随机选取84例患者做临床试验，探讨性别和治疗方法对该病的影响。变量赋值为：性别（X1，男=0，女=1），治疗方法（X2，传统疗法=0，新型疗法=1），疗效（Y，无效=1，有效=2，痊愈=3）。
代码语言：javascript复制df <- read.csv("例16-04.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
## ... ... ... ...
## 81 1 1 3
## 82 1 1 3
## 83 1 1 3
## 84 1 1 3
变为因子型：
代码语言：javascript复制# 因变量变为有序因子
df$Y <- factor(df$Y, levels = c(1,2,3),
labels = c("无效","有效","痊愈"),
ordered = T)
# 自变量变为无序因子
df$X1 <- factor(df$X1,levels = c(0,1),labels = c("男","女"))
df$X2 <- factor(df$X2,levels = c(0,1),labels = c("传统疗法","新型疗法"))
str(df)
## 'data.frame': 84 obs. of 3 variables:
## $ X1: Factor w/ 2 levels "男","女": 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 : Ord.factor w/ 3 levels "无效"<"有效"<..: 1 1 1 1 1 1 1 1 1 1 ...
使用MASS::polr拟合有序逻辑回归：
library(MASS)
fit <- polr(Y ~ X1 X2, data = df,Hess = TRUE,method = "logistic")
summary(fit)
## Call:
## polr(formula = Y ~ X1 X2, data = df, Hess = TRUE, method = "logistic")
##
## Coefficients:
## Value Std. Error t value
## X1女 1.319 0.5381 2.451
## X2新型疗法 1.797 0.4718 3.809
##
## Intercepts:
## Value Std. Error t value
## 无效|有效 1.8128 0.5654 3.2061
## 有效|痊愈 2.6672 0.6065 4.3979
##
## Residual Deviance: 150.0294
## AIC: 158.0294
结果也是没有给出P值，手动计算P值：
代码语言：javascript复制p <- pnorm(abs(coef(summary(fit))[, "t value"]),lower.tail = F)*2
p
## X1女 X2新型疗法 无效|有效 有效|痊愈
## 1.425572e-02 1.392807e-04 1.345300e-03 1.092866e-05
这次broom::tidy(fit)并没有直接给出P值：
broom::tidy(fit)
## # A tibble: 4 × 5
## term estimate std.error statistic coef.type
## <chr> <dbl> <dbl> <dbl> <chr>
## 1 X1女 1.32 0.538 2.45 coefficient
## 2 X2新型疗法 1.80 0.472 3.81 coefficient
## 3 无效|有效 1.81 0.565 3.21 scale
## 4 有效|痊愈 2.67 0.606 4.40 scale
OR值：
代码语言：javascript复制OR <- exp(coef(fit))
OR
## X1女 X2新型疗法
## 3.738765 6.033338
平行线检验(Brant-Wald test)：
代码语言：javascript复制brant::brant(fit)
## --------------------------------------------
## Test for X2 df probability
## --------------------------------------------
## Omnibus 1.83 2 0.4
## X1女 1.59 1 0.21
## X2新型疗法 0.01 1 0.94
## --------------------------------------------
##
## H0: Parallel Regression Assumption holds
P值>0.05，平行线检验通过，可以使用有序逻辑回归，通不过可以用多项逻辑回归。
模型整体的显著性检验：
代码语言：javascript复制# 先构建一个只有截距的模型
fit0 <- polr(Y ~ 1, data = df,Hess = TRUE,method = "logistic")
# 两个模型比较
anova(fit0, fit)
## Likelihood ratio tests of ordinal regression models
##
## Response: Y
## Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi)
## 1 1 82 169.9159
## 2 X1 X2 80 150.0294 1 vs 2 2 19.8865 4.80508e-05
P值＜0.01，模型是有意义的。
获取模型预测的类别：
代码语言：javascript复制pred <- predict(fit, df, type = "class")
head(pred)
## [1] 无效 无效 无效 无效 无效 无效
## Levels: 无效 有效 痊愈
获取模型预测的概率：
代码语言：javascript复制prob <- predict(fit, df, type = "probs") # 或者使用 fitted(fit)
head(prob)
## 无效 有效 痊愈
## 1 0.8597003 0.07536263 0.06493706
## 2 0.8597003 0.07536263 0.06493706
## 3 0.8597003 0.07536263 0.06493706
## 4 0.8597003 0.07536263 0.06493706
## 5 0.8597003 0.07536263 0.06493706
## 6 0.8597003 0.07536263 0.06493706
模型拟合优度的检验，这里使用卡方检验：
代码语言：javascript复制chisq.test(df$Y, pred)
##
## Pearson's Chi-squared test
##
## data: df$Y and pred
## X-squared = 14.246, df = 2, p-value = 0.0008065
计算伪R2：
代码语言：javascript复制DescTools::PseudoR2(fit, which = "all")
## McFadden CoxSnell Nagelkerke AldrichNelson VeallZimmermann
## 0.1170373 0.2108068 0.2429443 NA NA
## Efron McKelveyZavoina Tjur AIC BIC
## NA NA NA 158.0294131 167.7526803
## logLik logLik0 G2
## -75.0147065 -84.9579583 19.8865036
不仅给出了伪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/
