# RPD -- Example 11.17 -- Diagnostics: Linthurst Data

library(car)  # needed for leverage plots and variance inflation factors;
              # may need to install this:  install.packages("car")

linthurst <- read.table("linthurst.dat", row.names=1,
                        col.names=c("obsnum","loc","type","biomass",
                                    "salinity","pH","K","Na","Zn"))

linthurst.model <- lm(biomass ~ salinity + pH + K + Na + Zn, data=linthurst)

summary(linthurst.model)
anova(linthurst.model)

# Can get some basic diagnostic plots with:  plot(linthurst.model)


# Sec. 11.4

linthurst.fitted <- fitted(linthurst.model)        # predicted values
linthurst.resid <- residuals(linthurst.model)      # ordinary residuals
linthurst.standard <- rstandard(linthurst.model)   # standardized residuals
linthurst.student <- rstudent(linthurst.model)     # Studentized residuals
linthurst.h <- hatvalues(linthurst.model)          # leverages
linthurst.cook <- cooks.distance(linthurst.model)  # Cook's D
linthurst.dffits <- dffits(linthurst.model)        # DFFITS
linthurst.dfbetas <- dfbetas(linthurst.model)      # DFBETAS
linthurst.covratio <- covratio(linthurst.model)    # COVRATIO

print(data.frame(linthurst.fitted, linthurst.resid, linthurst.standard,
                 linthurst.cook), digits=3)
print(data.frame(linthurst.student, linthurst.h, linthurst.covratio,
                 linthurst.dffits), digits=3)
print(linthurst.dfbetas, digits=3)


# Sec. 11.4.1

plot(linthurst.fitted, linthurst.resid)  # plot of e versus Yhat
abline(h=0, lty=2)

dev.new()
qqnorm(linthurst.standard)  # normal probability plot of standardized residuals

dev.new()
plot(linthurst$salinity, linthurst.resid)
abline(h=0, lty=2)

dev.new()
plot(linthurst$pH, linthurst.resid)
abline(h=0, lty=2)

dev.new()
plot(linthurst$K, linthurst.resid)
abline(h=0, lty=2)

dev.new()
plot(linthurst$Na, linthurst.resid)
abline(h=0, lty=2)

dev.new()
plot(linthurst$Zn, linthurst.resid)
abline(h=0, lty=2)


  # _rescaled_ partial regression leverage plots

dev.new()
leverage.plots(linthurst.model, "(Intercept)", identify.points=FALSE)

dev.new()
leverage.plots(linthurst.model, "salinity", identify.points=FALSE)

dev.new()
leverage.plots(linthurst.model, "pH", identify.points=FALSE)

dev.new()
leverage.plots(linthurst.model, "K", identify.points=FALSE)

dev.new()
leverage.plots(linthurst.model, "Na", identify.points=FALSE)

dev.new()
leverage.plots(linthurst.model, "Zn", identify.points=FALSE)

  # setting identify.points=TRUE enables interactive point identification;
  # omitting the name of the term requests interactive menu selection


# Sec. 11.4.2

dev.new()
plot(linthurst.h, type="h")  # leverages
abline(h=0.267, lty=2)

dev.new()
plot(linthurst.cook, type="h")  # Cook's D
abline(h=0.0889, lty=2)

dev.new()
plot(linthurst.dffits, type="h")
abline(h=c(0.73,-0.73), lty=2)

dev.new()
plot(linthurst.covratio)
abline(h=c(1.4,0.6), lty=2)


# Collinearity Diagnostics

vif(linthurst.model)  # variance inflation factors