#-----------------------------------------------#
#--- Graphing Linear Regression Results in R ---#
#-----------------------------------------------#

# Here is a very basic regression example, generating four useful plots

# Let's make up some data with normal errors
x = seq(10,40,0.5)
y_ideal = x^2
N = length(x)
sd = 20
y = rnorm(N,y_ideal,sd)  # rnorm(n,mean,sd)

# perform the linear regression
model = lm(y~I(x^2))  # I(x^2) means include the x^2 term (only) in the model
summary(model)

# Now generate a series of four plots about the data and model
par(mfrow=c(2,2))

# plot the data and the best fit model
plot(y ~ x, cex = 1.3, main = "Data and Model")  
#now add the model line to the plot
lines(x, predict(model), col = "red")

# plot the residuals
library(MASS)  #for the studres function
esr = studres(model) # externally studentized residuals (from MASS library)
plot(esr~predict(model), main = "Residuals")

# Normal Probability Plot
qqnorm(esr)
qqline(esr, col = "red")

# Let's create a Williams Graph
norm_leverage = (length(esr)/2)*hatvalues(model)
plot(abs(esr) ~ norm_leverage, xlab = "leverage*n/p", ylab = "|esr|", main ="Williams Graph")

# add the critical lines to the Williams Graph
N = length(esr)
p = length(model$coefficients)
alpha = 0.01
tinv = qt(alpha/2/N,N-p-1)
Grubbs_crit = (N-1)/sqrt(N)*sqrt(tinv^2/(N-3+tinv^2))
abline(h=Grubbs_crit, col="red") # draw a red horizontal line at Gcrit
abline(v=2, col="red") # draw a red vertical line at 2






# More on plotting

plot(y ~ x, cex = 1.3)  # cex = 1.3 makes the data points 1.3X bigger

# bonus topic: LOWESS = Locally Weighted Scatterplot Smoothing
# this creates a nonparametric "smooth" line through the data
# it is usefull for seeing trends in the data (before modeling)
lines(lowess(x,y), col="blue")

# here is an enhanced plotting function, where the LOWESS line and/or
# a straight regression line can be automatically plotted as well
library(car)
scatterplot(y ~ x, smoother="FALSE", reg.line="FALSE")

# Other plotting systems in R are lattice and ggplot