复杂反函数预测值的置信区间

Question

我试图在某些预测值周围获得 95% 的置信区间，但我无法实现。

基本上，我估计的增长曲线是这样的：

set.seed(123)
dat=data.frame(size=rnorm(50,10,3),age=rnorm(50,5,2))
S <- function(t,ts,C,K) ((C*K)/(2*pi))*sin(2*pi*(t-ts))
sommers <- function(t,Linf,K,t0,ts,C)
  Linf*(1-exp(-K*(t-t0)-S(t,ts,C,K)+S(t0,ts,C,K)))
model <- nls(size~sommers(age,Linf,K,t0,ts,C),data=dat,
             start=list(Linf=10,K=4.7,t0=2.2,C=0.9,ts=0.1))

我有独立的尺寸测量值，我想预测年龄。所以函数的反函数，不是很直接，我是这样计算的：

model.out=coef(model)
S.out <- function(t) 
  ((model.out[[4]]*model.out[[2]])/(2*pi))*sin(2*pi*(t-model.out[[5]]))
sommers.out <- function(t) 
  model.out[[1]]*(1-exp(-model.out[[2]]*(t-model.out[[3]])-S.out(t)+S.out(model.out[[3]])))
inverse = function (f, lower = -100, upper = 100) {
  function (y) uniroot((function (x) f(x) - y), lower = lower, upper = upper)[1]
}
sommers.inverse = inverse(sommers.out, 0, 25)
x= sommers.inverse(10)  #this works with my complete dataset, but not with this fake one

虽然这工作正常，但我需要知道围绕此估计值 (x) 的置信区间 (95%)。对于线性模型，例如 "predict(... confidence=)"。我还可以 bootstrap 该函数以某种方式获取与参数关联的分位数（未找到方法），然后使用这些参数的极值来计算可预测的最大值和最小值。但这看起来并不是真正的好方法....

如有任何帮助，我们将不胜感激。

回答后编辑：

所以这有效（在 Ben Bolker 的书中有解释，请参阅答案）：

vmat = mvrnorm(1000, mu = coef(mfit), Sigma = vcov(mfit)) 
dist = numeric(1000) 
for (i in 1:1000) {dist[i] = sommers_inverse(9.938,vmat[i,])} 
quantile(dist, c(0.025, 0.975))

根据我提供的相当糟糕的假数据，这当然很糟糕。但是在真实数据上（我在重新创建时遇到问题），这没问题！

Answer 1

除非我弄错了，否则您将不得不使用常规（参数）自举或称为 "population predictive intervals" 的方法（例如，参见 chapter 7 of Bolker 2008 的第 5 节），这假设您的参数的抽样分布是多元正态分布。但是，我认为你可能有更大的问题，除非我在调整它时以某种方式搞砸了你的模型......

生成数据（注意随机数据实际上可能不利于测试您的模型 - 见下文...）

set.seed(123)
dat <- data.frame(size=rnorm(50,10,3),age=rnorm(50,5,2))
S <- function(t,ts,C,K) ((C*K)/(2*pi))*sin(2*pi*(t-ts))
sommers <- function(t,Linf,K,t0,ts,C)
    Linf*(1-exp(-K*(t-t0)-S(t,ts,C,K)+S(t0,ts,C,K)))

绘制数据和初始曲线估计值：

plot(size~age,data=dat,ylim=c(0,16))
agevec <- seq(0,10,length=1001)
lines(agevec,sommers(agevec,Linf=10,K=4.7,t0=2.2,ts=0.1,C=0.9))

我在使用 nls 时遇到了麻烦，所以我使用了 minpack.lm::nls.lm，它稍微稳健一些。（这里还有其他选项，例如计算导数和提供梯度函数，或使用 AD Model Builder 或 Template Model Builder，或使用 nls2 包。）

对于nls.lm我们需要一个函数returns残差:

sommers_fn <- function(par,dat) {
   with(c(as.list(par),dat),size-sommers(age,Linf,K,t0,ts,C))
}
library(minpack.lm)
mfit <- nls.lm(fn=sommers_fn,
           par=list(Linf=10,K=4.7,t0=2.2,C=0.9,ts=0.1),
       dat=dat)
coef(mfit)
##        Linf           K          t0           C          ts 
##  10.6540185   0.3466328   2.1675244 136.7164179   0.3627371

这是我们的问题：

plot(size~age,data=dat,ylim=c(0,16))
lines(agevec,sommers(agevec,Linf=10,K=4.7,t0=2.2,ts=0.1,C=0.9))
with(as.list(coef(mfit)), {
     lines(agevec,sommers(agevec,Linf,K,t0,ts,C),col=2)
     abline(v=t0,lty=2)
     abline(h=c(0,Linf),lty=2)
})

这种拟合，反函数的结果会极度不稳定，因为反函数是多对一的，反函数的个数值敏感地取决于参数值...

sommers_pred <- function(x,pars) {
    with(as.list(pars),sommers(x,Linf,K,t0,ts,C))
}
sommers_pred(6,coef(mfit))  ## s(6)=9.93

sommers_inverse <- function (y, pars, lower = -100, upper = 100) {
    uniroot(function(x) sommers_pred(x,pars) -y, c(lower, upper))$root
}
sommers_inverse(9.938, coef(mfit))  ## 0.28

如果我选择间隔非常，我可以得到正确答案...

sommers_inverse(9.938, coef(mfit), 5.5, 6.2)

也许您的模型在使用更真实的数据时会表现得更好。我希望如此...

复杂反函数预测值的置信区间

confidence interval around predicted value from complex inverse function

r

confidence-interval