if( !require('nlme')){
    install.packages('nlme', dependencies=TRUE,  repos='http://cran.rstudio.com/')
}
library(nlme) # needed for gnls

if( !require('sandwich')){
    install.packages('sandwich', dependencies=TRUE,  repos='http://cran.rstudio.com/')
}
library(sandwich)

library(lmtest)

source('../../BookCode/nlrwr/R/SSexp.R')

# Load in the data needed to source the code from the book:
#
load(file='../../BookCode/nlrwr/data/RGRcurve.rda')
load(file='../../BookCode/nlrwr/data/sockeye.rda')
load(file='../../BookCode/nlrwr/data/exp1.rda')
load(file='../../BookCode/nlrwr/data/exp2.rda')

# Needed functions for the transform-both-sides method:
#
source('../../BookCode/nlrwr/R/boxcox.nls.R')

# Pit the transform both sides model to the sockeye data:
#
source('../../BookCode/nlrwr/scripts/chap6.R')

# Lets verfy that after the transformation we have an adequate model fit:
#
spawnersVals = with(sockeye, seq(min(spawners), max(spawners), length.out=100))

recruits_predicted_m1 = predict(sockeye.m1, newdata=data.frame(spawners=spawnersVals))

recruits_predicted_m2 = predict(sockeye.m2, newdata=data.frame(spawners=spawnersVals))
recruits_predicted_m2 = ( -0.2 * recruits_predicted_m2 + 1 )^(-1/0.2) # I have to invert the empirical BCox transforamtion to get the predictions

#postscript("../../WriteUp/Graphics/Chapter6/chap_6_prob_2_data_N fits.eps", onefile=FALSE, horizontal=FALSE)
plot(sockeye$spawners, sockeye$recruits, xlab='spawners',  ylab='recruits')
lines(spawnersVals, recruits_predicted_m1, col='red')
lines(spawnersVals, recruits_predicted_m2, col='blue')
legend('topleft', c('data','Ricker model','tbs box-cox model'),  lty=c(1, 1, 1), col=c('black', 'blue','red'))
grid()
#dev.off()

# Lets verfy that after the transformation we have normal looking residuals:
#
res = sockeye$recruits[-12] - ( -0.2 * fitted.values(sockeye.m2) + 1 )^(-1/0.2) # residuals
sigma = sd(res)

#postscript("../../WriteUp/Graphics/Chapter6/chap_6_prob_2_std_residual_plot.eps", onefile=FALSE, horizontal=FALSE)
standardRes = res/sigma
qqnorm(standardRes, main='')
abline(a=0, b=1)
grid()
#dev.off()

print(shapiro.test(standardRes))