我想在R中运行多项式logit，并使用了两个库nnet和mlogit，它们会产生不同的结果并报告不同类型的统计信息。我的问题是：nnet报告的系数和标准误差与mlogit报告的系数和标准误差之间差异的根源是什么？
我想使用Latex将结果报告到stargazer文件。这样做时，会有一个权衡的问题：
如果使用mlogit的结果，则会得到所需的统计信息，例如psuedo R平方，但是，输出为长格式（请参见下面的示例）。
如果我使用来自nnet的结果，则格式与预期的一样，但它会报告我不感兴趣的统计信息（例如AIC），但不包括psuedo R squared。
我想在使用mlogit时以nnet的格式报告stargazer的统计信息。
这是一个可重现的示例，提供三种选择：
library(mlogit)
df = data.frame(c(0,1,1,2,0,1,0), c(1,6,7,4,2,2,1), c(683,276,756,487,776,100,982))
colnames(df) <- c('y', 'col1', 'col2')
mydata = df
mldata <- mlogit.data(mydata, choice="y", shape="wide")
mlogit.model1 <- mlogit(y ~ 1| col1+col2, data=mldata)
编译时的tex输出是我认为不希望的“长格式”：
现在，使用
nnet：library(nnet)
mlogit.model2 = multinom(y ~ 1 + col1+col2, data=mydata)
stargazer(mlogit.model2)
给出tex输出：
这是我想要的“宽”格式。注意不同的系数和标准误差。
最佳答案
据我所知，有三个R包可以估算多项式逻辑回归模型：mlogit，nnet和globaltest（来自Bioconductor）。我在这里不考虑mnlogit包，它是mlogit的更快，更有效的实现。
以上所有软件包均使用不同的算法，对于较小的样本，这些算法会给出不同的结果。对于中等大小的样本（使用n <- 100尝试），这些差异将消失。
考虑以下来自James Keirstead's blog的数据生成过程：
n <- 40
set.seed(4321)
df1 <- data.frame(x1=runif(n,0,100), x2=runif(n,0,100))
df1 <- transform(df1, y=1+ifelse(100 - x1 - x2 + rnorm(n,sd=10) < 0, 0,
ifelse(100 - 2*x2 + rnorm(n,sd=10) < 0, 1, 2)))
str(df1)
'data.frame': 40 obs. of 3 variables:
$ x1: num 33.48 90.91 41.15 4.38 76.35 ...
$ x2: num 68.6 42.6 49.9 36.1 49.6 ...
$ y : num 1 1 3 3 1 1 1 1 3 3 ...
table(df1$y)
1 2 3
19 8 13
这三个程序包估计的模型参数分别为：
library(mlogit)
df2 <- mlogit.data(df1, choice="y", shape="wide")
mlogit.mod <- mlogit(y ~ 1 | x1+x2, data=df2)
(mlogit.cf <- coef(mlogit.mod))
2:(intercept) 3:(intercept) 2:x1 3:x1 2:x2 3:x2
42.7874653 80.9453734 -0.5158189 -0.6412020 -0.3972774 -1.0666809
#######
library(nnet)
nnet.mod <- multinom(y ~ x1 + x2, df1)
(nnet.cf <- coef(nnet.mod))
(Intercept) x1 x2
2 41.51697 -0.5005992 -0.3854199
3 77.57715 -0.6144179 -1.0213375
#######
library(globaltest)
glbtest.mod <- globaltest::mlogit(y ~ x1+x2, data=df1)
(cf <- glbtest.mod@coefficients)
1 2 3
(Intercept) -41.2442934 1.5431814 39.7011119
x1 0.3856738 -0.1301452 -0.2555285
x2 0.4879862 0.0907088 -0.5786950
mlogit的globaltest命令在不使用参考结果类别的情况下拟合模型，因此可以按以下方式计算常用参数：(glbtest.cf <- rbind(cf[,2]-cf[,1],cf[,3]-cf[,1]))
(Intercept) x1 x2
[1,] 42.78747 -0.5158190 -0.3972774
[2,] 80.94541 -0.6412023 -1.0666813
关于三个软件包中参数的估计，在
mlogit::mlogit中详细解释了here中使用的方法。在
nnet::multinom中，模型是一个没有隐藏层，没有偏置节点和softmax输出层的神经网络；在我们的例子中，有3个输入单元和3个输出单元：nnet:::summary.nnet(nnet.mod)
a 3-0-3 network with 12 weights
options were - skip-layer connections softmax modelling
b->o1 i1->o1 i2->o1 i3->o1
0.00 0.00 0.00 0.00
b->o2 i1->o2 i2->o2 i3->o2
0.00 41.52 -0.50 -0.39
b->o3 i1->o3 i2->o3 i3->o3
0.00 77.58 -0.61 -1.02
最大条件似然是
multinom中用于模型拟合的方法。多项式对数模型的参数在
globaltest::mlogit中使用最大似然估计，并与等效对数线性模型和泊松似然一起使用。该方法在here中描述。对于由
multinom估计的模型，可以轻松计算McFadden的伪R平方，如下所示：nnet.mod.loglik <- nnet:::logLik.multinom(nnet.mod)
nnet.mod0 <- multinom(y ~ 1, df1)
nnet.mod0.loglik <- nnet:::logLik.multinom(nnet.mod0)
(nnet.mod.mfr2 <- as.numeric(1 - nnet.mod.loglik/nnet.mod0.loglik))
[1] 0.8483931
此时，使用
stargazer为由mlogit::mlogit估计的模型生成报告，该报告与multinom的报告尽可能相似。基本思想是将multinom创建的对象中的估计系数和概率替换为mlogit的相应估计。# Substitution of coefficients
nnet.mod2 <- nnet.mod
cf <- matrix(nnet.mod2$wts, nrow=4)
cf[2:nrow(cf), 2:ncol(cf)] <- t(matrix(mlogit.cf,nrow=2))
# Substitution of probabilities
nnet.mod2$wts <- c(cf)
nnet.mod2$fitted.values <- mlogit.mod$probabilities
结果如下：
library(stargazer)
stargazer(nnet.mod2, type="text")
==============================================
Dependent variable:
----------------------------
2 3
(1) (2)
----------------------------------------------
x1 -0.516** -0.641**
(0.212) (0.305)
x2 -0.397** -1.067**
(0.176) (0.519)
Constant 42.787** 80.945**
(18.282) (38.161)
----------------------------------------------
Akaike Inf. Crit. 24.623 24.623
==============================================
Note: *p<0.1; **p<0.05; ***p<0.01
现在，我正在处理最后一个问题：如何在上述
stargazer输出中可视化loglik，伪R2和其他信息。关于r - R中的多项式logit:mlogit与nnet，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/43623076/