mlogit 包无法恢复合成混合 logit 模型

机器算法验证 r 混合模式 mlogit
2022-04-10 14:27:21

我正在从以下合成混合效应模型生成数据，用于代理选择运输模式的效用：hj

Uhj=βpriceXh,pricej+βtimeXh,timej+γhXh,busj+ϵhjγh∼N(0,10)

其中的代理的价格和时间。是总线，则等于 1 。Xh,pricej,Xh,timejhjXh,busjj

真正的模型有。βprice=βtime=−1

对于每个代理，我为每种交通方式（红色巴士、火车、汽车）生成价格和时间协变量（）和中对每个代理的随机效应进行采样，然后从相应的 logit 模型中对运输选择进行采样。Xh,pricej,Xh,timejγhN(0,10)

这是生成选择数据的 R 代码：

library(mlogit)
set.seed(1987)
n.samples <- 10000
re.sd     <- 10 # random effect standard deviation

## generate price, time covariates
mode       <- c('red.bus', 'train', 'car')
is.bus     <- c(1, 0, 0)
price.data <- rnorm(3*n.samples) + c(2,5,10)
time.data  <- rnorm(3*n.samples) + c(10, 5, 2)
obs          <- data.frame(mode=mode,
                           price=price.data,
                           time=time.data,
                           bus=is.bus)


## generate random effect + choice data
obs$agent.id <- rep(1:n.samples, each=3)
obs$choice   <- NA
obs$util     <- NA
for( a.id in unique(obs$agent.id) ) {
    xx.sub <- obs[obs$agent.id == a.id,]
    obs$util[obs$agent.id == a.id]   <- xx.sub$price*-1 + xx.sub$time*-1 + xx.sub$bus*rnorm(1,sd=re.sd)
    uu <- obs$util[obs$agent.id == a.id]
    p.vec <- exp(uu)/sum(exp(uu))
    obs$choice[obs$agent.id == a.id] <- rmultinom(1, 1, p.vec)==1
}


logit.data <- mlogit.data(obs,
                          shape   = "long",
                          choice  = "choice",
                          varying = which( colnames(obs) %in% c('price', 'time', 'bus') ),
                          alt.var = 'mode')

以下是生成的数据集的前几行：

> head(logit.data)
             mode      price     time bus agent.id choice       util
1.red.bus red.bus  0.7113747 9.700200   1        1  FALSE -21.306406
1.train     train  5.2730015 5.908244   0        1   TRUE -11.181246
1.car         car 12.8485015 2.256558   0        1  FALSE -15.105060
2.red.bus red.bus  1.7021472 9.642485   1        2   TRUE  -7.677141
2.train     train  5.2443204 5.671528   0        2  FALSE -10.915848
2.car         car 10.2686018 2.250377   0        2  FALSE -12.518979

我试图用mlogit包拟合正确指定的模型，但我发现估计是错误的（特别是随机效应的标准偏差）：

m.mixed <- mlogit(choice ~ price + time + bus | 0,
                  data=logit.data,
                  rpar= c(bus = 'n'),
                  R = 300, halton = NA)

summary(m.mixed)
> summary(m.mixed)

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data, 
    rpar = c(bus = "n"), R = 300, halton = NA)

Frequencies of alternatives:
    car red.bus   train 
 0.1317  0.4084  0.4599 

bfgs method
8 iterations, 0h:1m:55s 
g'(-H)^-1g =   594 
last step couldn't find higher value 

Coefficients :
        Estimate Std. Error  t-value  Pr(>|t|)    
price  -1.253720   0.026698 -46.9592 < 2.2e-16 ***
time   -1.329888   0.032389 -41.0603 < 2.2e-16 ***
bus     1.402522   0.318802   4.3993 1.086e-05 ***
sd.bus 19.455180   2.107521   9.2313 < 2.2e-16 ***

真实sd.bus值为 10，但mlogit估计接近 20。

为什么 mlogit 不能恢复真实模型？

2个回答

您似乎遇到了一个不幸的优化参数组合，特别是关于 Halton 伪随机序列，它可能有一个错误。BFGS 似乎会提前停止R = 300，但不会出现其他明显更小或更大的值。幸运的是，在这种情况下，您不需要大的值R（或 Halton）。

在我最初运行您的代码时，我得到了与您相同的结果，运行时统计数据表明 BFGS 需要 8 次迭代才能收敛。然后我R将函数调用更改为等于 30：

> m.mixed <- mlogit(choice ~ price + time + bus | 0,
+                   data=logit.data,
+                   rpar= c(bus = 'n'),
+                   R = 30, halton = NA)
> 
> summary(m.mixed)

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data, 
    rpar = c(bus = "n"), R = 30, halton = NA)

Frequencies of alternatives:
    car red.bus   train 
 0.1317  0.4084  0.4599 

bfgs method
22 iterations, 0h:1m:28s 
g'(-H)^-1g = 3.83E-07 
gradient close to zero 

Coefficients :
        Estimate Std. Error  z-value Pr(>|z|)    
price  -0.988457   0.026861 -36.7989   <2e-16 ***
time   -0.990255   0.032661 -30.3195   <2e-16 ***
bus    -0.118121   0.227826  -0.5185   0.6041    
sd.bus 10.369252   0.846377  12.2513   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -8743.9

请注意，不仅系数估计值正确，而且考虑到它们的标准误差，在合理接近实际值的意义上，运行时间是R = 300案例的 40%，尽管需要 22 次 BFGS 迭代而不是 8 次。

