运行 mle2 函数时出错 (bbmle)

Error when running mle2 function (bbmle)

当 运行 R 中 bbmle 包中的 mle2() 函数时,我收到以下错误:

some parameters are on the boundary: variance-covariance calculations based on Hessian may be unreliable

我想了解这是因为我的数据有问题还是因为正确调用函数有问题。不幸的是,我无法 post 我的真实数据,所以我使用了相同样本大小的类似工作示例。

我使用的自定义 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

我不完全理解你的问题的背景,但是:

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

### 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 论点,即 returns 负对数似然当 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))
}