[R] Plot residuals against standard normal distribution

[Mike Miller]

On Sat, 7 Feb 2015, Mike Miller wrote:

> res <- residuals( model )
>
> resStd <- ( res - mean( res, na.rm=TRUE ) ) / sd( res, na.rm=TRUE )
>
> Another issue is how to make the theoretical quantiles for the normal 
> distribution.  There are a few methods:
>
> https://www.statsdirect.com/help/data_preparation/normal_scores.htm
>
> I usually use Blom:
>
> r <- rank( resStd )
> c <- 3/8
> N <- sum( !is.na( resStd ) )
> resNorm <- qnorm( ( r - c ) / ( N - 2*c + 1 ) )
> resNorm[ is.nan( resNorm ) ] <- NA
>
> Then you could plot it directly:
>
> plot(resNorm, resStd)
>
> When we use qqnorm() in R, it looks like R is using a Blom method with 
> c=1/2 instead of c=3/8.  I believe Blom recommended 3/8 and programs 
> that offer Blom normal scores use c=3/8.


I don't know if that was off-track because the OP was asking about 
density, but he also was asking about checking the distribution of 
residuals, so maybe this is appropriate.

I should add, if you don't mind using R's c=1/2, you can get the normal 
scores very quickly this way:

resNorm <- qqnorm( residuals( model ), plot.it=FALSE )$x

Apparently, 11 years ago R was using c=3/8 in qqnorm(), so I guess it 
changed.  Nordheim, Clayton and Yandell wrote about it in this document 
dated September 9, 2003:

https://www.stat.wisc.edu/~yandell/st571/R/append8.pdf

It is definitely using c=1/2 today.  I don't know where that is 
documented.

When I do a QQ-plot of uniform p-values, I like to add a confidence region 
with these lmits:

   qb95 <- qbeta(.95,1:N,N+1-(1:N))
   qb05 <- qbeta(.05,1:N,N+1-(1:N))

If we have N observations from a normal distribution with unknown mean and 
variance, can we create some kind of analogous region on a qqnorm() kind 
of plot?  It seems like there should be a way to get at least an 
approximate result using beta and t distributions, probably building on 
the qbeta() code above.

Mike