if( !require('investr') ){
    install.packages('investr', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(investr) # has the function plotFit

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

if( !require('nls2')){
    install.packages('nls2', dependencies=TRUE,  repos='http://cran.rstudio.com/')
}
library(nls2) # has the fitting function nls2

plot(lettuce$conc, lettuce$weight, xlab='concentration', ylab='weight')
grid()


brainCousensModel = function(predictor, b, d, e, f){
    ( d + f * predictor )/( 1 + (predictor / e)^b )
}

# Try to get a good starting set of parameters for this model using a grid search:
#
grid.brainCousens = expand.grid(
    list(
        d=seq(1.25, 1.75, length.out=20),
        f=seq(-0.1, 0.1, length.out=20),
        b=seq(-0.1, 0.5, length.out=20),
        e=seq(1, 100, length.out=20)))

lettuce.nls2gsm = nls2(weight ~ brainCousensModel(conc, b, d, e, f), data=lettuce, start=grid.brainCousens, algorithm='brute-force') # gsm = g(rid)s(earch)m(odel)
print('grid search parameters:')
print(coefficients(lettuce.nls2gsm))

lettuce.nls2 = nls2(weight ~ brainCousensModel(conc, b, d, e, f), data=lettuce, start=lettuce.nls2gsm, control=nls.control(maxiter= 100)) # refine the above grid searched initial guess
print('refined parameters:')
print(coefficients(lettuce.nls2))

# Now that we have a model, lets study the fitting assumptions following the steps discussed in this chapter of the book:
#
#postscript("../../WriteUp/Graphics/Chapter5/chap_5_prob_4.eps", onefile=FALSE, horizontal=FALSE)
par(mfrow=c(2, 1))
plotFit(lettuce.nls2gsm, main='lettuce data and fit(initial parameter estimates)')
grid()

plotFit(lettuce.nls2, main='lettuce data and fit(refined parameter estimates)')
grid()
par(mfrow=c(1, 1))
#dev.off()

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

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

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

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

print(shapiro.test(standardRes))

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