[R] Confidence Bands in nonlinear regression using optim and maximum likelihood

[Peter Dalgaard]

Cristian Montes wrote:
>  Hello,
> 
> I am trying to plot confidence bands on the mean and prediction bands for the following
> nonlinear regression, using maximum likelihood via optim.  A toy example with data and
> code of what I am trying to accomplish is:
> 
> VOL<-c(0.01591475, 1.19147935 ,6.34102460, 53.68809287, 91.90143074, 116.21397007,
> 146.41843056, 215.64535337, 256.53149673, 315.73609232)
> Age <-c(1.622222,  2.833333 , 3.927778,  7.150000,  8.166667,
> 8.752778 , 9.444444, 10.797222 ,11.725000, 12.775000)
> 
> with that I fit the following likelihood function
> 
> loglik <- function(param, x, y)
>     {
>     b0     <- param[1]
>     b1     <- param[2]
>     sigma  <- param[3]
> 
>     res    <- dnorm(y, b0*exp(b1/x), sigma, log= T)
>     .value <- sum(res)
>     }
> 
> 
> regression <- optim(par     = c(800,-2,.2),
>                     fn      = loglik,
>                     method  = "L-BFGS-B",
>                     lower   = c(-Inf, -Inf, 1e-10),
>                     control = list(fnscale = -1),
>                     x       = Age,
>                     y       = VOL,
>                     hessian = T)
> 
> confint.par  <- sqrt(diag(solve(regression$hessian *-1)))

Um, that's SE, not confidence intervals.

> 
> now that I have confidence intervals for each parameter, how can a draw confidence
> bands for the mean regression line?

You can use the delta method. You need to establish the gradient of the
expected mean with respect to the parameters, then pre- and
post-multiply it (as in g'Vg) onto the variance-covariance matrix of the
parameters.

All of this is based on linear approximation theory. Doing it via direct
likelihood profiling is considerably harder.


> 
> 
> Any help will be appreciated.
> 
> Cristián Montes
> Site Productivity Division Head
> Bioforest S.A.
> Chile.
> 
> 	[[alternative HTML version deleted]]
> 
> 
> 
> ------------------------------------------------------------------------
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com