Multinomial Logit Model (MNL) 模型R语言nnet包multinom函数实现实例

最近做项目涉及到要使用multinomial logit model (MNL) 模型。看了一堆文献讲mnl，但是没有给什么具体能上手的实例，就算有也是一笔带过，打算找一些使用R 语言来实现mnl模型的例子，在模仿和实践中慢慢理解。

Multinomial Logit Model又有很多其它说法，诸如Multinomial Logistic Regression等等。

本文的实例来自两篇文章。

［1］R Data Analysis Examples: Multinomial Logistic Regression： http://www.ats.ucla.edu/stat/r/dae/mlogit.htm

［2］How to: Multinomial regression models in R ：https://www.r-bloggers.com/how-to-multinomial-regression-models-in-r/

第一篇 R Data Analysis Examples: Multinomial Logistic Regression

第一篇是UCLA的idre机构网站中，关于R语言实现 Multinomial Logistic Regression 的教程

Multinomial logistic regression被用于输出结果为 nominal variables 的建模。

本文使用了一下的包，请确保你能载入这些包，如果你没有安装，可以使用语句：install.packages("packagename"), 或者如果你使用的包的版本太低，可以使用语句： update.packages() .

require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)

Version info: Code for this page was tested in R version 3.1.1 (2014-07-10)
On: 2015-12-17
With: reshape2 1.4.1; ggplot2 1.0.1; nnet 7.3-10; foreign 0.8-65; knitr 1.10.5

Multinomial Logistic Regression的例子

例1: 人们的职业选择结果可能会被父母的职业和他们自己的教育水平所影响。我们可以研究某个人的职业选择和他的教育水平、父母的职业之间的关系。而人们所选择的职业多种多样，不是只有一种或者两种。

例2：一个生物学家啃呢个会对短吻鳄所选择的食物感兴趣。成年的短吻鳄可能与幼年的短吻鳄的食物偏好不同。所以我们这里的有各种各样的食物作为选择结果，该结果是被短吻鳄的形体大小和环境变量所影响。

例3：进入高等教育的学生对于选择什么样的学习项目类型有三种选择：普通项目，针对工作的项目以及学术型的项目。他们的选择可能被他们的写作成绩和社会经济地位影响。

数据描述

我们的数据分析例子使用的是第三个例子，使用 hsbdemo 数据。首先是先读入数据。

ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")

该数据包括200个学生的选择的项目类型(prog，三种类型 categorical variable), 他们的社会地位（ses 三种地位 categorical variable），写作分数（write， a continuous variable）。读完数据后，我们可以使用一些语句来对我们的数据建立一个初步的概念和感觉。

with(ml, table(ses, prog))

##         prog
## ses      general academic vocation
##   low         16       19       12
##   middle      20       44       31
##   high         9       42        7

从结果中可以看出，社会经济地位中等的学生最多，低等的最少。社会经济地位高的学生中，绝大多数都是选择学术型的program 。

with(ml, do.call(rbind, tapply(write, prog, function(x) c(M = mean(x), SD = sd(x)))))

##             M   SD
## general  51.3 9.40
## academic 56.3 7.94
## vocation 46.8 9.32

从结果中可以看出，选择学术型项目的学生的写作成绩平均分最高，且波动最小。选择职业型项目的学生写作成绩平均分最低。

Multinomial Logistic Regression 方法

以下的方法中我们使用的nnet 包中的multinom 函数。其实在R语言的其它包中也有其它函数可以实现MNL方法(诸如 mlogit)。我们选择multinom 函数的原因是，它不需要将数据reshape(而mlogit就需要)。

在运行我们的模型之前，选择一个合适的参照组（a reference group）非常重要。我们可以使用relevel 来调换outcome variable 的等级顺序。值得注意的是， multinom 包不包含对于回归系数的p-value的计算。

ml$prog2 <- relevel(ml$prog, ref = "academic")
test <- multinom(prog2 ~ ses + write, data = ml)

## # weights:  15 (8 variable)
## initial  value 219.722458
## iter  10 value 179.982880
## final  value 179.981726
## converged

summary(test)

## Call:
## multinom(formula = prog2 ~ ses + write, data = ml)
##
## Coefficients:
##          (Intercept) sesmiddle seshigh   write
## general         2.85    -0.533  -1.163 -0.0579
## vocation        5.22     0.291  -0.983 -0.1136
##
## Std. Errors:
##          (Intercept) sesmiddle seshigh  write
## general         1.17     0.444   0.514 0.0214
## vocation        1.16     0.476   0.596 0.0222
##
## Residual Deviance: 360
## AIC: 376

