r - confidence interval around predicted value from complex inverse function -

i'm trying 95% confidence interval around predicted values, not capable of achieving this.

basically, estimated growth curve 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, predict age. therefore, inverse of function, not straightforward, calculated 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 complete dataset, not fake 1 

although works fine, need know confidence interval (95%) around estimate (x). linear models there example "predict(... confidence=)". bootstrap function somehow quantiles associated parameters (didn't find how), use extremes of calculate maximum , minimum values predictable. doesn't way of doing this....

any appreciated.

edit after answer:

so worked (explained in book of ben bolker, see answer):

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

on rather bad fake data gave, works of course rather horrible. on real data (which have problem recreating), ok!

unless i'm mistaken, you're going have use either regular (parametric) bootstrapping or method called either "population predictive intervals" (e.g., see section 5 of chapter 7 of bolker 2008), assumes sampling distributions of parameters multivariate normal. however, think may have bigger problems, unless i've somehow messed model in adapting ...

generate data (note random data may bad testing 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 data , 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 nls used minpack.lm::nls.lm, more robust. (there other options here, e.g. calculating derivatives , providing gradient function, or using ad model builder or template model builder, or using nls2 package.)

for nls.lm need function returns 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 kind of fit, results of inverse function going extremely unstable, inverse function many-to-one, number of inverse values depending sensitively on 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 pick interval very can correct answer ...

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

maybe model better behaved more realistic data. hope ...