yalestats:

#' minitabQQ function
#'
#'This function allows you to visualize the normal probability plots from Minitab.Input two different lists representing variables that you want to compare, and the function takes care of the rest.
#'@param x A list of values that will serve as your explanatory "x" variable.
#'@param y A list of values that will serve as your dependent "y" variables.
#'
#'@keywords minitabQQ
#'@export
minitabQQ <- function(x,y){ #In use, the function will take a single argument: your OLS standardized residuals
  modelResiduals <- rstandard(lm(y~x))
  x1 <- sort(modelResiduals) #Sort the residuals low to high to use as the x variable
  y1 <- qnorm(ppoints(length(x1))) #Make a vector of probabilities of the same length (n) as your residuals, based on a standard normal distribution. Then apply this to the quantile function that answers the question, what is the z-score of the pth quantile (percentile from ppoints) of the normal distribution? The qnorm function expects a mean and sd to be specified. When they are not specified the default is mean=0 and sd=1, which is the case for our standardized residuals. These values will be the y ordinate values on the plot, since these are the z-scores that correspond with our standardized residual values.
  plot(x1,y1, type="n", axes=F, ylab="Percent", xlab="Residual", main="Normal Probability Plot") #Create a plot but suppress the content
  axis(1);box() #Add x axis and a bounding box
  probs <- c(0.01, 0.05, seq(0.1, 0.9, by = 0.1), 0.95, 0.99) #Create a controlled vector of probabilities to use in labeling. Include small and large values to enable better visualization.
  qprobs <- qnorm(probs) #Create z-scores based on standard normal distribution with mean 0 and SD of 1. Answers the question, what is the z-score of the pth quantile of the normal distribution? These are the locations at which we will place labels
  axis(2, at = qprobs, labels = 100 * probs) #plot the y axis labels, scaled to 0-100
  points(x1,y1, pch=16, col="red", cex=1.2) #Add the actual residuals to plot
  #Manually compute linear model equation
  x2 <- quantile(x1, c(0.25, 0.75)) #Create two percentile locations based on x data
  y2 <- qnorm(c(0.25, 0.75)) #Create two percentile locations based on standard normal distribution
  slope <- diff(y2)/diff(x2) #Change in y over change in x = slope
  intercept <- y2[1] - slope * x2[1] #Compute intercept for the model
  abline(intercept, slope, col = "blue") #Plot line using intercept and slope
}