### Solved: Recreate Minitab Normal Probability Plot in R

I'm taking a stats class where everything is done in Minitab, which is Teh Suck because Minitab only runs on Windows and I No Haz Windows. I've been successfully using R as a substitute, which has turned out to be a non-issue for the most part. However...

Minitab produces a very specific type of probability plot, the "Normal Probability Plot", which has no near-equivalent in R. Thus I embarked on a journey to try to recreate the damn thing. There are lots of threads all over the intertoobs that solve parts of the problem (which turns out to be fairly involved), but I couldn't find anything which puts all the pieces together into a single, ready-made solution. What follows is a pretty close recreation; I just had to move on to the rest of my life before I could get the legend right.

Presented in the hopes that it'll save other people the effort:

minitab_normal_prob_plot <- function(data, x_label) { # The labels for the y-axis, corresponding to percentiles y_axis_labels = c(1,5,10,20,30,40,50,60,70,80,90,95,99) # Lengths, mean, and sd of data n = length(data) my_mean = mean(data) my_sd = sd(data) ### Set up the y-axis values # Translate labels to decimal percentages percentages = y_axis_labels / 100 # Convert percentages to z-values and shift so that all values are >= 0 y_axis_points = qnorm(percentages) y_shift = y_axis_points[1] y_axis_points = y_axis_points - y_shift # The minimum and maximum y values y_min = y_axis_points[1] y_max = y_axis_points[length(y_axis_points)] ### Calculate the main data set # x and y values per http://en.wikipedia.org/wiki/Normal_probability_plot. x_data_points = sort(data) data_percents = c() for(i in 1:n) { if (i == 1) { data_percents[i] = 1 - 0.5^(1/n) } else if (i == n) { data_percents[i] = 0.5^(1/n) } else { data_percents[i] = (i - 0.3175)/(n+0.365) } } ### Trend line calculation # Project a line represented expected distribution values on assumption # that data is normal. trend_x0 = qnorm(percentages[1], mean = my_mean, sd = my_sd) trend_x1 = qnorm( percentages[length(percentages)], mean = my_mean, sd = my_sd ) # Convert percents to z-values and shift as before y_data_points = qnorm(data_percents) - y_shift ### Set up the envelope # Stolen from # http://stackoverflow.com/questions/3929611/recreate-minitab-normal-probability-plot library(MASS) qprobs<-qnorm(percentages) fd<-fitdistr(data, "normal") #Maximum-likelihood Fitting of Univariate Dist from MASS xp_hat<-fd$estimate[1]+qprobs*fd$estimate[2] #estimated perc. for the fitted normal v_xp_hat<- fd$sd[1]^2+qprobs^2*fd$sd[2]^2+2*qprobs*fd$vcov[1,2] #var. of estimated perc xpl<-xp_hat + qnorm(0.025)*sqrt(v_xp_hat) #lower bound xpu<-xp_hat + qnorm(0.975)*sqrt(v_xp_hat) #upper bound ### Set up the x-axis x_min = min(c(data, trend_x0, trend_x1, xpl, xpu)) x_max = max(c(data, trend_x0, trend_x1, xpl, xpu)) ### Plot it all # Data set. Points plotted twice due to keep them from getting clobbered by # white rectangle. par(bg = "beige") plot( x_data_points, y_data_points, xlim = c(x_min, x_max), ylim = c(y_min, y_max), axes = FALSE, ylab = "Percent", xlab = x_label, pch = 16, col = "red", main = paste("Probability Plot of", x_label,"\nNormal - 95% CI") ) rect(par("usr")[1], par("usr")[3], par("usr")[2], par("usr")[4], col = "white") points(x_data_points, y_data_points, pch = 16, col = "red") # Trend line segments(trend_x0, y_min, trend_x1, y_max, col = "blue") # Lower and upper bounds lines(xpl, y_axis_points, col = "blue") lines(xpu, y_axis_points, col = "blue") # Y-axis gridlines for (i in 1:length(y_axis_points)) { abline(h = y_axis_points[i], col = "gray", lty = 2) } # Axes axis(1) axis(2, at = y_axis_points, labels = y_axis_labels) # Box and x-grid box() grid(ny = NA, lty = 2) # Legend legend( "topright", c( paste("Mean", my_mean, sep = " "), paste("StDev", my_sd, sep = " ") ), bg = "white" ) }