z <- summary(test)$coefficients/summary(test)$standard.errors
z

##          (Intercept) sesmiddle seshigh write
## general         2.45    -1.202   -2.26 -2.71
## vocation        4.48     0.612   -1.65 -5.11

#2-tailed z test
p <- (1 - pnorm(abs(z), 0, 1))*2
p

##          (Intercept) sesmiddle seshigh    write
## general     1.45e-02     0.229  0.0237 6.82e-03
## vocation    7.30e-06     0.541  0.0989 3.18e-07

我们首先可以看到，即使我们已经将multinom运行赋值给了test，模型运行后仍然得到了一些结果。运行结果包括 some iteration history 和最终的负log-likelihood 179.981726. 这一结果乘以2后即为模型summary结果中的 Residual Deviance:360，它可以用来和 nested models 来比较，但是本例中我们不做比较。

模型summary 后的结果一块是系数，一块是 standard errors。每一个块中，都有一行与模型等式对应。在系数块中，第一行是general program 与academic program 的比较，第二行代表 vocation program 与 academic program 的比较（academic program 是我们的参考组reference group。如果我们把第一行的系数记为b1_, 第二行的系数记为b2_, 我们可以写出我们的模型等式。

ln(P(prog=general)P(prog=academic))=b10+b11(ses=2)+b12(ses=3)+b13write

ln(P(prog=vocation)P(prog=academic))=b20+b21(ses=2)+b22(ses=3)+b23write

一种选择结果类别的概率与基准选择的概率的比被称为相对危险度（relative risk）, 有时候只是为了描述回归参数，我们也叫做odds. 相对危险度是对线性方程等号右边部分的取幂，取幂后的回归系数就是自变量每单位变化的相对危险度。下面我们开对模型的系数取幂来观察是如何变化的。

## extract the coefficients from the model and exponentiate
exp(coef(test))

##          (Intercept) sesmiddle seshigh write
## general         17.3     0.587   0.313 0.944
## vocation       184.6     1.338   0.374 0.893

“write”变量增加一个单位，普通项目vs学术项目的相对危险风险比（the relative risk ratio）是0.9437
“ses”社会地位变量从1变为3时，普通项目VS学术项目的相对风险比是0.3126

你也可以使用预测的概率来帮助你理解这个模型。你可以使用 fitted 函数来得到模型的拟合值（估计值）。（在这里和原数据比对了一下，感觉不是很精确呀。）

head(pp <- fitted(test))

##   academic general vocation
## 1    0.148   0.338    0.513
## 2    0.120   0.181    0.699
## 3    0.419   0.237    0.345
## 4    0.173   0.351    0.476
## 5    0.100   0.169    0.731
## 6    0.353   0.238    0.409

然后，如果你想检测我们变量中某个变量的改变对预测结果的影响，可以创建一个小的数据组，改变其中一个变量，其它变量保持不变。首先，我们让"write" 变量保持不变，检测社会地位的改变对预测值的影响。

dses <- data.frame(ses = c("low", "middle", "high"),write = mean(ml$write))
predict(test, newdata = dses, "probs")

##   academic general vocation
## 1    0.440   0.358    0.202
## 2    0.478   0.228    0.294
## 3    0.701   0.178    0.121

另一种理解模型的方式是在三种不同的社会地位下，连续改变"write"值，并对该社会地位中的预测值取平均。

dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41),write = rep(c(30:70), 3))## store the predicted probabilities for each value of ses and write
pp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))## calculate the mean probabilities within each level of ses
by(pp.write[, 3:5], pp.write$ses, colMeans)

## pp.write$ses: high
## academic  general vocation
##    0.616    0.181    0.203
## ------------------------------------------------------
## pp.write$ses: low
## academic  general vocation
##    0.397    0.328    0.275
## ------------------------------------------------------
## pp.write$ses: middle
## academic  general vocation
##    0.426    0.201    0.373

有时候，一组图能过很好地表达大量的信息。使用我们上面的“pp.write” 对象，我们可以绘制在不同社会地位下，预测值与写作分数的关系的图。

## melt data set to long for ggplot2
lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")
head(lpp) # view first few rows

##   ses write variable probability
## 1 low    30 academic      0.0984
## 2 low    31 academic      0.1072
## 3 low    32 academic      0.1165
## 4 low    33 academic      0.1265
## 5 low    34 academic      0.1370
## 6 low    35 academic      0.1483

## plot predicted probabilities across write values for
## each level of ses facetted by program type
ggplot(lpp, aes(x = write, y = probability, colour = ses)) +geom_line() +facet_grid(variable ~ ., scales="free")