运行 mle2 函数时出错 (bbmle)

Question

从 R 中的 bbmle 包运行 mle2() 函数时收到以下错误：

一些参数在边界上：基于 Hessian 的方差-协方差计算可能不可靠

我正在尝试了解这是否是由于我的数据问题或正确调用函数的问题。不幸的是，我无法发布我的真实数据，因此我使用了相同样本量的类似工作示例。

我使用的自定义dAction 函数是一个softmax 函数。优化必须有上限和下限，所以我使用的是 L-BFGS-B 方法。

library(bbmle)
set.seed(3939)

### Reproducible data
dat1 <- rnorm(30, mean = 3, sd = 1)
dat2 <- rnorm(30, mean = 3, sd = 1)
dat1[c(1:3, 5:14, 19)] <- 0
dat2[c(4, 15:18, 20:22, 24:30)] <- 0

### Data variables
x <- sample(1:12, 30, replace = TRUE)
pe <- dat1
ne <- dat2

### Likelihood
dAction <- function(x, a, b, t, pe, ne, log = FALSE) {
  u <- exp(((x - (a * ne) - (b * pe)) / t))
  prob <- u / (1 + u)

  if(log) return(prob) else return(-sum(log(prob)))
}

### Fit
fit <- mle2(dAction,
            start = list(a = 0.1, b = 0.1, t = 0.1),
            data = list(x = x, pe = pe, ne = ne),
            method = "L-BFGS-B",
            lower = c(a = 0.1, b = 0.1, t = 0.1),
            upper = c(a = 10, b = 1, t = 10))

Warning message:
In mle2(dAction, start = list(a = 0.1, b = 0.1, t = 0.1), data = list(x = x,  :
  some parameters are on the boundary: variance-covariance calculations based on Hessian may be unreliable

这是summary() 的结果：

summary(fit)
Maximum likelihood estimation

Call:
mle2(minuslogl = dAction, start = list(a = 0.1, b = 0.1, t = 0.1), 
    method = "L-BFGS-B", data = list(x = x, pe = pe, ne = ne), 
    lower = c(a = 0.1, b = 0.1, t = 0.1), upper = c(a = 10, b = 1, 
        t = 10))

Coefficients:
  Estimate Std. Error z value Pr(z)
a      0.1         NA      NA    NA
b      0.1         NA      NA    NA
t      0.1         NA      NA    NA

-2 log L: 0.002048047 

Warning message:
In sqrt(diag(object@vcov)) : NaNs produced

以及置信区间的结果

confint(fit)
Profiling...

  2.5 %    97.5 %
a    NA 1.0465358
b    NA 0.5258828
t    NA 1.1013322

Warning messages:
1: In sqrt(diag(object@vcov)) : NaNs produced
2: In .local(fitted, ...) :
  Non-positive-definite Hessian, attempting initial std err estimate from diagonals

Answer 1

我不完全了解您的问题的背景，但是：

这个问题（是否是一个真正的问题在很大程度上取决于我不理解的上述上下文）与您的限制有关。如果我们在没有约束的情况下进行拟合：

### Fit
fit <- mle2(dAction,
            start = list(a = 0.1, b = 0.1, t = 0.1),
            data = list(x = x, pe = pe, ne = ne))
            ## method = "L-BFGS-B",
            ## lower = c(a = 0.1, b = 0.1, t = 0.1),
            ## upper = c(a = 10, b = 1, t = 10))

我们得到的系数低于您的界限。

coef(fit)
         a          b          t 
0.09629301 0.07724332 0.02405173 

如果这是正确的，则至少有一个约束将处于活动状态（即，当我们符合下限时，我们的参数中至少有一个会达到边界 - 事实上，这就是全部）。当拟合在边界上时，用于计算置信区间（Wald 区间）的最简单机器不起作用。但是，这不会影响您在上面报告的配置文件置信区间估计值。这些是正确的 - 下限报告为 NA，因为置信下限位于边界处（如果您愿意，可以将它们替换为 0.1）。

如果您没想到最佳拟合会出现在边界上，那么我不知道发生了什么，可能是数据问题。

你的对数似然函数没有错，但它有点令人困惑，因为你有一个 log 参数，它在 log=FALSE（默认）和 log=TRUE 时返回负对数似然。在我意识到这一点之前，我重写了函数（我还尽可能在对数尺度上进行计算，使其在数值上更加稳定）。

dAction <- function(x, a, b, t, pe, ne) {
  logu <- (x - (a * ne) - (b * pe)) / t
  lprob <- logu - log1p(exp(logu))
  return(-sum(lprob))
}