title: Extract prediction band from lme fit
body: <p>I have following model </p>
  
  <pre><code>x  &lt;- rep(seq(0, 100, by=1), 10)
  y  &lt;- 15 + 2*rnorm(1010, 10, 4)*x + rnorm(1010, 20, 100)
  id &lt;- NULL
  for(i in 1:10){ id &lt;- c(id, rep(i,101)) }
  dtfr &lt;- data.frame(x=x,y=y, id=id)
  library(nlme)
  with(dtfr, summary(     lme(y~x, random=~1+x|id, na.action=na.omit)))
  model.mx &lt;- with(dtfr, (lme(y~x, random=~1+x|id, na.action=na.omit)))
  pd       &lt;- predict( model.mx, newdata=data.frame(x=0:100), level=0)
  with(dtfr, plot(x, y))
  lines(0:100, predict(model.mx, newdata=data.frame(x=0:100), level=0), col="darkred", lwd=7)
  </code></pre>
  
  <p><img src="https://i.stack.imgur.com/atxBI.png" alt="enter image description here"></p>
  
  <p>with <code>predict</code> and <code>level=0</code> i can plot the mean population response. How can I extract and plot the 95% confidence intervals / prediction bands from the nlme object for the whole population?</p>