还需要 22 次迭代，运行时间增加了R = 10030% 多一点，达到初始情况所需的一半，并且系数估计与运行中的基本相同R = 30：

Coefficients :
        Estimate Std. Error  z-value Pr(>|z|)    
price  -0.989468   0.026877 -36.8144   <2e-16 ***
time   -0.991518   0.032679 -30.3410   <2e-16 ***
bus    -0.108744   0.226858  -0.4793   0.6317    
sd.bus 10.316022   0.841978  12.2521   <2e-16 ***

完全避免默认的 Halton 序列也解决了这个问题：

> m.mixed <- mlogit(choice ~ price + time + bus | 0,
+                   data=logit.data,
+                   rpar= c(bus = 'n'),
+                   R = 300, halton = NULL)
> 
> summary(m.mixed)

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data, 
    rpar = c(bus = "n"), R = 300, halton = NULL)

 *** blah blah blah ***     

Coefficients :
        Estimate Std. Error  z-value Pr(>|z|)    
price  -0.988426   0.026859 -36.8006   <2e-16 ***
time   -0.990478   0.032655 -30.3314   <2e-16 ***
bus    -0.143370   0.229167  -0.6256   0.5316    
sd.bus 10.382572   0.847424  12.2519   <2e-16 ***

这种情况下的运行时间仍然比 Halton 情况下低 10%，并且 BFGS 需要 23 次迭代才能收敛。8 iterations在这种情况下，这R = 300绝对是一个异常值。

然而，设置R = 400并生成“正确”结果再次打破了估计，BFGS 需要 12 次迭代，估计为.halton = NAR = 299sd.bus = 15.248...

编辑：

我还尝试了具有 4000 个样本的不同种子，但R = 300仍然halton = NA产生了不好的结果，甚至比原来的情况更糟。重新参数化调用以指定 Halton prime 和 drop 参数会产生不稳定的结果；prime=11工作，但prime=29惨遭失败R=300。然后我浏览了mlogitR 代码（感谢您找到它，@khoda！），但 Halton 序列代码可以正常工作。

多个其他测试，结合评论中提到的 OP 测试，让我得出结论，Halton 在一维情况下不能始终如一地工作，至少对于这个问题。在实际实践中，当真实参数无法与估计值进行比较时，有必要在调用中尝试几种不同的参数化，halton并检查结果的一致性（以及对数似然的值，我怀疑） . 完全避免 Halton 并为随机数生成器指定一个递增的值序列，直到获得稳定的估计是一种替代方案，除了运行时考虑之外，它也可能是可行的。RmlogitR

我相信 BFGS 的实施是罪魁祸首。我的前两个线索是：

  1. mlogit()使用参数method='bhhh'而不是默认值调用会bfgs产生更准确的估计。
  2. 当我从 获得不准确的估计时bfgs，优化器的停止条件是last step couldn't find higher value，这表明 BFGS 步骤不是上升方向。

我遵循了L-BFGS-B FORTRAN SUBROUTINES FOR LARGE BOUND BOUND CONSTRAINED OPTIMIZATION中的方法，我在此引用：

如果线搜索在对目标函数进行 20 次评估后无法找到目标值足够低的点，我们会得出结论，当前方向没有用。 在这种情况下，所有校正向量都被丢弃，并沿最陡下降方向重新开始迭代

我更新了mlogit.optim()函数中的 ，mlogit/R/mlogit.tools.R以便如果在线搜索中未找到上升步骤，则将逆 Hessian 的 BFGS 近似重置为恒等式。我将最大重置次数限制为 10（在这些测试中，我从未达到最大值）。更新如下所示（此复制/粘贴中的第一行未更改）：

    # eval the function and compute the gradient and the hessian
    x <- eval(f, parent.frame())
    if (is.null(x)){
        if(method == 'bfgs' && num.bfgs.reset < 10) {
            num.bfgs.reset <- num.bfgs.reset + 1
            Hm1            <- diag(nrow(Hm1))
            x              <- oldx
            next # try again
        } else {
            ## x is null if steptol is reached
            code = 3
            break
        }
    }

运行与以前相同的模拟，我得到了更准确的标准偏差估计：

Call:
mlogit(formula = choice ~ price + time + bus | 0, data = logit.data, 
    rpar = c(bus = "n"), R = 300, halton = NA, method = "bfgs")

Frequencies of alternatives:
    car red.bus   train 
 0.1317  0.4084  0.4599 

bfgs method
21 iterations, 0h:3m:50s 
g'(-H)^-1g = 1.8E-06 
successive function values within tolerance limits 

Coefficients :
        Estimate Std. Error  z-value Pr(>|z|)    
price  -0.989864   0.026878 -36.8276   <2e-16 ***
time   -0.991994   0.032682 -30.3530   <2e-16 ***
bus    -0.118717   0.228376  -0.5198   0.6032    
sd.bus 10.357554   0.847985  12.2143   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -8743.2

random coefficients
    Min.   1st Qu.     Median       Mean  3rd Qu. Max.
bus -Inf -7.104781 -0.1187173 -0.1187173 6.867347  Inf

如果我更改R、生成新数据集等，我会得到类似的结果。

其它你可能感兴趣的问题
上一篇PMF 的值可以大于 1 吗？ 下一篇统计中的符号（参数/估计器/估计）