load(file='../../BookCode/nlrwr/data/O.mykiss.rda')

mask = complete.cases(O.mykiss)
O.mykiss = O.mykiss[mask, ]

# Plot the data to get an idea of what it looks like:
#
plot(O.mykiss$conc, O.mykiss$weight, xlab='conc', ylab='weight')
grid()

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

fit1 = nls(weight ~ SSexp(conc, b, y0), data=O.mykiss)
print(summary(fit1))

# Now that we have a model, lets study the fitting assumptions following the steps discussed in this chapter of the book:
#
concGrid = with(O.mykiss, seq(min(conc), max(conc), length.out=100))
weightPred_m1 = predict(fit1, newdata=data.frame(conc=concGrid))

#postscript("../../WriteUp/Graphics/Chapter6/chap_6_prob_4_data_N_fits.eps", onefile=FALSE, horizontal=FALSE)
plot(O.mykiss$conc, O.mykiss$weight, xlab='conc', ylab='weight', main='data and fit')
lines(concGrid, weightPred_m1, type='l')
grid()
#dev.off()

# Residual plot:
#
plot(fitted(fit1), residuals(fit1), xlab='fitted values',  ylab='residuals')
abline(h=0, col='black')
grid()

# F-test:
#
n = dim(O.mykiss)[1]
nc = length(unique(O.mykiss$conc))
if( nc <n ){ # we must have repeated measurements
    print('F-test')
    full_model = lm(weight ~ as.factor(conc), data=O.mykiss)
    print(anova(fit1, full_model))
}

# Variance homogeneity:
#
plot(fitted(fit1), abs(residuals(fit1)), xlab='fitted values',  ylab='abs residuals')
grid()

# Normality of the residuals:
#
standardRes = residuals(fit1)/summary(fit1)$sigma
qqnorm(standardRes, main='')
abline(a=0, b=1)
grid()

print(shapiro.test(standardRes))

# Independence:
#
#plot(residuals(fit1), c(residuals(fit1)[-1], NA),  xlab='residuals', ylab='lagged residuals')
#grid()
acf(residuals(fit1))

# Try to transform-both-sides using a logarithmic transformation:
#
fit2 = lm(log(weight) ~ conc, data=O.mykiss)

# Map back to:
#
# weight = y0 exp(conc/b)
#
cs = as.double(coefficients(fit2))
y0 = exp(cs[1])
b = 1/cs[2]

TBS_model = function(conc, b, y0){
    weight = y0 * exp(conc/b)
    weight
}

weightPred_m2 = TBS_model(concGrid, b, y0)

#postscript("../../WriteUp/Graphics/Chapter6/chap_6_prob_4_data_N_fits.eps", onefile=FALSE, horizontal=FALSE)
plot(O.mykiss$conc, O.mykiss$weight, xlab='conc', ylab='weight', main='data and fit')
lines(concGrid, weightPred_m1, type='l', col='blue')
lines(concGrid, weightPred_m2, type='l', col='green')
grid()
#dev.off()