Estimating variance in linear mixed models
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