Linear model diagnostic plots i.e. hacking plot.lm(...)

plot.lm.fun.R

##
## Code was compiled by Paul Julian
## contact info: pjulian@ufl.edu
##
## Example of how to hack plot.lm

#devtools::install_github("SwampThingPaul/AnalystHelper")
library(AnalystHelper)

#for color maps
library(viridisLite)
library(wesanderson)

## Functions

## Use AnalystHelper::axis_fun(...) or run this function
axis_fun=function(side,at,at2,labels,cex.axis=1,line=-0.25,lwd=1,maj.tcl=-0.6,min.tcl=-0.3,las=1,axisLine=0){
  axis(side,line=line,at=at,labels=labels,las=las,tcl=maj.tcl,lty=0,cex.axis=cex.axis);
  axis(side,at=at,labels=F,las=las,tcl=maj.tcl,lwd=lwd,line=axisLine);
  axis(side,at=at2,labels=F,tcl=min.tcl,lwd=lwd,line=axisLine)
}


qq.function=function(y){
  n=length(y)
  r=order(order(y))
  if(n>10){p=(r-1/2)/n}else{p=(r-3/8)/(n+1/4)}
  qqnorm.val=qnorm(p)
  return(qqnorm.val)
}

# using R data in the ?plot.lm example for demonstration
lm.SR <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings)

par(family="serif",mar=c(2,3,1,0.5),oma=c(2,1,0.25,0.75));
layout(matrix(1:4,2,2))
plot(lm.SR)


par(family="serif",mar=c(2,3,1,0.5),oma=c(2,1,0.25,0.75));
layout(matrix(1:4,2,2))

cols=rainbow(length(lm.SR$fitted.values))
ylim.val=c(-10,10);by.y=5;ymaj=seq(ylim.val[1],ylim.val[2],by.y);ymin=seq(ylim.val[1],ylim.val[2],by.y/2)
xlim.val=c(5,17);by.x=2;xmaj=seq(xlim.val[1],xlim.val[2],by.x);xmin=seq(xlim.val[1],xlim.val[2],by.x/2)
plot(lm.SR$fitted.values,lm.SR$residuals,ylim=ylim.val,xlim=xlim.val,type="n",axes=F,ylab=NA,xlab=NA)
abline(h=ymaj,v=xmaj,lty=3,col="grey80")
abline(h=0)
with(lm.SR,points(fitted.values,residuals,pch=21,bg=adjustcolor(cols,0.5),col="grey",lwd=0.1,cex=1.25))
fit.vals=lm.SR$fitted.values
res.vals=lm.SR$residuals
with(lowess(fit.vals,res.vals),lines(x,y,col="red",lwd=2))
axis_fun(1,line=-0.5,xmaj,xmin,xmaj);axis_fun(2,ymaj,ymin,format(round(ymaj,2)));box(lwd=1)
mtext(side=2,line=2.5,"Residuals")
mtext(side=1,line=1.75,"Fitted Values")

cols=viridisLite::viridis(length(lm.SR$fitted.values))
ylim.val=c(-2,3);by.y=1;ymaj=seq(ylim.val[1],ylim.val[2],by.y);ymin=seq(ylim.val[1],ylim.val[2],by.y/2)
xlim.val=c(-2,2);by.x=1;xmaj=seq(xlim.val[1],xlim.val[2],by.x);xmin=seq(xlim.val[1],xlim.val[2],by.x/2)
rstd=rstandard(lm.SR)
qq.x=qq.function(lm.SR$residuals)
plot(rstd~qq.x,ylim=ylim.val,xlim=xlim.val,type="n",axes=F,ylab=NA,xlab=NA)
abline(h=ymaj,v=xmaj,lty=3,col="grey80")
points(qq.x,rstd,pch=21,bg=adjustcolor(cols,0.5),col="grey",lwd=0.1,cex=1.25)
abline(0,1,lty=3)
axis_fun(1,line=-0.5,xmaj,xmin,xmaj);axis_fun(2,ymaj,ymin,ymaj);box(lwd=1)
mtext(side=2,line=2.5,"Standardized Residuals")
mtext(side=1,line=1.75,"Theoretical Quantiles")

cols=wesanderson::wes_palette("Zissou1",length(lm.SR$fitted.values),"continuous")
ylim.val=c(0,2);by.y=1;ymaj=seq(ylim.val[1],ylim.val[2],by.y);ymin=seq(ylim.val[1],ylim.val[2],by.y/2)
xlim.val=c(5,17);by.x=2;xmaj=seq(xlim.val[1],xlim.val[2],by.x);xmin=seq(xlim.val[1],xlim.val[2],by.x/2)
plot(sqrt(abs(rstd))~fit.vals,ylim=ylim.val,xlim=xlim.val,type="n",axes=F,ylab=NA,xlab=NA)
abline(h=ymaj,v=xmaj,lty=3,col="grey80")
points(fit.vals,sqrt(abs(rstd)),pch=21,bg=adjustcolor(cols,0.5),col="grey",lwd=0.1,cex=1.25)
with(lowess(fit.vals,sqrt(abs(rstd))),lines(x,y,col="red",lwd=2))
axis_fun(1,line=-0.5,xmaj,xmin,xmaj);axis_fun(2,ymaj,ymin,ymaj);box(lwd=1)
mtext(side=2,line=2,expression(sqrt("Standardized Residuals")))
mtext(side=1,line=1.75,"Fitted Values")

cols=viridisLite::magma(length(lm.SR$fitted.values))
ylim.val=c(-2,3);by.y=1;ymaj=seq(ylim.val[1],ylim.val[2],by.y);ymin=seq(ylim.val[1],ylim.val[2],by.y/2)
xlim.val=c(0,0.5);by.x=0.1;xmaj=seq(xlim.val[1],xlim.val[2],by.x);xmin=seq(xlim.val[1],xlim.val[2],by.x/2)
lev=hatvalues(lm.SR)
plot(rstd~lev,ylim=ylim.val,xlim=xlim.val,type="n",axes=F,ylab=NA,xlab=NA)
abline(h=ymaj,v=xmaj,lty=3,col="grey80")
abline(h=0)
points(lev,rstd,pch=21,bg=adjustcolor(cols,0.5),col="grey",lwd=0.1,cex=1.25)
with(lowess(lev,rstd),lines(x,y,col="red",lwd=2))
inf=lm.influence(lm.SR)
hh=seq(min(range(inf)[1],range(inf)[2]/100),xlim.val[2]+(xlim.val[2]*0.5),length.out=101)
hh=hh[hh>0]
crit.cook=0.5;
cl.h=sqrt(crit.cook*length(coef(lm.SR))*(1-hh)/hh)
lines(hh,cl.h,lty=2,col=2)
lines(hh,-cl.h,lty=2,col=2)
axis_fun(1,line=-0.5,xmaj,xmin,format(xmaj));axis_fun(2,ymaj,ymin,ymaj);box(lwd=1)
mtext(side=2,line=2,"Standardized Residuals")
mtext(side=1,line=1.75,"Leverage")

##
## End