diagANOVA <- function(m,raw.res=FALSE,tests=FALSE) # Perform diagnostics on an ANOVA { # Set up graophing window par(mfrow = c(2,2)) # Make 2 rows and columns in the figure window if (raw.res) { r = residuals(m) # raw residuals # Check: normality hist(r,freq=FALSE,main="Density Plot",xlab="Raw Residuals") lines(density(r)) boxplot(r,main="Boxplot",ylab="Raw Residuals") xy=qqnorm(r,main="Normal Plot",ylab="Raw Residuals") qqline(r) # Check: homogeneity of variance of raw res plot(fitted(m),r,main="Res Vs. Fits",ylab="Raw Residuals") # raw residuals abline(0,0) } else { r = resid(m,type="pearson") # Latest version of R requires this #r = resid(m,type="p") # standardized residuals # Check: normality of st res hist(r,freq=FALSE,main="Density Plot",xlab="St Residuals") lines(density(r)) boxplot(r,main="Boxplot",ylab="St Residuals") xy=qqnorm(r,main="Normal Plot",ylab="St Residuals") qqline(r) # Check: homogeneity of variance of st res plot(fitted(m),r,main="Res Vs. Fits",ylab="St Residuals") # raw residuals 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 ks.test(r,"pnorm") } # Set plotting window back to normal par(mfrow = c(1,1)) }