[R] Linear Mixed Model set-up

[Ben Bolker]

Maggie Neff <maggie.neff <at> gmail.com> writes:

> 
> Hello,
> 
> I have some data covering contaminant concentrations in fish over a time
> period of ~35 years.  Each year, multiple samples of fish were taken (with
> varying sample sizes each year).  Ultimately, I want an estimation of the
> variance between years, and the variance within years + random effects.  I
> used a linear mixed model to estimate these variances, but after reading a
> number of different references and examples, I am still unclear as to
> whether I have set up the model correctly to obtain these values.
> 
> I've used the *lme* function as follows - the example here is on an
> abbreviated version of my data set:
> 
  [snip]

> > contaminant<-fish$CONTAMINANT
> > year<-fish$YEAR
> > mod<-lme(contaminant~year,random=~1|year,data=data)
> > varcomp(mod,cum=FALSE)
>       year     Within
> 0.02695566 0.05758531
> attr(,"class")
> [1] "varcomp"
> 
> Thanks in advance for your help - I very new to formula-building in R.
> 

  This looks reasonable.  Note that your across-year variation
is not the raw year-to-year variation, but the variation around
the linear trend.

library(nlme)
mod <- lme(contaminant~year,random=~1|year,data=fish)
fishagg <- aggregate(contaminant~year,data=fish,FUN=mean)
pframe <- data.frame(year=1970:1984)

par(bty="l",las=1)
with(fish,plot(contaminant~year))
with(fishagg,points(contaminant~year,pch=3,col=2))
with(data.frame(pframe,
           contaminant=predict(mod,newdata=pframe,level=0)),
     lines(year,contaminant,col=4))
with(data.frame(pframe,
           contaminant=predict(mod,newdata=pframe,level=1)),
     lines(year,contaminant,col=4,lty=2,type="b"))

## don't know why predict() with newdata doesn't work
## on this model ...
with(data.frame(year=fish$year,
           contaminant=predict(mod2,level=0)),
     lines(year,contaminant,col=5))
with(data.frame(year=fish$year,
           contaminant=predict(mod2,level=1)),
     lines(year,contaminant,col=5,lty=2,type="b"))

library(ape)
varcomp(mod)
varcomp(mod2)


The r-sig-mixed-models mailing list is probably best for
mixed model questions ...