Multinomial Logit Model (MNL) 模型R语言nnet包multinom函数实现实例
来源:互联网 发布:做标准件网络销售需要 编辑:程序博客网 时间:2024/06/07 00:27
最近做项目涉及到要使用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.errorsz
## (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 testp <- (1 - pnorm(abs(z), 0, 1))*2p
## (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_, 我们可以写出我们的模型等式。
一种选择结果类别的概率与基准选择的概率的比被称为相对危险度(relative risk), 有时候只是为了描述回归参数,我们也叫做odds. 相对危险度是对线性方程等号右边部分的取幂,取幂后的回归系数就是自变量每单位变化的相对危险度。下面我们开对模型的系数取幂来观察是如何变化的。
## extract the coefficients from the model and exponentiateexp(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
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
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
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 writepp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))## calculate the mean probabilities within each level of sesby(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
## melt data set to long for ggplot2lpp <- 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 typeggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~ ., scales="free")
- Multinomial Logit Model (MNL) 模型R语言nnet包multinom函数实现实例
- multinomial logit model 多项 Logit 模型
- R语言利用nnet包训练神经网络模型
- 多项式事件模型(multinomial event model)
- R语言实现SOM(自组织映射)模型(三个函数包+代码)
- LDA主题模型参数推导-dirichlet-multinomial model posterior
- R 语言 RFM 模型实现
- R语言实现 logistic 模型
- R语言实现RMF模型
- logit和logistic模型
- 基于R中的神经网络包(nnet)做分类(高级篇)
- R语言中的模型修正函数update
- R语言 caret包 findCorrelation()函数用法
- R语言之函数与包
- RFM模型及R语言实现
- RFM模型及R语言实现
- KMV模型的R语言实现
- R语言 数据挖掘 R包与函数
- 二叉树前序、中序、后序遍历相互求法
- Leetcode 331. Verify Preorder Serialization of a Binary Tree (Medium) (cpp)
- 欢迎使用CSDN-markdown编辑器
- 欢迎使用CSDN-markdown编辑器
- 从人眼到色彩空间
- Multinomial Logit Model (MNL) 模型R语言nnet包multinom函数实现实例
- Windows消息对Edit控件的处理
- 用R做中文文本分析--用R进行文本挖掘与分析:分词、画词云
- Mac OS X,下载并安装ant
- 第40篇 WebRTC(三)
- UDP连接调用connect()函数
- 【20.73%】【CF 716B】Complete the Word
- 【物联网(IoT)开发】Arduino IDE的工具>开发板菜单中找到我的开发板型号怎么办?
- 读过的文章 2016-09-18