R help - Aug 2010 - Confidence Bands in nonlinear regression using optim and maximum likelihood

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


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


Any help will be appreciated.

Cristián Montes
Site Productivity Division Head
Bioforest S.A.
Chile.

	[[alternative HTML version deleted]]

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