### R code from vignette source 'fishgrowth.Rnw' ################################################### ### code chunk number 1: fishgrowth.Rnw:10-11 ################################################### options(show.signif.stars=FALSE,digits=4,width=80) ################################################### ### code chunk number 2: one ################################################### loc <- "http://www.stat.umn.edu/~sandy/SMBassWBlong.txt" data <- read.table(loc, header=TRUE)[, c(3, 4, 5, 6, 9, 10, 11, 12)] library(car) some(data, n=5) ################################################### ### code chunk number 3: fishgrowth.Rnw:35-37 ################################################### library(xtable) xtable(xtabs(~ yearclass + age, data), digits=0) ################################################### ### code chunk number 4: fishgrowth.Rnw:62-65 ################################################### data$age <- factor(data$age) data$year <- factor(data$year) m1 <- lm(inc~ 0 + age + year + gear, data=data) ################################################### ### code chunk number 5: fishgrowth.Rnw:120-122 ################################################### library(lme4) print(n2 <- lmer(inc ~ 0 + age + gear + (1|year) + (1|fish),data), corr=FALSE) ################################################### ### code chunk number 6: fishgrowth.Rnw:141-145 ################################################### out<- rbind(t(ranef(n2)$year), sqrt(attr(ranef(n2, postVar=TRUE)$year, "postVar")[1,1, ])) rownames(out) <- c("Estimate", "SE") out ################################################### ### code chunk number 7: fig6 ################################################### plot(density(ranef(n2)$fish[, 1]), xlab="Increment", main="Fish effects, n2") ################################################### ### code chunk number 8: fishgrowth.Rnw:168-171 ################################################### print(n3 <-lmer(inc~0+age+gear+(1|year)+(1|year:age)+ (1|fish),data),corr=FALSE) anova(n2,n3) ################################################### ### code chunk number 9: fishgrowth.Rnw:182-186 ################################################### out1 <- rbind(t(ranef(n3)$year), sqrt(attr(ranef(n3,postVar=TRUE)$year,"postVar")[1,1,])) rownames(out1) <- c("Estimate", "SE") out1 ################################################### ### code chunk number 10: fishgrowth.Rnw:190-194 ################################################### coefplot <- function(xaxis,coefs,ses,col="black",...) { require(plotrix) plotCI(xaxis,coefs,ses,...) lines(xaxis,coefs,lty=2,col=col)} ################################################### ### code chunk number 11: fig7 ################################################### m1.coefs <- summary(m1)$coef par(mfrow=c(1,2)) coefplot(1:7, fixef(n2)[1:7], sqrt(diag(vcov(n2))[1:7]), xlab="Age", pch="N", ylab="Increment", ylim=c(1.0, 1.55)) coefplot(.1+(1:7), fixef(n3)[1:7], sqrt(diag(vcov(n3))[1:7]), pch="I", add=TRUE, col="red", scol="red") coefplot(1:7, m1.coefs[1:7,1], m1.coefs[1:7,2], xlab="Age", pch="F", ylab="Increment", col="blue", scol="blue", add=TRUE) legend("bottomleft", legend=c("NoInt", "Int", "Fixed"), col=c("black", "red", "blue"), pch=c("N","I","F")) # year plot coefplot(1983:1989,out[1, ], out[2, ], pch="N", xlab="Year", ylab="Increment", ylim=c(-.5, .3)) coefplot(.1+(1983:1989),out1[1, ],out1[2, ], add=TRUE, pch="I", col="red", scol="red") coefplot(1983:1989, c(0,m1.coefs[8:13, 1]), c(0,m1.coefs[8:13, 2]), pch="F", xlab="Age", ylab="Increment", col="blue", scol="blue", ylim=c(1.0, 1.55), add=TRUE) legend("bottomleft", legend=c("NoInt", "Int", "Fixed"), col=c("black", "red", "blue"), pch=c("N", "I", "F"))