Extract prediction band from lme fit

by

Warning: Read this thread on r-sig-mixed models before doing this. Be very careful when you interpret the resulting prediction band.

From r-sig-mixed models FAQ adjusted to your example:

set.seed(42)
x <- rep(0:100,10)
y <- 15 + 2*rnorm(1010,10,4)*x + rnorm(1010,20,100)
id<-rep(1:10,each=101)

dtfr <- data.frame(x=x ,y=y, id=id)

library(nlme)

model.mx <- lme(y~x,random=~1+x|id,data=dtfr)

#create data.frame with new values for predictors
#more than one predictor is possible
new.dat <- data.frame(x=0:100)
#predict response
new.dat$pred <- predict(model.mx, newdata=new.dat,level=0)

#create design matrix
Designmat <- model.matrix(eval(eval(model.mx$call$fixed)[-2]), new.dat[-ncol(new.dat)])

#compute standard error for predictions
predvar <- diag(Designmat %*% model.mx$varFix %*% t(Designmat))
new.dat$SE <- sqrt(predvar) 
new.dat$SE2 <- sqrt(predvar+model.mx$sigma^2)

library(ggplot2) 
p1 <- ggplot(new.dat,aes(x=x,y=pred)) + 
geom_line() +
geom_ribbon(aes(ymin=pred-2*SE2,ymax=pred+2*SE2),alpha=0.2,fill="red") +
geom_ribbon(aes(ymin=pred-2*SE,ymax=pred+2*SE),alpha=0.2,fill="blue") +
geom_point(data=dtfr,aes(x=x,y=y)) +
scale_y_continuous("y")
p1

plot

More Related Contents:

  1. Getting Warning: ” ‘newdata’ had 1 row but variables found have 32 rows” on predict.lm
  2. How does predict.lm() compute confidence interval and prediction interval?
  3. How to get coefficients and their confidence intervals in mixed effects models?
  4. lme4::lmer reports “fixed-effect model matrix is rank deficient”, do I need a fix and how to?
  5. Working with neuralnet in R for the first time: get “requires numeric/complex matrix/vector arguments”
  6. Why does the output of of a linear mixed model using lme4 show one level of a factor but not another?
  7. how to debug “factor has new levels” error for linear model and prediction [duplicate]
  8. Is there a dplyr equivalent to data.table::rleid?
  9. Capitalize the first letter of both words in a two word string
  10. `levels
  11. Creating arbitrary panes in ggplot2
  12. Getting path of an R script
  13. Difference: “Compile PDF” button in RStudio vs. knit() and knit2pdf()
  14. Rotating x axis labels in R for barplot
  15. How do I convert a factor into date format?
  16. Pass a vector of variable names to arrange() in dplyr
  17. How subset a data frame by a factor and repeat a plot for each subset?
  18. change background and text of strips associated to multiple panels in R / lattice
  19. Check whether values in one data frame column exist in a second data frame
  20. adding default values to item x group pairs that don’t have a value (df %>% spread %>% gather seems strange)
  21. Milliseconds in POSIXct Class
  22. heatmap with values (ggplot2)
  23. How to choose variable to display in tooltip when using ggplotly?
  24. Using different scales as fill based on factor
  25. Conditional binary join and update by reference using the data.table package
  26. predict.lm() with an unknown factor level in test data
  27. How to read csv file in R where some values contain the percent symbol (%)
  28. Is there a way to output text to the R console in color
  29. Split a string every 5 characters
  30. Using a static (prebuilt) PDF vignette in R package

Leave a Comment