[R] Confidence intervals for non-lm curve
Irene Mantzouni
ima at aqua.dtu.dk
Thu Jun 26 15:47:46 CEST 2008
Hi all!
I would like to estimate confidence intervals for a non lm model.
For example, I use a mixed model of the form:
md=lme(y~x1+I(x1^2)+x2 ...)
Parameters x1+I(x1^2) are fixed effects and I would like to plot the predicted (partial) curve corresponding to these ones, along with 90% CI bands.
Thus, I simulate (partial) predictions:
Lx=c()
Ux=c()
Mx=c()
z=c()
for (j in 1:length(x1)){
x=x1[j]
for(i in 1:2000){
s1 <- rnorm(1,.026,.027) # mean and sd estimated by lme for x1
s2 <- rnorm(1,-.01,.005) # mean and sd estimated by lme for I(x1^2)
z[i] <- s1*x+s2*(x^2)
}
Lx[j]=quantile(z,.05)
Ux[j]=quantile(z,.95)
Mx[j]=mean(z)
}
And then plot vectors Lx, Mx and Ux for lower, mean and upper curves, respectively.
Is this approach correct?
Any alternatives?
Thank you!!
More information about the R-help
mailing list