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

[Ista Zahn]

Hi Andrew and Maggie,

On Wednesday, February 01, 2012 03:31:05 PM Andrew Miles wrote:
> 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)

The residual variance from the lmer model is the within-year variance, no?

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