###################################################
### chunk number 1: 
###################################################
options(show.signif.stars=FALSE,digits=4)


###################################################
### chunk number 2: 
###################################################
loc <- "http://www.stat.umn.edu/~sandy/SMBassWB.txt"
wdata <- read.table( url(loc),header=TRUE)
wdata <- wdata[wdata$yearcap-wdata$agecap>1982,]
growthFrame <- function(data,idvars=1:8,incvars=9:22,increments=TRUE,
    yearcap="yearcap",agecap="agecap",fishid=NULL) {
 require(reshape)
 if(!increments) { # translate radii to increments
   data[,incvars] <- t(apply(cbind(0,data[,incvars]),1,diff))}
 data$FishID <- if(is.null(fishid)) 1:dim(data)[1] else data[[fishid]]
 data$YearClass <- data[[yearcap]]-data[[agecap]]
 p <- dim(data)[2]
 ids <- c(idvars,p-1,p)
 data1 <- melt(data,id=ids,measure.vars=incvars,variable_name="Age")
 names(data1)[length(ids)+2] <- c("Increment")
 data1 <- data1[!is.na(data1$Increment),]
 data1$Age <- as.numeric(data1$Age)
 data1$Year <- factor(data1[[yearcap]]-data1[[agecap]]+data1$Age-1)
 data1$Age  <- factor(data1$Age)
 data1 <- data1[order(data1$FishID),]
 data1$FishID <- factor(data1$FishID)
 data1}
wdata1 <- growthFrame(wdata,idvars=c(1:7,20),incvars=8:19,increments=FALSE)


###################################################
### chunk number 3: 
###################################################
m1 <- lm(Increment ~ 0+ Age + Year + gear, data=wdata1)
summary(m1)


###################################################
### chunk number 4: 
###################################################
m2 <- update(m1, ~.+ Age:Year)
anova(m1,m2)
AIC(m1,m2)


###################################################
### chunk number 5: 
###################################################
loc <- "http://www.stat.umn.edu/~sandy/BeautifulsAll.txt"
hdata <- read.table(url(loc),header=TRUE)
hdata1 <- growthFrame(hdata)


###################################################
### chunk number 6: 
###################################################
h1 <- lm(Increment ~ 0+ Age + Year, hdata1)
h2 <- update(h1, ~.+Age:Year)
anova(h1,h2)
AIC(h1,h2)


###################################################
### chunk number 7: 
###################################################
library(lme4)
n1 <- lmer(Increment~0+Age+gear+(1|Year)+(1|FishID),wdata1)
print(n1,corr=FALSE)


###################################################
### chunk number 8: 
###################################################
n1a <- lmer(Increment~0+Age+gear+(1|Year),wdata1)
anova(n1a,n1,thest="Chisq")


###################################################
### chunk number 9: 
###################################################
hn1 <- lmer(Increment ~ 0+Age + (1|Year) + (1|FishID), hdata1)
print(hn1,corr=FALSE)


###################################################
### chunk number 10: 
###################################################
print(hn2 <- update(hn1, ~.+as.numeric(Year)),corr=F)


###################################################
### chunk number 11: 
###################################################
hn3 <- lmer(Increment ~ 0+Age + (1|Year) + (1|FishID)+(1|Year:Age),hdata1)
anova(hn1,hn3)


###################################################
### chunk number 12: 
###################################################
loc <- "http://www.stat.umn.edu/~sandy/SMBassWB.txt"
wdata <- read.table( url(loc),header=TRUE)


###################################################
### chunk number 13:  eval=FALSE
###################################################
## wdata <- read.table("SMBassWB.txt",header=TRUE)


###################################################
### chunk number 14: 
###################################################
wdata <- wdata[wdata$yearcap-wdata$agecap>1982,]
wdata[c(1,5,50,100,150),]


###################################################
### chunk number 15: 
###################################################
wdata1 <- growthFrame(wdata,idvars=c(1:7,20),incvars=8:19,increments=FALSE)
wdata1[c(1,5,50,100,150),]


###################################################
### chunk number 16: 
###################################################
fixef(n1)


###################################################
### chunk number 17: 
###################################################
age.n1 <- fixef(n1)[1:7]


###################################################
### chunk number 18: 
###################################################
(age.se.n1 <- sqrt(diag(vcov(n1))[1:7]))


###################################################
### chunk number 19: 
###################################################
(year.n1 <- ranef(n1)$Year[,1])


###################################################
### chunk number 20: 
###################################################
fish.n1 <- ranef(n1)$FishID[,1]


###################################################
### chunk number 21: figden
###################################################
plot(density(fish.n1),main="Fish effects", xlab="Increment (mm)",las=1)


###################################################
### chunk number 22: 
###################################################
var.year.n1 <- attr(ranef(n1,postVar=TRUE)[["Year"]],"postVar")


