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

plot(Indometh$time, Indometh$conc, xlab='time', ylab='conc')
grid()

Indometh.1.m1 = nls(conc ~ SSbiexp(time, A1, lrc1, A2, lrc2), data=Indometh, subset=Subject==1)
print(summary(Indometh.1.m1))

# 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_5.eps", onefile=FALSE, horizontal=FALSE)
plotFit(Indometh.1.m1, main='Indometh data and fit(Subject==1)')
grid()
#dev.off()
    
# Residual plot:
#
plot(fitted(Indometh.1.m1), residuals(Indometh.1.m1), xlab='fitted values', ylab='residuals')
abline(h=0, col='black')
grid()

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

# Variance homogeneity:
#
plot(fitted(Indometh.1.m1), abs(residuals(Indometh.1.m1)), xlab='fitted values', ylab='abs residuals')
grid()

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

print(shapiro.test(standardRes))

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