Effects Mixed Models for Fish Growth
Sanford Weisberg, [e-mail] October 18, 2008
Much like tree rings on trees, many fish preserve a record of their growth history in annular rings on fish scales and other bony parts. The number of rings provides the age of the fish at capture, and the increment between the rings is a measure of annual growth. We use data from four samples (two years, two gear types) of walleye on West Bearskin Lake, Minnesota. > loc <- "http://www.stat.umn.edu/~sandy/SMBassWBlong.txt" > data <- read.table(url(loc), header = TRUE) > names(data) [1] "species" [7] "lencap" [13] "fage" "lake" "radcap" "fyear" "gear" "age" "yearcap" "inc" "fish" "agecap" "yearclass" "year"
> data$fage <- factor(data$age) > data$fyear <- factor(data$year) > library(xtable) The variables year and age are, respectively, the year of capture and age at capture. The yearclass is the year in which the fish was born, and gear is the gear used for fishing, either T for trapnetting done in the spring or E for electrofishing done in the fall. Each row of the data file has information about one increment, so for example a six year old fish would contribute 6 lines to the data file. There were 430 fish in the study. The tables below give counts of the number of increments. > xtable(xtabs(~yearclass + age, data), digits = 0) 1 177 32 1 18 150 39 3 2 177 32 1 18 131 23 0 3 177 32 1 16 84 0 0 4 177 32 0 7 0 0 0 5 177 24 0 0 0 0 0 6 119 14 0 0 0 0 0 7 48 0 0 0 0 0 0
1983 1984 1985 1986 1987 1988 1989
Weisberg (1993, Can J. Fisheries, 50, 1229) proposed a fixed-effects linear model for growth based on annular increments. We define ycka to be the a-th annular increment for the k-th fish from year-class c. Then the basic fixed-effects model can be described by the equation ycka = a + c+a-1 + gck + ecka 1 (1)
where a is the intrinsic growth or annular increment for a fish in the a-th year of life, is the environmental effect for year = c + a - 1, which is the year in which a fish in year class c will be of age a. The ecka were assumed to be independent errors with mean zero and common variance 2 . Model (1) does not include an intercept, which allows an estimate to be obtained for a for each a. is a parameter to be estimated and gck is a dummy variable equal to one if gear type used to catch the k-th fish in year-class c is T and equal to zero if the gear type is E. This model is easily fit with lm: > m1 <- lm(inc ~ 0 + fage + fyear + gear, data = data) The fixed-effects growth model has a number of defects. First, the estimates of the a and the are not unique. For any constant m, adding m to each of the ^a and simultaneously subtracting m from all the would give a new set of estimates that gives the same fitted ^ values as those produced above by R. R produces unique estimates by setting = 0 for ^ equal to the earliest year in the data; SAS would produce unique estimates by setting to zero the environmental effect with equal to the last year in the data, and other programs might do something else entirely. The result of this is to make the estimates depend on the years of sampling, and so comparison of estimates between locations or over disconnected time periods may be impossible. However, in comparisons between years or between ages within a single study the arbitrary constant m will cancel so these comparisons are meaningful. Second, the assumption that the increments are independently distributed, with common constant variance was originally made for computational convenience, but can hardly be expected to reflect the reality of biological growth. All fish within a year are likely to be correlated, as environmental effects are likely to impact all fish in that year. Similarly, increments on the same fish will in most populations be (positively) correlated. Fish that grow poorly in year would typically be more likely to grow poorly in year + 1, while relatively larger fish in a year-class are more likely to continue to exhibit more growth in succeeding years.
A mixed-effects model
Each fish contributes more than one measurement, and so a random effect between fish should be expected. Similarly, treating years as a random effect, at least as a baseline model, can make sense in problems with data from many years. We now assume that the random environmental effects h are random draws from a normal 2 distribution with mean zero and variance h . Similarly, let fck be the random effect with zero 2 mean and variance f for the k-th fish of year-class c. The mixed model is ycka = a + hc+a-1 + gck + fck + ecka (2)
The R fit is given in Table 1 This appears to be the appropriate baseline model for these data, as it includes fixed intrinsic age effects, predicted random effects for years, and also models within-fish correlation by allowing each fish to have its own overall level of growth that applies to all increments on that fish. The estimated within fish correlation is .0502/(.0502 + .0268 + 2 2 0.0825) = 0.314. The estimate of f is larger than the estimate of h , but smaller than the estimate of 2 . Apparently most of the variation in these data is between fish and the random error not otherwise modeled in (2). Not addressed yet are the "estimates" of the year random effects, which are of primary interest to the fish ecologists. These can be recovered using the rather obscure commands shown in Table 1. We will present more on the estimated year and age effects later. Here is a reasonable graph that summarizes the fish effects, namely a density estimate. The syntax of the ranef command is unfortunately very obscure. 2
> library(lme4) > print(n2 <- lmer(inc ~ 0 + fage + gear + (1 | fyear) + (1 | fish), + data), corr = FALSE) Linear mixed model fit by REML Formula: inc ~ 0 + fage + gear + (1 | fyear) + (1 | fish) Data: data AIC BIC logLik deviance REMLdev 1164 1224 -571 1142 1104 Random effects: Groups Name Variance Std.Dev. fish (Intercept) 0.0502 0.224 fyear (Intercept) 0.0268 0.164 Residual 0.0825 0.287 Number of obs: 1710, groups: fish, 420; fyear, 7 Fixed effects: Estimate Std. Error t value fage1 1.2724 0.0669 19.02 fage2 1.3082 0.0666 19.64 fage3 1.2592 0.0686 18.34 fage4 1.2110 0.0720 16.82 fage5 1.1370 0.0727 15.64 fage6 1.2239 0.0751 16.30 fage7 1.1642 0.0872 13.35 gearT 0.1252 0.0289 4.34 > out <- rbind(t(ranef(n2)$fyear), sqrt(attr(ranef(n2, postVar = TRUE)$fyear, + "postVar")[1, 1, ])) > rownames(out) <- c("Estimate", "SE") > out 1983 1984 1985 1986 1987 1988 1989 Estimate 0.16837 0.11337 -0.26097 -0.18985 0.05770 0.07437 0.03701 SE 0.02512 0.02349 0.02345 0.02266 0.01837 0.01935 0.02498 Table 1: lme4 fit
3
> plot(density(ranef(n2)$fish[, 1]), xlab = "Increment", main = "Fish effects, n2")
Fish effects, n2
2.5 Density 0.0 0.5 1.0 1.5 2.0
-0.5
0.0 Increment
0.5
1.0
1
Mixed-effects with random effect interactions
A valid criticism of the fixed effects model is that if a yearage interaction is present, the model becomes too complex to be of much practical use. A mixed model with interactions is much simpler: ycka = a + hc+a-1 + (ih)a,c+a-1 + ck + fck + ecka (3)
2 This adds only one new variance term ih for interaction rather than a large number of fixed effects terms. The predictions for the year main effects are still valid and useful.
> print(n3 <- lmer(inc ~ 0 + fage + gear + (1 | fyear) + (1 | fyear:fage) + + (1 | fish), data), corr = FALSE) Linear mixed model fit by REML Formula: inc ~ 0 + fage + gear + (1 | fyear) + (1 | fyear:fage) + (1 | Data: data AIC BIC logLik deviance REMLdev 1135 1201 -556 1090 1111 Random effects: Groups Name Variance Std.Dev. fish (Intercept) 0.0498 0.2232 fyear:fage (Intercept) 0.0370 0.1923 fyear (Intercept) 0.0050 0.0707 Residual 0.0795 0.2820 Number of obs: 1710, groups: fish, 420; fyear:fage, 26; fyear, 7 Fixed effects: Estimate Std. Error t value fage1 1.3975 0.0884 15.81 fage2 1.3193 0.0929 14.20 fage3 1.1859 0.1027 11.55 fage4 1.1582 0.1259 9.20 fage5 1.1485 0.1485 7.73 fage6 1.2198 0.1508 8.09 fage7 1.1954 0.2074 5.76 gearT 0.1169 0.0287 4.07 4 fish)
> out1 <- rbind(t(ranef(n3)$fyear), sqrt(attr(ranef(n3, postVar = TRUE)$fyear, + "postVar")[1, 1, ])) > rownames(out1) <- c("Estimate", "SE") > out1 1983 1984 1985 1986 1987 1988 1989 Estimate 0.005973 0.01495 -0.01557 -0.07156 0.01926 0.03656 0.01040 SE 0.066428 0.06310 0.06225 0.05967 0.05705 0.05556 0.05509 Including the interaction decreases AIC by about 29, so the interaction model seems to match the data more closely than the no-interaction data. The year:age interaction is larger than the year variance, suggesting that for this particular data set the interactions are more important than previously claimed. The plot of the year effects shown above indicates that the apparently poor years of 1985 and 1986 may not have been as poor as suggested in the no-interaction model, and the years 1983 and 1984 might be better than previously reported. The standard errors for the year effects are larger in the interaction model as well. We can summarize these three models using a couple of graphs: > coefplot <- function(xaxis, coefs, ses, col = "black", ...) { + require(plotrix) + plotCI(xaxis, coefs, ses, ...) + lines(xaxis, coefs, lty = 2, col = col) + } > > > + > + > + + > + > + > + > + + > + m1.coefs <- summary(m1)$coef par(mfrow = c(1, 2)) coefplot(1:7, fixef(n2)[1:7], sqrt(dia...