Hello all,

I'm comparing the duration (duration) of two different behaviour 
(behaviourtype) for 10 individuals (indid) in two different years 
(year). One of the behaviour was not observed for some individuals in 
the first year, so the indid and year are only partially crossed (see 
end of mail for data distribution). I'm fitting the model:

lmer(duration~behaviourtype+(1|indid)+(1|year),REML=T,data=mydata)

I have only two levels for year so modeling it as a random effect may 
not be that important (suitable?).
The model summary is this:

Linear mixed model fit by REML
Formula: duration ~ behaviourtype + (1 | year) + (1 | indid)
    Data: mydata
   AIC  BIC logLik deviance REMLdev
  1708 1724 -848.8     1706    1698
Random effects:
  Groups   Name        Variance Std.Dev.
  indid    (Intercept)  13.543   3.6800
  year     (Intercept)  17.683   4.2051
  Residual             410.779  20.2677
Number of obs: 192, groups: indid, 10; year, 2

Fixed effects:
                       Estimate Std. Error t value
(Intercept)            32.8539     3.9592   8.298
behaviourtypeForaging   0.7121     3.0569   0.233

Correlation of Fixed Effects:
             (Intr)
bhvrtypFrgn -0.449


I was quite happy with this and wanted to perform a LR test with the 
model without the fixed effect (behaviourtype) - I gathered that one 
should rather use the likelihood estimated via ML (and that is what is 
done by the anova() function anyway), so for my 'personal interest' I 
fitted the model again with REML=F:

Linear mixed model fit by maximum likelihood
Formula: duration ~ behaviourtype + (1 | year) + (1 | indid)
    Data: mydata
   AIC  BIC logLik deviance REMLdev
  1716 1733 -853.2     1706    1699
Random effects:
  Groups   Name        Variance   Std.Dev.
  indid    (Intercept) 7.4845e-12 2.7358e-06
  year     (Intercept) 5.5096e+00 2.3473e+00
  Residual             4.2026e+02 2.0500e+01
Number of obs: 192, groups: indid, 10; year, 2

Fixed effects:
                       Estimate Std. Error t value
(Intercept)             32.228      2.814  11.453
behaviourtypeForaging    1.590      3.005   0.529

Correlation of Fixed Effects:
             (Intr)
bhvrtypFrgn -0.604


The variance estimates of the random effects are VERY different from the 
REML estimates.
Admittedly my understanding of the fitting process of lmer is rather 
small, and doing some google work or reading Zuur et al. or Gelman & 
Hill books did not improve my understanding of what is happening here (I 
may have to learn to read better). In particular, my questions for you guys:
1) is my model implementation right ?
2) are differences in ML and REML estimates affecting the LR test of the 
model against a 'null' model (no fixed effect) ?

Many thanks for sharing some of your most valuable resources (brain, time).
simon

PS/ This is the data distribution across behaviour, years, and the 10 
individuals:

, ,  = year1


               1  2  3  4  5  6  7  8  9 10
   behaviour1  2  1  3  2  4  0  3  1  0  7
   behaviour2  4  5  7  4  5 10  0  1  7  5

, ,  = year2


               1  2  3  4  5  6  7  8  9 10
   behaviour1 13  1  5  6  1 13  3  9  2 10
   behaviour2  4  4  8  7  6  3  9  6  8  3



I'm happy to make mine the fortune:

Getting flamed for asking dumb questions on a public mailing list is all part of growing up and being a man/woman.
    -- Michael Watson (in a discussion on whether answers on R-help should be more polite)
       R-help (December 2004)






	[[alternative HTML version deleted]]


