多分类logistic回归

在临床研究中，接触最多的是二分类数据，如淋巴癌是否转移，是否死亡，这些因变量最后都可以转换成二分类0与1的问题。然后建立二元logistic回归方程，可以得到影响因素的OR值。

那么如果遇到多分类变量，如何进行logistic回归呢？譬如临床疗效分为好，中，差，三类，或者根据指标进行分类，分为高，中，低三类，我用１、２、３代表作为因变量，进行logistic回归分析。

image.png

接下来，该文，主要介绍，如果因变量为三分类变量，如何进行回归分析及机器学习算法对三分类资料的处理。关于原理理论部分可参见；这里主要讲如何在R实现三分类回归，计算系数及p值与OR值

1.数据案例

这里主要用到DALEX包里面包含的HR数据，里面记录了职工在工作岗位的状态与年龄，性别，工作时长，评价及薪水有关。根据7847条记录来评估，如果一个职工属于男性，68岁，薪水及评价处于3等级，那么该职工可能会处于什么状态。

library(DALEX)

library(iBreakDown)

library(nnet)

library(questionr)

try(data(package="DALEX"))

data(HR)

HR= HR %>% as.tbl() %>%

mutate(evaluation=factor(evaluation),

salary=factor(salary))

HR

## GLM

fit = multinom(status ~ . , data = HR, probabilities = TRUE, model = TRUE)

summary(fit)

coef(fit)

> summary(fit)

Call:

multinom(formula = status ~ ., data = HR, model = TRUE, probabilities = TRUE)

Coefficients:

(Intercept) gendermale age hours evaluation3 evaluation4 evaluation5

ok -5.47276 0.03437426 0.002594237 0.08305463 -0.07275332 -0.06763166 -0.156932

promoted -13.10377 0.10391193 0.004277562 0.19697483 -0.11679839 3.49127986 3.290217

salary1 salary2 salary3 salary4 salary5

ok 1.543631 2.469598 2.413207 1.758516 -0.09513189

promoted 1.650777 2.498608 2.435680 1.790657 -0.01215312

Std. Errors:

(Intercept) gendermale age hours evaluation3 evaluation4

ok 0.2407953 0.06427342 0.002784445 0.003669856 0.07434565 0.1061762

promoted 0.3475852 0.08023646 0.003458588 0.004692886 0.11569077 0.1303817

evaluation5 salary1 salary2 salary3 salary4 salary5

ok 0.1077599 0.1199301 0.1227032 0.1219577 0.1212186 0.1371131

promoted 0.1302725 0.1456579 0.1490790 0.1482045 0.1486752 0.1624046

Residual Deviance: 10744.64

AIC: 10792.64

> coef(fit)

(Intercept) gendermale age hours evaluation3 evaluation4 evaluation5

ok -5.47276 0.03437426 0.002594237 0.08305463 -0.07275332 -0.06763166 -0.156932

promoted -13.10377 0.10391193 0.004277562 0.19697483 -0.11679839 3.49127986 3.290217

salary1 salary2 salary3 salary4 salary5

ok 1.543631 2.469598 2.413207 1.758516 -0.09513189

promoted 1.650777 2.498608 2.435680 1.790657 -0.01215312

我们构建了三元回归模型，以status中fired为参照，计算ok与promoted中各个因素的系数。

有了这些系数，我们就可以写出回归方程了，然后再计算各个因素对应的p值

如，这里的例子介绍了其他因素的系数，然后计算对因变量的方程here

image.png

P值

通过Anova函数，可以输出fit中影响因素的p值，其中hours，evaluation及salary有统计学意义。说明他们对员工在职影响很大。然后进一步计算or值。

需要借助questionr包中的odds.ratio函数。

> Anova(fit)

Analysis of Deviance Table (Type II tests)

Response: status

LR Chisq Df Pr(>Chisq)

gender 1.7 2 0.4299

age 1.7 2 0.4329

hours 3464.1 2 <2e-16 ***

evaluation 2390.2 6 <2e-16 ***

salary 1132.4 10 <2e-16 ***

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> round(odds.ratio(fit),2)

OR 2.5 % 97.5 % p

ok/(Intercept) 0.00 0.00 0.01 <2e-16 ***

ok/gendermale 1.03 0.91 1.17 0.59

ok/age 1.00 1.00 1.01 0.35

ok/hours 1.09 1.08 1.09 <2e-16 ***

ok/evaluation3 0.93 0.80 1.08 0.33

ok/evaluation4 0.93 0.76 1.15 0.52

ok/evaluation5 0.85 0.69 1.06 0.15

ok/salary1 4.68 3.70 5.92 <2e-16 ***

ok/salary2 11.82 9.29 15.03 <2e-16 ***

ok/salary3 11.17 8.79 14.19 <2e-16 ***

ok/salary4 5.80 4.58 7.36 <2e-16 ***

ok/salary5 0.91 0.69 1.19 0.49

promoted/(Intercept) 0.00 0.00 0.00 <2e-16 ***

promoted/gendermale 1.11 0.95 1.30 0.20

promoted/age 1.00 1.00 1.01 0.22

promoted/hours 1.22 1.21 1.23 <2e-16 ***

promoted/evaluation3 0.89 0.71 1.12 0.31

promoted/evaluation4 32.83 25.43 42.39 <2e-16 ***

promoted/evaluation5 26.85 20.80 34.66 <2e-16 ***

promoted/salary1 5.21 3.92 6.93 <2e-16 ***

promoted/salary2 12.17 9.08 16.29 <2e-16 ***

promoted/salary3 11.42 8.54 15.27 <2e-16 ***

promoted/salary4 5.99 4.48 8.02 <2e-16 ***

promoted/salary5 0.99 0.72 1.36 0.94

---

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

参考