###################################################
### chunk number 23: 
###################################################
(year.se.n1 <- sqrt(attr(ranef(n1,postVar=TRUE)[["Year"]],"postVar")[1,1,]))


###################################################
### chunk number 24: 
###################################################
growthCoefPlot <- function(xaxis,coefs,ses,col.line="black",...) {
require(plotrix)
plotCI(xaxis,coefs,ses,...)
lines(xaxis,coefs,lty=2,col=col.line)}


###################################################
### chunk number 25: fig4
###################################################
par(mfrow=c(1,2))
growthCoefPlot(1:7,age.n1,age.se.n1,xlab="Age (years)",
       ylab="Increments (mm)",las=1)
growthCoefPlot(as.numeric(levels(wdata1$Year)),year.n1,year.se.n1,
       xlab="Year",ylab="Increments (mm)",las=1)


###################################################
### chunk number 26: fig3
###################################################
n3 <- lmer(Increment~0+Age+gear+(1|Year)+ (1|Age:Year)+(1|FishID),wdata1)
year.se.n3 <- sqrt(attr(ranef(n3,postVar=TRUE)[["Year"]],"postVar")[1,1,])
par(mfrow=c(1,2))
growthCoefPlot(1:7,fixef(n1)[1:7],sqrt(diag(vcov(n1))[1:7]),xlab="Age (years)",
       ylab="Increment (mm)",ylim=c(1.0,1.55),las=1)
growthCoefPlot(.1+(1:7),fixef(n3)[1:7],sqrt(diag(vcov(n3))[1:7]),
       add=TRUE,col="black",scol="black",pch="x")
growthCoefPlot(1983:1989,ranef(n1)$Year[,1],year.se.n1,
      xlab="Year",ylab="Increment (mm)",las=1)
growthCoefPlot(.1+(1983:1989),ranef(n3)$Year[,1],year.se.n3,add=TRUE,
      col="black",scol="black",pch="x")


###################################################
### chunk number 27: 
###################################################
set.seed(12345)


###################################################
### chunk number 28: 
###################################################
(temps <- data.frame(Year=levels(wdata1$Year), Temp=62+3*rnorm(7)))


###################################################
### chunk number 29: 
###################################################
meanTemp <- apply(wdata1,1,
     function(x) temps[temps[,1]==as.numeric(x[["Year"]]),2])


###################################################
### chunk number 30: 
###################################################
n4 <- update(n1,~.+meanTemp)
print(n4,corr=FALSE)


###################################################
### chunk number 31: 
###################################################
am <- 7
ages <- 1:am
ya <- cumsum(fixef(n1)[ages])
V <- vcov(n1)[1:7,1:7]
C <- matrix(1,nrow=am,ncol=am)
for(i in 1:(am-1)){ for(j in (i+1):am) {C[i,j] <- 0}}
Vy <- C %*% V %*% t(C)
Vy <- .5 *(Vy + t(Vy))


###################################################
### chunk number 32: 
###################################################
vbf <- function(param,ages) {
   Linf <- param[1]
   K0 <- param[2]
   a0 <- param[3]
   Linf*(1-exp(-(log(2)/K0)*(ages-a0)))}
negllik <- function(param,R,ages) {
 fit <- vbf(param[1:3],ages)
 tau <- param[4]
 (am/2)*log(tau) + (1/(2*tau*tau))*
             sum( backsolve(R,ya-fit,transpose=T)^2 ) }


###################################################
### chunk number 33: 
###################################################
param <- 1.05*max(ya)
param[2] <- 0.5 * param
param[3] <- 0
param[4] <- 1
names(param) <- c("Linf", "K0", "a0", "tau")
param


###################################################
### chunk number 34: 
###################################################
ans <- optim(param,negllik,R=chol(Vy),age=ages,
   upper=c(Inf,Inf,0,10),lower=c(0,0,-Inf,0),method="L-BFGS-B",hessian=TRUE)
ans$par
sqrt(diag(solve(ans$hessian[-3,-3])))


###################################################
### chunk number 35: 
###################################################
ya <- cumsum(fixef(hn1))
am <- length(ya)
ages <- 1:am
V = vcov(hn1)[ages,ages]
C = matrix(1,nrow=am,ncol=am)
for (i in 1:(am-1)){for (j in (i+1):am) {C[i,j]=0}}
Vy = C %*% V %*% t(C)
Vy <- (Vy + t(Vy))/2
#Starting values
ans<-optim(param,negllik,R=chol(Vy),age=ages,
    upper=c(Inf,Inf,0,10),lower=c(0,0,-Inf,.1),method="L-BFGS-B",hessian=TRUE )
ans$par
sqrt(diag(solve(ans$hessian[-3,-3])))