Hi, Indeed I fail to add the last portion of the code,
confint.par <- cbind(confint.par * qnorm(.975), confint.par *qnorm(.025)) I will try the delta method, I was hoping there was an easier way of doing it. Thanks Peter for your prompt answered. Cristián Montes Site Productivity Division Head Bioforest S.A. Chile -----Mensaje original----- De: Peter Dalgaard [mailto:pda...@gmail.com] Enviado el: Martes, 03 de Agosto de 2010 02:48 a.m. Para: Cristian Montes CC: r-help@r-project.org Asunto: Re: [R] Confidence Bands in nonlinear regression using optim and maximum likelihood 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@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....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@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.
