如何在不反转 R 中的奇异 Hessian 矩阵的情况下获得置信区间？答案

【问题标题】：How do I get confidence intervals without inverting a singular Hessian matrix in R?如何在不反转 R 中的奇异 Hessian 矩阵的情况下获得置信区间？
【发布时间】：2011-02-10 17:22:36
【问题描述】：

我是一名使用最大似然法在 R 中研究流行病学模型的学生。我创建了负对数似然函数。这看起来有点恶心，但它是：

NLLdiff = function(v1, CV1, v2, CV2, st1 = (czI01 - czV01), st2 = (czI02 - czV02), st01 = czI01, st02 = czI02, tt1 = czT01, tt2 = czT02) { 
    prob1 = (1 + v1 * CV1 * tt1)^(-1/CV1)
    prob2 = ( 1 + v2 * CV2 * tt2)^(-1/CV2) 
    -(sum(dbinom(st1, st01, prob1, log = T)) + sum(dbinom(st2, st02, prob2, log = T)))
 }

第一行看起来如此糟糕的原因是因为它需要的大部分数据都是在那里输入的。例如，czI01 已经声明。我这样做只是为了让以后对该函数的调用不必都包含糟糕的向量。

然后我使用 mle2（库 bbmle）针对 CV1、CV2、v1 和 v2 进行了优化。这也有点恶心，看起来像：

ml.cz.diff = mle2 (NLLdiff, start=list(v1 = vguess, CV1 = cguess, v2 = vguess, CV2 = cguess), method="L-BFGS-B", lower = 0.0001)

现在，到这里为止一切正常。 ml.cz.diff 为我提供了值，我可以将其变成合理地适合我的数据的图。我也有几个不同的模型，可以得到 AICc 值来比较它们。但是，当我尝试获得关于 v1、CV1、v2 和 CV2 的置信区间时，我遇到了问题。基本上，我在 CV1 上得到了一个负界限，这是不可能的，因为它实际上代表了生物模型中的一个平方数以及一些警告。

有没有更好的方法来获得置信区间？或者，真的，一种获得在这里有意义的置信区间的方法？

我看到的情况是，巧合的是，我的 hessian 矩阵对于优化空间中的某些值是奇异的。但是，由于我正在优化 4 个变量并且没有过于广泛的编程知识，所以我无法想出一个不依赖于 hessian 的优化方法。我已经用谷歌搜索了这个问题 - 它表明我的模型很糟糕，但我正在重建之前完成的一些工作，这表明我的模型真的并不糟糕（我使用 ml.cz.diff 制作的图看起来像原始作品的图）。我还阅读了手册的相关部分以及 Bolker 的书R 中的生态模型。我还尝试了不同的优化方法，这导致运行时间更长但错误相同。 “SANN”方法没有在一个小时内完成运行，所以我没有等着看结果。

简而言之：我的置信区间很差。有没有一种相对简单的方法可以在 R 中修复它们？

我的向量是：

czT01 = c(5, 5, 5, 5, 5, 5, 5, 25, 25, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 50, 50)
czT02 = c(5, 5, 5, 5, 5, 10, 10, 10, 10, 10, 25, 25, 25, 25, 25, 50, 50, 50, 50, 50, 75, 75, 75, 75, 75)
czI01 = c(25, 24, 22, 22, 26, 23, 25, 25, 25, 23, 25, 18, 21, 24, 22, 23, 25, 23, 25, 25, 25)
czI02 = c(13, 16, 5, 18, 16, 13, 17, 22, 13, 15, 15, 22, 12, 12, 13, 13, 11, 19, 21, 13, 21, 18, 16, 15, 11)
czV01 = c(1, 4, 5, 5, 2, 3, 4, 11, 8, 1, 11, 12, 10, 16, 5, 15, 18, 12, 23, 13, 22)
czV02 = c(0, 3, 1, 5, 1, 6, 3, 4, 7, 12, 2, 8, 8, 5, 3, 6, 4, 6, 11, 5, 11, 1, 13, 9, 7)

我的猜测是：

v = -log((c(czI01, czI02) - c(czV01, czV02))/c(czI01, czI02))/c(czT01, czT02)
vguess = mean(v)
cguess = var(v)/vguess^2

也有可能我做的其他事情完全错误，但我的结果似乎合理，所以我没有抓住它。

【问题讨论】：

标签： r statistics confidence-interval

【解决方案1】：

您可以更改参数化，以便始终满足约束条件。将可能性重写为 ln(CV1) 和 ln(CV2) 的函数，这样您就可以确定 CV1 和 CV2 严格为正。

NLLdiff_2 = function(v1, lnCV1, v2, lnCV2, st1 = (czI01 - czV01), st2 = (czI02 - czV02), st01 = czI01, st02 = czI02, tt1 = czT01, tt2 = czT02) { 
prob1 = (1 + v1 * exp(lnCV1) * tt1)^(-1/exp(lnCV1))
prob2 = ( 1 + v2 * exp(lnCV2) * tt2)^(-1/exp(lnCV2)) 
-(sum(dbinom(st1, st01, prob1, log = T)) + sum(dbinom(st2, st02, prob2, log = T)))
 }

【讨论】：