R中的多项式logit : mlogit与nnet

Question

我想在R中运行一个多项式logit，并且使用了两个库nnet和mlogit，它们产生不同的结果并报告不同类型的统计数据。我的问题是：

1. nnet报告的系数和标准误差与mlogit报告的标准误差之间差异的原因是什么？
2. 我想使用Latex向stargazer文件报告我的结果。当这样做时，有一个问题的权衡：

- If I use the results from `mlogit` then I get the statistics I wish, such as psuedo R squared, however, the output is in long format (see example below). 
- If I use the results from `nnet` then the format is as expected, but it reports statistics that I am not interested in such as AIC, but does not include, for example, 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输出：

​

​

这是我想要的“宽”格式。注意不同的系数和标准误差。

Answer 1

据我所知，有三个R包允许估计多项式logistic回归模型：mlogit，nnet和globaltest (来自生物导体)。在这里，我不认为mnlogit包是一个更快、更有效的mlogit实现。

所有上述软件包都使用不同的算法，对于小样本，给出不同的结果。对于适度的样本大小，这些差异消失了(尝试使用n <- 100)。

考虑从詹姆斯·基尔斯特德的博客中获取的以下数据生成过程

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中使用的这里方法。

在nnet::multinom中，该模型是一个没有隐藏层、没有偏置节点和一个软最大输出层的神经网络；在我们的例子中，有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中用于模型拟合的方法。

用最大似然估计多项式logit模型的参数，并利用等效对数线性模型和泊松似然估计多项式logit模型的参数。该方法被描述为这里。

对于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

现在我正在研究最后一个问题:如何在上面的R2输出中可视化loglik、伪stargazer和其他信息。

Answer 2

如果您使用的是天文望远镜，您可以使用omit删除不需要的行或引用。下面是一个快速的例子，希望它能为你指明正确的方向。

注：我的假设是，您正在使用Rstudio和rmarkdown与针织品。

```{r, echo=FALSE}

图书馆(Mlogit)

df =data.frame(c(0，1，1，2，0)，c(1,6,7,4,2,2,1)，c(683,276,756,487,776,100,982))

(Df) <- c('y'，'col1'，'col2')

mydata = df

mldata <- mlogit.data(mydata，mldata= "y"，shape="wide")

mlogit.model1 <- mlogit(y ~1\ col1+col2，data=mldata)

mlogit.col1 <- mlogit(y ~1\x{e76f} col1，data = mldata)

mlogit.col2 <- mlogit(y ~1\x{e76f} col2，data = mldata)

# MLOGIT

```{r echo = FALSE, message = TRUE, error = TRUE, warning = FALSE, results = 'asis'}

图书馆(望星者)

天文望远镜(mlogit.model1，type = "html")

观星者(mlogit.col1，

      mlogit.col2,

      type = "html",

      omit=c("1:col1","2:col1","1:col2","2:col2"))

结果：

​

​

注意，第二个图像省略了1:col1，2:col2，1:col2和2:col2。

​

​