[R-sig-ME] Estimating variance in linear mixed models

[Andrew Miles]

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]]
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models