0

我是一名使用最大似然法在 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 上得到了一个负界限，这是不可能的，因为它实际上代表了生物模型中的一个平方数以及一些警告。

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

我看到发生的是，巧合的是，我的粗麻布矩阵对于优化空间中的某些值是奇异的。但是，由于我正在优化 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

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

4

1 回答 1

6

您可以更改参数化，以便始终满足约束条件。将可能性重写为 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)))
 }
于 2010-04-27T08:38:10.887 回答