title: Error in model prediction with lme
body: <p>I'm trying to make predictions from my lme-model:</p>
  
  <pre><code>mod&lt;-lme(height~direction+time+distnest+loop+twc+twc:direction,random=~1|bird_id/FT_no,data=dat,correlation=corAR1(0.5))      
  </code></pre>
  
  <p>Unfortunately I have a nested random and 5 fixed effects and do not really know how to handle them all.</p>
  
  <p>All fixed effects are numeric, except "direction" and "time", which are two-level factors ("land" or "sea"; "day" or "night").</p>
  
  <p>For model prediction I've tried the following:</p>
  
  <pre><code>newdat&lt;-data.frame(loop=seq(min(dat$loop),max(dat$loop),length=100),
                 direction=factor("land",levels(dat$direction)),
                 time_code=factor("1",levels(dat$time)),
                 distnest=mean(dat$distnest),
                 twc=mean(dat$twc))
  
  newdat$pred&lt;-predict(mod,newdata=newdat,level=0)
  plot(dat$loop,dat$height,pch=16,las=1,cex.lab=1.2)
  lines(newdat$loop,yhat,lwd=2) ### for plotting one of the fixed effects
  newdat$predse&lt;-predict(mod,newdat,se.fit=TRUE)$se.fit
  </code></pre>
  
  <p>But when using the se.fit=TRUE there is an error: "Error in predict.lme(mod3, newdat, se.fit = TRUE) : cannot evaluate groups for desired levels on 'newdata'"</p>
  
  <p>Does se.fit not work for lmes? Omitting the level=0 did not work. Is the code wrong?</p>
  
  <p>I also tried the code from glmm.wikidot:</p>
  
  <pre><code>newdat&lt;-expand.grid(direction=c("land","sea"),time_code=c("1","2"),loop30=c(0.08,1),distne    st=c(0.01,99.43),twc=c(-56.88744,57.93735))
  newdat$pred&lt;-predict(mod3,newdat,level=0)
  Designmat &lt;- model.matrix(eval(eval(mod3$call$fixed)[-2]), newdat[-ncol(newdat)]) 
  predvar &lt;- diag(Designmat %*% mod3$varFix %*% t(Designmat)) 
  newdat$SE &lt;- sqrt(predvar) 
  newdat$SE2 &lt;- sqrt(predvar+mod3$sigma^2)
  library(ggplot2) 
  pd &lt;- position_dodge(width=0.4) 
  g0 &lt;- ggplot(newdat,aes(x=loop30,y=pred,colour=direction))+ 
  geom_point(position=pd) 
  g0 + geom_linerange(aes(ymin=pred-2*SE,ymax=pred+2*SE), position=pd)
  ## prediction intervals 
  g0 + geom_linerange(aes(ymin=pred-2*SE2,ymax=pred+2*SE2), position=pd)
  </code></pre>
  
  <p>But the plot is completely wrong. Can anyone help me with the correct code? Thanks a lot in advance.</p>
  
  <p>Best wishes, Anna</p>