[R] Plot residuals against standard normal distribution

[Mike Miller]

On Mon, 2 Feb 2015, Mikael Olai Milhøj wrote:

> I'm having trouble trying to plot the density of the residuals against 
> the standard normal distribution N(0,1). (I'm trying to see if my 
> residuals are well-behaved).
>
> I know hwo to calculate the standardized residuals (I guess that there 
> may be a simple way using a R function) and then plot this by using the 
> density function
>
> y<-(model$residuals-mean(model$residuals))/sd(model$residuals)
> plot(density(y))
>
> But I don't know how to add the N(0,1) curve. Any suggestions? Thanks in 
> advance

This isn't good for you?

qqnorm( residuals( model ) )

Residuals usually have a zero mean, but I guess you can center, anyway, 
and if there could be NAs, you will need to deal with them:

res <- residuals( model )

resStd <- ( res - mean( res, na.rm=TRUE ) ) / sd( res, na.rm=TRUE )

qqnorm( resStd )

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.

Best,
Mike

-- 
Michael B. Miller, Ph.D.
Minnesota Center for Twin and Family Research
Department of Psychology
University of Minnesota
http://scholar.google.com/citations?user=EV_phq4AAAAJ