#
# Epage > 137
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

if(is.null(version$language) == FALSE){ 
  require(alr3)
}else{
  library(alr3)
}

if(is.null(version$language) == FALSE) data(walleye)
attach(walleye)

postscript("../../WriteUp/Graphics/Chapter11/prob_3_scatterplot.eps", onefile=FALSE, horizontal=FALSE)
plot( age, length, pch=period ) 
dev.off() 

# Note: we compare this data (by period) in the same way in which the "three sources of methionine"
# example in the book was done 
#

# walleye: with three periods
tdata <- walleye 
# create the indicators for the categories of P
tdata$P1 <- tdata$P2 <- tdata$P3 <- rep(0,dim(tdata)[1])
tdata$P1[tdata$period==1] <- 1
tdata$P2[tdata$period==2] <- 1
tdata$P3[tdata$period==3] <- 1

# use the same procedure as in prob_2.R to get estimates for the initial guess at the parameters
# 
LInfinity0     <- 1.05*max(length) # this is an initial guess for L_{\infty}
minusKTMinusT0 <- log( 1 - length / LInfinity0 ) # this is the response of the linear fit 

mf <- lm( minusKTMinusT0 ~ age ) # fit a linear model and extract \beta_{0} and \beta_1
b0 <- coefficients(mf)[1]
b1 <- coefficients(mf)[2]

# here are our initial guesses at the parameters K and t_0: 
K0  <- -b1
t00 <- -b0/b1 

# fit the various possible models: 
# 
# common regressions
m0 <- nls( length ~ L * ( 1 - exp(-K*(age-t0)) ),
           data=tdata,start=list(L=LInfinity0,K=K0,t0=t00))
# most general
m1 <- nls( length ~ ( P1*( L1 * ( 1 - exp(-K1*(age-t01)) ) ) +
                      P2*( L2 * ( 1 - exp(-K2*(age-t02)) ) ) +
                      P3*( L3 * ( 1 - exp(-K3*(age-t03)) ) ) ), 
           data=tdata,start= list(L1=LInfinity0,K1=K0,t01=t00,
                                  L2=LInfinity0,K2=K0,t02=t00,
                                  L3=LInfinity0,K3=K0,t03=t00) )
# common intercept
m2 <- nls(length ~ ( P1*( L * ( 1 - exp(-K1*(age-t01)) ) ) +
                     P2*( L * ( 1 - exp(-K2*(age-t02)) ) ) +
                     P3*( L * ( 1 - exp(-K3*(age-t03)) ) ) ), 
           data=tdata,start= list(L=LInfinity0,K1=K0,t01=t00,K2=K0,t02=t00,K3=K0,t03=t00) )

# common intercept-rate model
m3 <- nls( length ~ ( P1*( L * ( 1 - exp(-K*(age-t01)) ) ) +
                      P2*( L * ( 1 - exp(-K*(age-t02)) ) ) +
                      P3*( L * ( 1 - exp(-K*(age-t03)) ) ) ), 
           data=tdata,start= list(L=LInfinity0,K=K0,t01=t00,t02=t00,t03=t00) )

# common intercept-origin model
m4 <- nls( length ~ ( P1*( L * ( 1 - exp(-K1*(age-t0)) ) ) +
                      P2*( L * ( 1 - exp(-K2*(age-t0)) ) ) +
                      P3*( L * ( 1 - exp(-K3*(age-t0)) ) ) ), 
           data=tdata,start= list(L=LInfinity0,K1=K0,K2=K0,K3=K0,t0=t00) )

anova(m0,m1) # the more specific model is better 

attach(tdata)

postscript("../../WriteUp/Graphics/Chapter11/prob_3_scatterplot.eps", onefile=FALSE, horizontal=FALSE)

plot(age,length,pch=period)
xx <- seq(0,max(age),length=100)

# plot predictions for the models above:
# 
#lines( xx, predict(m0,data.frame(age=xx)),lty=1) # common regression 
lines( xx, predict(m1,data.frame(age=xx,P1=rep(1,100), P2=rep(0,100),P3=rep(0,100))),lty=2) # most general Period==1 
lines( xx, predict(m1,data.frame(age=xx,P1=rep(0,100), P2=rep(1,100),P3=rep(0,100))),lty=3) # most general Period==2 
lines( xx, predict(m1,data.frame(age=xx,P1=rep(0,100), P2=rep(0,100),P3=rep(1,100))),lty=4) # most general Period==3 

dev.off()