Hi everyone,

I have a problem with the model prediction of my lme model. The model is:

mod<-lme(height~direction+time+distnest+loop+twc+twc:direction,random=~1|bird_id/FT_no,data=dat,correlation=corAR1(0.5))

The summary looks as follows:

Linear mixed-effects model fit by REML
  Data: b
        AIC      BIC    logLik
   876.3957 919.9384 -427.1979

Random effects:
  Formula: ~1 | bird_id
         (Intercept)
StdDev:   0.1366816

  Formula: ~1 | FT_no %in% bird_id
          (Intercept)  Residual
StdDev: 0.0001120659 0.7654412

Correlation Structure: AR(1)
  Formula: ~1 | bird_id/FT_no
  Parameter estimate(s):
       Phi
0.4649851
Fixed effects: log(height + 11) ~ direction + time_code + distnest + loop30 +      twc + twc:direction
                       Value  Std.Error  DF   t-value p-value
(Intercept)       2.5815572 0.20754492 360 12.438547  0.0000
directionsea     -0.3241188 0.10184748 360 -3.182394  0.0016
time_code2       -0.3646015 0.14809643 360 -2.461919  0.0143
distnest          0.0159491 0.00288208 360  5.533870  0.0000
loop30            0.5399952 0.19264292 360  2.803089  0.0053
twc              -0.0004175 0.00269104 360 -0.155147  0.8768
directionsea:twc  0.0060729 0.00393018 360  1.545194  0.1232
  Correlation:
                  (Intr) drctns tm_cd2 dstnst loop30 twc
directionsea     -0.334
time_code2       -0.082 -0.185
distnest         -0.456  0.237  0.055
loop30           -0.841  0.089 -0.009  0.176
twc               0.066 -0.004  0.003  0.005 -0.077
directionsea:twc -0.027 -0.014 -0.004 -0.012  0.031 -0.646

Standardized Within-Group Residuals:
         Min          Q1         Med          Q3         Max
-3.18912725 -0.56260570 -0.08412658  0.38321599  4.71633567

Number of Observations: 394
Number of Groups:
            bird_id FT_no %in% bird_id
                  8                 28


A short excerpt of my data:

>head(dat)
         lat      lon height method  project   bird_id       date day month
2  53.67150 7.098743      2 Logger WINDBIRD HJ14_2013 21.05.2013  21     5
3  53.64794 7.145718     12 Logger WINDBIRD HJ14_2013 21.05.2013  21     5
4  53.62625 7.206992     56 Logger WINDBIRD HJ14_2013 21.05.2013  21     5
8  53.59216 7.235380     35 Logger WINDBIRD HJ14_2013 21.05.2013  21     5
9  53.56353 7.245068     39 Logger WINDBIRD HJ14_2013 21.05.2013  21     5
10 53.53348 7.266665     25 Logger WINDBIRD HJ14_2013 21.05.2013  21     5
    year hour minute second     time speed FT_no distflown speed_poin distnest
2  2013    7     13      1 07:13:01  36.0     1      1.81      19.06     1.81
3  2013    7     18     43 07:18:43  44.9     1      4.06      42.77     5.87
4  2013    7     24     20 07:24:20  44.9     1      4.72      50.39    10.53
8  2013    7     46     56 07:46:56  34.0     1      2.61      27.75    14.44
9  2013    7     52     35 07:52:35  48.7     1      3.25      34.50    17.23
10 2013    7     58     13 07:58:13  42.9     1      3.64      38.73    20.74
    timediff trcourse direction trdur time_sec  sex course winddir
2       339     0.84       sea   339    25981 <NA>    141     230
3       342     2.73       sea   681    26323 <NA>    128     230
4       337     4.93      land  1018    26660 <NA>    140     230
8       338     6.97      land  2374    28016 <NA>    205     230
9       339     8.48      land  2713    28355 <NA>    115     230
10      338    10.17      land  3051    28693 <NA>    150     230
    windspeed_ms windspeed winddiff seg_15    loop15 seg_30 loop30 colony prec
2           4.7     16.92       89      1 1.0000000      1   0.99  Juist    0
3           4.7     16.92      102      1 1.0000000      1   0.99  Juist    0
4           4.7     16.92       91      2 0.9753846      1   0.99  Juist    0
8           4.7     16.92       25      3 1.0000000      2   0.99  Juist    0
9           4.7     16.92      115      4 0.9934518      2   0.99  Juist    0
10          4.7     16.92       81      4 0.9934518      2   0.99  Juist    0
    cloudcover prec_f                 cet time_code winddiff2       twc
2           8      0 2013-05-21 07:13:01         1       -89  8.632196
3           8      0 2013-05-21 07:18:43         1      -102  1.718830
4           8      0 2013-05-21 07:24:20         1       -90 -7.581406
8           8      0 2013-05-21 07:46:56         1       -25 16.771152
9           8      0 2013-05-21 07:52:35         1      -115 -5.512702
10          8      0 2013-05-21 07:58:13         1       -80 -1.867752




For plotting the output for each relevant fixed effect I first of all used "predict",

newdat<-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<-predict(mod,newdata=newdat,levels=0)
plot(dat$loop,dat$height,pch=16,las=1,cex.lab=1.2)
lines(newdat$loop,yhat,lwd=2)
newdat$predse<-predict(mod,newdat,se.fit=TRUE)$se.fit

but got this error message:Error in predict.lme(mod3, newdata = newdat, levels = 0) :
   cannot evaluate groups for desired levels on 'newdata'


Did I make a mistake in my predict code?

Alternatively I followed the help from glmm.wikidot.com/faq:

newdat<-expand.grid(direction=c("land","sea"),time=c("1","2"),loop=c(0.08,1),distnest=c(0.01,99.43),twc=c(-56.88744,57.93735))
newdat$pred<-predict(mod,newdat,level=0)
Designmat <- model.matrix(eval(eval(mod$call$fixed)[-2]), 
newdat[-ncol(newdat)])
predvar <- diag(Designmat %*% mod$varFix %*% t(Designmat))
newdat$SE <- sqrt(predvar)
newdat$SE2 <- sqrt(predvar+mod$sigma^2)

library(ggplot2)
pd <- position_dodge(width=0.4)
g0 <- ggplot(newdat,aes(x=loop,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)


But the plot seems more than wrong. Can anyone help me with my problem?

Thanks a lot in advance.

Best wishes,
Anna


	[[alternative HTML version deleted]]

