I'm trying to get confidence intervals on my predictions for a non-linear mixed model in nlme, using resampling of parameter values. I got a result but would like to know if I'm going in the right direction.

I here use the 'Loblolly' data from the nlme help to illustrate what I have done.

library(effects) library(nlme) library(MASS) fm1 <- nlme(height ~ SSasymp(age, Asym, R0, lrc), data = Loblolly, fixed = Asym + R0 + lrc ~ 1, random = Asym ~ 1, start = c(Asym = 103, R0 = -8.5, lrc = -3.3)) xvals=seq(min(Loblolly$age),max(Loblolly$age),length.out=100) nresamp=1000 pars.picked = mvrnorm(nresamp, mu = fixef(fm1), Sigma = vcov(fm1)) # pick new parameter values by sampling from multivariate normal distribution based on fit yvals = matrix(0, nrow = nresamp, ncol = length(xvals)) for (i in 1:nresamp) { yvals[i,] = sapply(xvals,function (x) SSasymp(x,pars.picked[i,1], pars.picked[i,2], pars.picked[i,3])) } quant = function(col) quantile(col, c(0.025,0.975)) # 95% percentiles conflims = apply(yvals,2,quant) # 95% confidence intervals 

Now that I have my confidence limits I create a graph:

meany = sapply(xvals,function (x) SSasymp(x,fixef(fm1)[[1]], fixef(fm1)[[2]], fixef(fm1)[[3]])) par(cex.axis = 2.0, cex.lab=2.0) plot(0, type='n', xlim=c(3,25), ylim=c(0,65), axes=F, xlab="age", ylab="height"); axis(1, at=c(3,1:5 * 5), labels=c(3,1:5 * 5)) axis(2, at=0:6 * 10, labels=0:6 * 10) for(i in 1:14) { data = subset(Loblolly, Loblolly$Seed == unique(Loblolly$Seed)[i]) lines(data$age, data$height, col = "red", lty=3) } lines(xvals,meany, lwd=3) lines(xvals,conflims[1,]) lines(xvals,conflims[2,]) 

Here's the plot: