Regression and point-specific p-values (using R for explanation)

Question

Consider the following (in R):

library(MASS)
plot(stack.loss~Air.Flow,data=stackloss)
regression <- rlm(stack.loss~Air.Flow,data=stackloss)
abline(regression)

For each of the points I would like to to test against the null hypothesis that it really is located on the line (rather than where it really is) - in essence quantifying the degree to which these are outliers with respect to the regression model.

Possibly using hopeless search term combinations, I have been unable to identify appropriate methodology. Would it be valid to use for each residual a one-sample t-test (ar equivalent) against the residual population to establish such a measure?

Edit: Further reading seems to indicate that "Tolerance Interval" is what I might be looking for - the "tolerance" R package provides calculations of that. I am however puzzled by the apparent contradictory nature of such intervals for regressions as defined by different publications. In the R package case the tolerance bands seem to track the regression line parallely (see Fig. 12 and below example), while alternative sources such as this operate with bands that are reminiscent of confidence intervals in shape (see Fig. 4). The latter seems more intuitive, but other's opinions would be appreciated.

For each point in a sample I am actually looking for "which x/95% tolerance band is running through this point" ...

Here's how usage of "tolerance" might look like (adding to above code) for 95%/95% 2-sided nonparametric regression tolerance bounds:

library(tolerance)
tol.bounds <- npregtol.int(x=stackloss$Air.Flow,        
    y=stackloss$stack.loss, y.hat=regression$fit, side=2, alpha=0.05, P=0.95, 
    method="WILKS")
lines(tol.bounds$x,y=tol.bounds$"2-sided.lower",col="red")
lines(tol.bounds$x,y=tol.bounds$"2-sided.upper",col="red")

Answer 1

From your comment I think it would be most usefull to look at the IWLS weights:

library(ggplot2)
p <- ggplot(stackloss,aes(x=Air.Flow,y=stack.loss)) + 
  geom_point(aes(colour=regression$w),size=5) + 
  geom_smooth(method="rlm",se=FALSE) +
  scale_colour_gradient(low="red",high="black",name="Huber weights")
print(p)

Answer 2

The most straightforward way to "quantify outlierishness" would be through the examination of the residuals themselves. (Eg. a qq-plot) But OK, let's see what you want. So a p-value is "the probability, given a null hypothesis for the probability distribution of the data, that the outcome would be as extreme as, or more extreme than, the observed outcome" (MacKay Chapt. 37).

So what is your null hypothesis usually? For starters one would assume a expected value of 0 which would mean that there is no bias in the estimator and that secondly the standard deviation of the residuals is $\sigma$, $\sigma$ being your residual standard error (1.915 in your case). So what we want is really the chance of observing a residual as large as the one we observed given the probability distribution we assumed for the data (where the data are now the residuals themselves). Thus assuming $\epsilon \sim N(0,\sigma^2)$ we can do something like : as.numeric(round( digits= 5,1- pnorm( mean=0, sd=1.915, q= abs(regression_residuals)))) and get a p-value (or z-value in this case to be a bit more exact). But here is the "catch", this is all rubbish. It is rubbish because we have used "robust regression" which explicitly encodes the assumption that our residuals are not normally distributed and most possibly exhibit heteroskedasticity (Wikipedia - Robust Regression).

So back square one. As some commenters have already noted: Are asking the correct question? I think you are not, at least not totally. If want to find if indeed some point really is "bluntly outlierish" and that "this is a good thing" then you might as well use a linear model, encode your normality assumptions by using lm() and then look at the outliers . Clearly those residuals will be strongly correlated with the Huber weights (regression$w) you are having here, and with robust residuals you are looking at right now also. I believe you need to work a bit on the formulation of your model. What exactly are you asking? Here is a brief explanation of how robust regression works against normal linear regression. Are you sure you want to use robust regression? Given the assumptions you are will to built in your model you can answer the question that stem from them, not the other way around. :)