Use a robust additive model to fit the data and then explore its residuals.
An "additive model" approximates the value $y_{ij}$ at row $i,$ column $j$ as the sum of a row term $\rho_i$ and a column term $\gamma_j,$ with any difference known as the "residual" $r_{ij}.$ Thus,
$$r_{ij} = y_{ij} - (\rho_i + \gamma_j).$$
Relatively large residuals will crop up in any cell that is inconsistent with the typical values suggested by its row and column.
A "robust" method estimates the row and column terms in a way that is insensitive to values that would yield large residuals. There are many robust methods. For instance, John Tukey described a simple "median polish" in his EDA book in which you alternately estimate the terms as the medians of their rows and columns, make suitable adjustments, and repeat until the fit works well. I like this, but it's not readily available in software (because convergence is not guaranteed in all cases).
To fit the first example, I used the rlm
("robust linear model") function in the MASS
library for R
.
"Exploration" is usually most effective when done visually, by plotting suitable graphics.
Simple, standard graphical tools can help identify unusually large residuals. For instance, here is a boxplot of these 20 residuals that clearly shows the two extreme values--and just how extreme they are:
For two-way tables like your examples, map the residuals.
(For reference I have posted the original data values.) By making small-sized residuals nearly transparent, the largest ones really stand out in this map: the large negative residual (corresponding to the original value of 1 in row 1, column D) and the large positive residual (corresponding to the original value of 9 in row 4, column C).
To see how insightful this can be, apply it to your second example. You will be surprised at the entries with the most extreme residuals! (The value of 6 in the third row is not an outlier. That's because it is perfectly consistent with the large values that occur in column A.)
Here is the R
code.
# X <- data.frame(y = c(
# 10, 9, 8, 8, 5,
# 7, 3, 3, 2, 3,
# 6, 1, 1, 1, 1),
# Row = factor(rep(1:3, each=5), levels=3:1),
# Column = rep(c("A","B","C","D","E"), 3))
X <- data.frame(y = c(
10, 10, 10, 1, 10,
10, 9, 8, 8, 10,
1, 2, 2, 1, 2,
1, 1, 9, 1, 1),
Row = factor(rep(1:4, each=5), levels=4:1),
Column = rep(c("A","B","C","D","E"), 4))
library(MASS)
fit <- rlm(y ~ Column + Row - 1, X)
summary(fit)
# plot(fit)
X$Residual <- residuals(fit)
boxplot(X$Residual, horizontal=TRUE)
library(ggplot2)
library(scales)
ggplot(X) +
geom_raster(aes(Column, Row, fill=Residual, alpha=abs(Residual))) +
geom_text(aes(Column, Row, label=y), color="Gray") +
scale_fill_gradient2(low=muted("blue"), high=muted("red")) +
ggtitle("First example")