Estimating variance in linear mixed models

Hello,

This is both an R and stats question.  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). The trend is modeled using a linear regression on ln[contaminant]
and year.  The ultimate goal of this project is to do a power analysis for
a monitoring program, and to do this, I need an estimation of both random
within-year variation and random between-year variation.

I used a linear mixed model to estimate these variances.  My questions is
whether the function as I've set it up will return the estimates that I'm
looking for, and if so, which values within the output reflect those
estimates.  If someone knows of a better/alternative way to estimate these
values, that would also be useful.

fish<-read.csv("data.csv",header=TRUE)
fish

  SPECIES YEAR CONT
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

year<-fish$YEAR
cont<-fish$CONT
reg3<-lme(year~cont,random=~1|year,method="REML")
reg3

Linear mixed-effects model fit by REML
 Data: NULL
 Log-restricted-likelihood: 1243.336
 Fixed: year ~ cont
(Intercept)         cont
1.978213e+03 4.947785e-14

Random effects:
Formula: ~1 | year
       (Intercept)    Residual
StdDev:     4.23609 1.18049e-13

Number of Observations: 55
Number of Groups: 9

	[[alternative HTML version deleted]]

A linear multilevel model like you are using can be used to estimate the
between-year variance, like so (note I'm using the lmer function, not lme,
which I am less familiar with):

library(lme4)
mod=lmer(cont ~ 1 + (1|year), data=fish)

The Random Effects for (Intercept) will give you the variance of the
estimates of contaminants in each year.  However, with only 9 different
years (and hence 9 different estimates of contaminants), the random effect
might not be very reliable.  You can get model predicted estimates for each
of the years like so: coef(mod)

I can't immediately think of a way to estimate within-year variance using a
MLM, but you could get estimates of the within-year variance doing
something like this: tapply(cont, year, var)

This should break down cont by year, and then estimate the variance in each
year.

I hope this helps.

Andrew Miles

On Feb 1, 2012, at 1:54 PM, Maggie Neff wrote:

Hello,

This is both an R and stats question.  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). The trend is modeled using a linear regression on ln[contaminant]
and year.  The ultimate goal of this project is to do a power analysis for
a monitoring program, and to do this, I need an estimation of both random
within-year variation and random between-year variation.

I used a linear mixed model to estimate these variances.  My questions is
whether the function as I've set it up will return the estimates that I'm
looking for, and if so, which values within the output reflect those
estimates.  If someone knows of a better/alternative way to estimate these
values, that would also be useful.

fish<-read.csv("data.csv",header=TRUE)
fish

  SPECIES YEAR CONT

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

year<-fish$YEAR
cont<-fish$CONT
reg3<-lme(year~cont,random=~1|year,method="REML")
reg3

Linear mixed-effects model fit by REML

 Data: NULL
 Log-restricted-likelihood: 1243.336
 Fixed: year ~ cont

(Intercept)         cont
1.978213e+03 4.947785e-14

Random effects:
Formula: ~1 | year

       (Intercept)    Residual

StdDev:     4.23609 1.18049e-13

Number of Observations: 55
Number of Groups: 9

	[[alternative HTML version deleted]]