r - 比较 multinom 和 stan multi logit 回归

Question

我正在尝试使用 stan 来理解/重现调查结果.但是，我被困在某个地方。我是否使用了错误的 stan 模型？

library(nnet);library(rstan);library(dplyr);library(tidyr)

#set up data 
n <- 100
set.seed(1)
dat <- data.frame(DV=factor(sample(letters[1:5], n, replace=T)),
       x1 = rnorm(n, mean=4.5, sd=1.3),
       x2 = sample(c(1:5), prob=c(0.035, 0.167, 0.083, 0.415, 0.298)),
       x3 = sample(c(0,1), prob=c(.51, .49)),
       x4 = round(rnorm(n, mean=48, sd=15),0))

library(nnet)
f <- as.formula("DV ~ x1 + x2 + x3 + x4")
out <- multinom(f, data=dat)
#summary(out)

#Store output to compare later on:
out.multinom <- tidyr::gather(as.data.frame(coef(out)), values) %>%  
  mutate(option=rep(row.names(coef(out)), 5)) %>%
  mutate(coef= paste0(option, ":", values))%>%
  dplyr::select(-option, -values)%>%
  rename(multinom = value)

到现在为止还挺好。现在将估计值与 stan 中的估计值进行比较:

该模型是从 here 复制粘贴的:

stan_model <- "
 data {
 int K;
 int N;
 int D;
 int y[N];
 matrix[N, D] x;
 }
 parameters {
 matrix[D, K] beta;
 }
 model {
 matrix[N, K] x_beta = x * beta;

 to_vector(beta) ~ normal(0, 2);

 for (n in 1:N)
 y[n] ~ categorical_logit(x_beta[n]);
 }

"

rstan_options(auto_write = TRUE)
options(mc.cores = 4)

M <- model.matrix(f, dat)

#data for stan
datlist <- list(N=nrow(M),                     #nr of obs
                K=length(unique(dat[,1])),     #possible outcomes
                D=ncol(M),                     #dimension of predictor matrix
                x=M,                           #predictor matrix
                y=as.numeric(dat[,1]))

 #estimate model
 b.out <- stan(model_code=stan_model, 
          data=datlist,
          iter = 1000,
          chains = 4,
          seed = 12591,
          control = list(max_treedepth = 11))


 res <- summary(b.out, par="beta", probs=.5)$summary %>% as.data.frame

 #store
 out.stan <- data.frame(beta=rep(colnames(M), each=5), 
                   value.stan = res[,1]) %>%
  mutate(option=rep(levels(dat$DV), length.out=25))%>%
  mutate(coef= paste0(option, ":", beta))%>%
  dplyr::select(-option, -beta)

#compare
merge(out.multinom, out.stan, by="coef", all.y=T)


            coef     multinom   value.stan
1  a:(Intercept)           NA -0.532803345
2           a:x1           NA  0.017020378
3           a:x2           NA  0.227622393
4           a:x3           NA -0.001617129
5           a:x4           NA  0.011291841
6  b:(Intercept) -0.308050243 -0.794106266
7           b:x1  0.314860536  0.353267471
8           b:x2 -0.305248243 -0.094612371
9           b:x3 -0.181849471 -0.203225779
10          b:x4 -0.002589588  0.007308861
11 c:(Intercept)  1.241939293  0.391090113
12          c:x1 -0.265908390 -0.230931524
13          c:x2 -0.121426457  0.113370970
14          c:x3 -0.486321891 -0.496869965
15          c:x4  0.004659122  0.019329331
16 d:(Intercept)  1.655236959  0.767339620
17          d:x1 -0.332715090 -0.291936815
18          d:x2 -0.159596712  0.072145756
19          d:x3  0.132897149  0.145568093
20          d:x4  0.002003899  0.017336138
21 e:(Intercept)  0.970658209  0.253888841
22          e:x1 -0.057914885 -0.017299638
23          e:x2 -0.501888386 -0.306445142
24          e:x3  0.515320410  0.517075558
25          e:x4  0.009115124  0.023669295

由于 estimands 不同，有些事情不是它应该的样子。我几乎可以肯定，我只是简单地复制并粘贴 stan 就错过了一些东西。代码。它是什么？ (我不是在谈论 multinom 使用 a 作为引用类别的事实)。

问题不应该在于数据集的大小，因为我使用了更多数据(这个例子是为 stackoverflow 准备的)。

谢谢你的帮助。

Answer 1

从输出来看，它看起来像 multinom被调用的函数使用 K-1系数以使模型可识别。这需要 a系数隐式为零。如果减去 a来自其他系数的系数，b通过 e ，您会得到大致相同的结果。结果不会完全相同，因为 Stan 在这里为您提供了后验方法和 multinom可能是给你别的东西。有使用 K 的讨论与 K-1在 Stan 用户指南中对多项式进行编码，就在您在“可识别性”标题下链接的同一部分中。