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

plot(USPop$year, USPop$population, xlab='year', ylab='population')
grid()

# Use a self starter function for this functional form:
#
USPop.m1 = nls(population ~ SSlogis(year, Asym, xmid, scal), data=USPop)
print(summary(USPop.m1))


# Convert the coefficients in the functional form given by the self-starter function into the functional form given in this problem:
#
beta1 = as.double(coefficients(USPop.m1)[1])
beta2 = as.double(coefficients(USPop.m1)[2] /coefficients(USPop.m1)[3])
beta3 = 1/as.double(coefficients(USPop.m1)[3])
print('Estimates of beta coefficients:')
print(c(beta1, beta2, beta3))

# Now that we have a model, lets study the fitting assumptions following the steps discussed in this chapter of the book:
#
plotFit(USPop.m1, main='USPop data and fit')
grid()

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

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

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

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

print(shapiro.test(standardRes))

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