[R-sig-ME] Estimating variance in linear mixed models
Andrew Miles
rstuff.miles at gmail.com
Wed Feb 1 21:31:05 CET 2012
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
More information about the R-sig-mixed-models
mailing list