mlogit 的 vglm() 和 multinomial() 结果的巨大差异
Huge difference in result of vglm() and multinomial() for mlogit
我正在为 iris
数据集做多项逻辑回归模型,
library(VGAM)
mlog1 <- vglm(Species ~ ., data=iris, family=multinomial())
coef(mlog1)
系数为:
(Intercept):1 (Intercept):2 Sepal.Length:1 Sepal.Length:2 Sepal.Width:1
34.243397 42.637804 10.746723 2.465220 12.815353
Sepal.Width:2 Petal.Length:1 Petal.Length:2 Petal.Width:1 Petal.Width:2
6.680887 -25.042636 -9.429385 -36.060294 -18.286137
然后我使用 multinom()
函数并做同样的事情:
library(nnet)
mlog2 <- multinom(Species ~ ., data=iris)
系数:
Coefficients:
(Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
versicolor 18.69037 -5.458424 -8.707401 14.24477 -3.097684
virginica -23.83628 -7.923634 -15.370769 23.65978 15.135301
这两个结果好像差距很大?我哪里做错了?我怎样才能修复它们并获得相似的结果?
差距由两个因素造成:(1)VGAM
中的multinomial()
族默认选择参考作为响应因子的最后一级,而multinom()
中的nnet
始终以第一层为参考。 (2) 鸢尾花数据中的物种类别可以线性分离,因此导致系数非常大,标准误差也很大。当对数似然几乎没有进一步改变时,数值优化究竟停止在哪里,在不同的实现中可能会有所不同,但实际上是无关紧要的。
作为一个没有分离的例子,考虑一个学校选择回归模型,该模型基于 AER
包中德国社会经济小组(1994-2002)的数据:
data("GSOEP9402", package = "AER")
library("nnet")
m1 <- multinom(school ~ meducation + memployment + log(income) + log(size),
data = GSOEP9402)
m2 <- vglm(school ~ meducation + memployment + log(income) + log(size),
data = GSOEP9402, family = multinomial(refLevel = 1))
然后,两个模型导致基本相同的系数:
coef(m1)
## (Intercept) meducation memploymentparttime memploymentnone
## Realschule -6.366449 0.3232377 0.4422277 0.7322972
## Gymnasium -22.476933 0.6664295 0.8964440 1.0581122
## log(income) log(size)
## Realschule 0.3877988 -1.297537
## Gymnasium 1.5347946 -1.757441
coef(m2, matrix = TRUE)
## log(mu[,2]/mu[,1]) log(mu[,3]/mu[,1])
## (Intercept) -6.3666257 -22.4778081
## meducation 0.3232500 0.6664550
## memploymentparttime 0.4422720 0.8964986
## memploymentnone 0.7323156 1.0581625
## log(income) 0.3877985 1.5348495
## log(size) -1.2975203 -1.7574912
我正在为 iris
数据集做多项逻辑回归模型,
library(VGAM)
mlog1 <- vglm(Species ~ ., data=iris, family=multinomial())
coef(mlog1)
系数为:
(Intercept):1 (Intercept):2 Sepal.Length:1 Sepal.Length:2 Sepal.Width:1
34.243397 42.637804 10.746723 2.465220 12.815353
Sepal.Width:2 Petal.Length:1 Petal.Length:2 Petal.Width:1 Petal.Width:2
6.680887 -25.042636 -9.429385 -36.060294 -18.286137
然后我使用 multinom()
函数并做同样的事情:
library(nnet)
mlog2 <- multinom(Species ~ ., data=iris)
系数:
Coefficients:
(Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
versicolor 18.69037 -5.458424 -8.707401 14.24477 -3.097684
virginica -23.83628 -7.923634 -15.370769 23.65978 15.135301
这两个结果好像差距很大?我哪里做错了?我怎样才能修复它们并获得相似的结果?
差距由两个因素造成:(1)VGAM
中的multinomial()
族默认选择参考作为响应因子的最后一级,而multinom()
中的nnet
始终以第一层为参考。 (2) 鸢尾花数据中的物种类别可以线性分离,因此导致系数非常大,标准误差也很大。当对数似然几乎没有进一步改变时,数值优化究竟停止在哪里,在不同的实现中可能会有所不同,但实际上是无关紧要的。
作为一个没有分离的例子,考虑一个学校选择回归模型,该模型基于 AER
包中德国社会经济小组(1994-2002)的数据:
data("GSOEP9402", package = "AER")
library("nnet")
m1 <- multinom(school ~ meducation + memployment + log(income) + log(size),
data = GSOEP9402)
m2 <- vglm(school ~ meducation + memployment + log(income) + log(size),
data = GSOEP9402, family = multinomial(refLevel = 1))
然后,两个模型导致基本相同的系数:
coef(m1)
## (Intercept) meducation memploymentparttime memploymentnone
## Realschule -6.366449 0.3232377 0.4422277 0.7322972
## Gymnasium -22.476933 0.6664295 0.8964440 1.0581122
## log(income) log(size)
## Realschule 0.3877988 -1.297537
## Gymnasium 1.5347946 -1.757441
coef(m2, matrix = TRUE)
## log(mu[,2]/mu[,1]) log(mu[,3]/mu[,1])
## (Intercept) -6.3666257 -22.4778081
## meducation 0.3232500 0.6664550
## memploymentparttime 0.4422720 0.8964986
## memploymentnone 0.7323156 1.0581625
## log(income) 0.3877985 1.5348495
## log(size) -1.2975203 -1.7574912