confidence interval around predicted value from complex inverse function

Question

I'm trying to get a 95% confidence interval around some predicted values, but am not capable of achieving this.

Basically, I estimated a growth curve like this:

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))

I have independent size measurements, for which I would like to predict the age. Therefore, the inverse of the function, which is not very straightforward, I calculated like this:

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

Although this works fine, I need to know the confidence interval (95%) around this estimate (x). For linear models there is for example "predict(... confidence=)". I could also bootstrap the function somehow to get the quantiles associated with the parameters (didn't find how), to then use the extremes of those to calculate the maximum and minimum values predictable. But that doesn't really look like the good way of doing this....

Any help would be greatly appreciated.

EDIT after answer:

So this worked (explained in the book of Ben Bolker, see answer):

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))

On the rather bad fake data I gave, this works of course rather horrible. But on the real data (which I have a problem recreating), this is ok!

Answer 1

Unless I'm mistaken, you're going to have to use either regular (parametric) bootstrapping or a method called either "population predictive intervals" (e.g., see section 5 of chapter 7 of Bolker 2008), which assumes that the sampling distributions of your parameters are multivariate Normal. However, I think you may have bigger problems, unless I've somehow messed up your model in adapting it ...

Generate data (note that random data may actually bad for testing your model - see below ...)

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 the data and the initial curve estimate:

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))

I had trouble with nls so I used minpack.lm::nls.lm, which is slightly more robust. (There are other options here, e.g. calculating the derivatives and providing the gradient function, or using AD Model Builder or Template Model Builder, or using the nls2 package.)

For nls.lm we need a function that returns the residuals:

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 

Here's our problem:

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)
})

With this kind of fit, the results of the inverse function are going to be extremely unstable, as the inverse function is many-to-one, with the number of inverse values depending sensitively on the parameter values ...

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

If I pick my interval very carefully I can get back the correct answer ...

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

Maybe your model will be better behaved with more realistic data. I hope so ...