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)

source('SSlogLogistic.R')

plot(heartrate$pressure, heartrate$rate, xlab='pressure', ylab='rate')
grid()

heartrate.m1 = nls(rate ~ SSlogLogistic(pressure, b, c, d, e), data=heartrate)
print(summary(heartrate.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_3.eps", onefile=FALSE, horizontal=FALSE)
plotFit(heartrate.m1, main='heartrate data and fit')
grid()
#dev.off()

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

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

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

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

print(shapiro.test(standardRes))

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