diagRES <- function(m,res="n",num.panes=2,tests=FALSE)
# Perform diagnostics on residuals from a statistical model
# The input RES specifies the type of residuals to display and possibly test.
# 	RES = "r" indicate raw residuals (data - predicted)
#	RES = "n" (only works for nlme models) indicate residuals normalized by the serial correlation structure of the model (required 
#               to assess model fit of ARIMA models)
#	RES = "p" indicate pearson or standardized residuals, the default outputted by R's nlme and lme4 models
{
# Set up graphing window
if (num.panes==4)
{par(mfrow = c(2,2)) # Make 2 rows and columns in the figure window
} else {
  par(mfrow=c(1,2))
}

if (res=="r")
{
r = residuals(m)   # raw residuals
lab = "Raw Res"
}   # if "r"
else
{
if (res=="n")
{
r = resid(m,type="normalized")
lab = "Normalized Res"
}  # if "n"
else 
{
if (res=="p")
{r = resid(m,type="pearson")   # Latest version of R requires this
#r = resid(m,type="p")  # standardized residuals
lab = "Pearson Res"
}   # if "p"
}   # else "n"
}   # else "p"


# Check: normality 
if (num.panes==4)
{hist(r,freq=FALSE,main="Density Plot",xlab=lab)
lines(density(r))
boxplot(r,main="Boxplot",ylab=lab)
}
xy=qqnorm(r,main="Normal Plot",ylab=lab)
qqline(r)
# Check: homogeneity of variance of raw res
plot(fitted(m),r,main="Res Vs. Fits",ylab=lab)
abline(0,0)

if (tests)
{
# Check the correlation in the qq-plot
print(sprintf('In this sample of size n=%d, correlation of the residuals in the qq-plot is r=%f',length(r),cor(xy$y,xy$x)))
print("In the following table, if r < critical.r, then the qq-plot suggests the residuals are not normal:")
n=c(5,10,15,20,25,30,40,50,60,75)
critical.r=c(.832,.88,.911,.929,.941,.949,.96,.966,.971,.976)
print(data.frame(n,critical.r))  

# Shapiro Wilks test
print(shapiro.test(r))

# Kolmogorov test - too sensitive
print(ks.test(r,"pnorm"))
}   # if tests

  # Set plotting window back to normal
  par(mfrow = c(1,1)) 
  
}