R help - Jan 2012 - Linear Mixed Model set-up

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:
> fish<-read.csv("data.csv",header=TRUE)
> fish
   SPECIES YEAR CONTAMINANT
1  Walleye 1970        2.83
2  Walleye 1970        2.56
3  Walleye 1970        2.83
4  Walleye 1970        2.56
5  Walleye 1970        2.77
6  Walleye 1970        2.56
7  Walleye 1970        2.64
8  Walleye 1970        2.22
9  Walleye 1970        2.56
10 Walleye 1970        2.40
11 Walleye 1975        1.59
12 Walleye 1975        1.53
13 Walleye 1975        2.16
14 Walleye 1975        1.60
15 Walleye 1975        2.16
16 Walleye 1976        2.03
17 Walleye 1976        1.97
18 Walleye 1976        1.95
19 Walleye 1976        2.36
20 Walleye 1976        1.82
21 Walleye 1976        1.99
22 Walleye 1977        1.06
23 Walleye 1977        2.00
24 Walleye 1977        1.97
25 Walleye 1977        2.00
26 Walleye 1977        1.99
27 Walleye 1977        1.95
28 Walleye 1977        2.10
29 Walleye 1977        2.29
30 Walleye 1977        2.20
31 Walleye 1979        1.90
32 Walleye 1979        1.98
33 Walleye 1979        2.00
34 Walleye 1979        2.11
35 Walleye 1980        1.92
36 Walleye 1980        2.00
37 Walleye 1980        1.98
38 Walleye 1980        2.25
39 Walleye 1981        1.22
40 Walleye 1981        1.36
41 Walleye 1981        1.48
42 Walleye 1981        1.86
43 Walleye 1981        1.41
44 Walleye 1982        1.25
45 Walleye 1982        1.10
46 Walleye 1982        1.28
47 Walleye 1982        1.28
48 Walleye 1982        1.77
49 Walleye 1982        1.59
50 Walleye 1982        1.61
51 Walleye 1982        1.55
52 Walleye 1984        1.25
53 Walleye 1984        1.41
54 Walleye 1984        1.50
55 Walleye 1984        1.39
> 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.

	[[alternative HTML version deleted]]

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 ...