Although the second order mixed differences
$$D_{12}f(x_1, x_2, y_1, y_2) = f(x_1,y_1)+f(x_2,y_2)-f(x_1,y_2)-f(x_2,y_1)$$
provide useful information, by themselves they do not have enough power to discriminate between slightly noisy additive functions and non-additive functions.
The trick is to compare $D_{12},$ suitably transformed to have good statistical properties, to some measure of the spread of $f,$ without assuming anything at all about $f$. This answer gives the details and provides evidence that this procedure can work well even for small grids and relatively large amounts of noise.
A function $f$ of two variables $x$ and $y$ is "additively separable" ("additive" for short) when there exist functions $g$ and $h$ such that $f(x,y) = g(x)+h(y).$ A "noisy" additive function includes independent additive random errors $\varepsilon(x,y)$, so that the full model is
$$f(x,y) = g(x) + h(y) + \varepsilon(x,y).$$
Suppose we have observed $f(x_i, y_j)$ for all $mn$ combinations of $m$ distinct values $x_1, x_2, \ldots, x_m$ and $n$ distinct values $y_1, y_2, \ldots, y_n.$ (In full generality--assuming absolutely nothing about $f$, not even whether it is continuous--it matters not what the values $x_i$ or $y_j$ are, because they can be subsumed in the definitions of $g$ and $h$, respectively. Furthermore, the ordering of the $x_i$ and $y_j$ are immaterial: the definition of additivity makes no use of that.)
The first idea is that additive functions $f$ can be detected by comparing $D_{12}$ to zero, because
$$\eqalign{
D_{12}f(x_1, x_2, y_1, y_2) &= g(x_1)+h(y_1)+\varepsilon(x_1,y_1) + g(x_2)+h(y_2)+\varepsilon(x_2,y_2) \\
&- g(x_1) - h(y_2) - \varepsilon(x_1,y_2) - g(x_2) - h(y_1) - \varepsilon(x_2, y_1) \\
&=\varepsilon(x_1,y_1) +\varepsilon(x_2,y_2)- \varepsilon(x_1,y_2)- \varepsilon(x_2, y_1)\\
&= D_{12}\varepsilon(x_1, x_2, y_1, y_2),
}$$
which is the sum of four independent errors. Provided these "noise terms" tend to be small, all mixed differences ought to be small when applied to additive functions. When applied to non-additive functions, the mixed differences ought to be greater.
Difficulties arise for several reasons:
Because the orders of the $x_i$ and $y_j$ do not matter, one can compute mixed differences for all $m(m-1)/2$ distinct pairs of $x$'s and all $n(n-1)/2$ distinct pairs of $y$'s, even though only $(m-1)(n-1)$ of them are (linearly) independent. Which differences should be used? (I propose using them all.)
Because adding and subtracting values of a function throughout its domain can cause relatively large values to be combined, the variation of a function's values can swamp any signal in the mixed differences.
It is tempting instead to fit a model of the form
$$\mathbb{E}(f(x_i, y_j)) = \alpha_i + \beta_j$$
and then study its residuals--but is there any way that can help and, if so, how should the model be fit? Can we do this in a robust fashion that is resistant not only to outlying values of the $\varepsilon(x_i,y_j)$ but also to possibly extreme variation in $g$ and $h$ themselves?
A solution I have hit upon after much experimentation addresses these issues by constructing a set of paired values $(a_k, b_k)$ derived from the observations of $f$. The indexes $k$ designate $2\times 2$ blocks in the data and therefore correspond to all tuples $(i_1, i_2, j_1, j_2)$ for which $1 \le i_1 \lt i_2 \le m$ and $1 \le j_1 \lt j_2 \le n$. For such a tuple
My proposal is to normalize all the values of $a_k$ and $b_k$ (computed for all possible blocks) relative to the maximum value of $a_k$ and look at their scatterplot. For an additive function with no error, the $b_k$ are always zero and the scatterplot will be perfectly horizontal with no correlation. When there is error, suppose the typical variance of the $\varepsilon(x,y)$ is $\sigma^2$. Then the variance of the mixed differences will be $4\sigma^2$ while the $a_k$ will extend from near $0$ to the full range of $f$. The square root of the absolute values of those first differences therefore will typically be around $2\sigma$ or so; provided this is small compared to the range of $f$, the scatterplot ought to be perfectly level--at least for the larger values of $a_k$. (For the smaller values, $b_k$ will be constrained to be small no matter what.) This is the second key idea and it seems to work.
This reasoning yields an additivity diagnostic plot: the normalized scatterplot of $(a_k, b_k)$ for noisy additive functions will
Tend to be near the horizontal axis and
Rise from the origin to a horizontal value well below $1$.
By contrast, the scatterplot for non-additive functions will tend to produce larger values of $b_k$ when the $a_k$ are large, because large ranges of functional values ought to exhibit relatively large deviations from additivity. Consequently for many non-additive functions we expect the diagnostic plot to be nearly diagonal and linear, rising from the origin $(0,0)$ nearly to the opposite corner at $(1,1).$
Pictures are worth more than these words, so let's look at some examples. In the figure I have constructed three functions whose values range approximately from $-1$ to $1$, sampled them on a $7$ by $7$ grid covering the square $[-1,1]\times [-1,1]$, generated a fixed normally-distributed noise vector of $\varepsilon(x_i, y_j)$ (with unit variance), and added various multiples $\sigma = 0, 0.1, 0.3$ of this noise vector to the three samples. Each line of the figure corresponds to one function. It displays the additivity diagnostic plots for the three versions--one with no noise and two with increasing noise--then shows the functional values on a grid, and finally shows the noisiest values on a grid. The noise in that latter version has a standard deviation equal to $1/6$ the original range of the function, which is pretty big. Notice, too, that a $7$ by $7$ grid has only $(7-1)(7-1) = 36$ degrees of freedom (independent mixed differences), which is a pretty small amount of data concerning these functions: all in all, this is a fairly difficult test of the procedure.
The red lines are lowess smooths of the scatterplots.
It should be evident which of these functions is additive: only the first one, $f(x,y) = (x^2+y^2)/2.$ The other two, $xy^3$ and $\cos(\pi x) \sin(2\pi y),$ are decidedly non-additive. They are distinguished from the first in all three additivity diagnostic plots, which closely follow the preceding descriptions. In short, because the plots on the first row (a) are low in height, (b) exhibit low correlation, and (c) have smooths that level off for larger values of $a_k$ (although the one for $\sigma=0.3$ is marginal), they signal that the first function may be additive. Because the smoothed plots on the second and third rows all closely approximate the major diagonal of the square, they signal non-additivity of the data.
The visual similarities among the noisy versions of these functions (in the rightmost column) attest to the utility of this diagnostic plot.
Here is the R
code used to construct the figure. Do not try to apply it to large grids of data! Because all $2 \times 2$ blocks of the grid are evaluated, the computational effort and storage requirements are both $O(m^2 n^2),$ which is pretty bad. For instance, a medium-sized grid of (say) $m=50$ and $n=100$ has $6063750$ blocks to evaluate! (For anyone lucky enough to have such a rich dataset, consider subsampling the blocks rather than exhaustively evaluating all of them.)
#
# Compute all rectangular differences.
# More generally, `fun` will be applied to all 2 x 2 sub-blocks of a.
#
diff.rect <- function(z, a=matrix(c(1,-1,-1,1), 2), fun) {
if (missing(fun)) fun <- function(x) sum(x*a)
m <- dim(z)[1]; n <- dim(z)[2]
ii <- combn(1:m, 2); jj <- combn(1:n, 2)
sapply(1:dim(jj)[2], function(j) {
sapply(1:dim(ii)[2], function(i) fun(z[ii[,i], jj[,j]]))
})
}
#
# Compare ranges to differences.
#
additive.test <- function(z, ...) {
z.diff <- diff.rect(z)
z.range <- diff.rect(z, fun=function(x) diff(range(x)))
scale <- max(z.range)
z.range <- z.range / scale
z.diff <- sqrt(abs(z.diff) / scale)
plot(z.range, z.diff, cex=0.7, col="#40404080", ylim=c(0,1), ...)
abline(h = median(z.diff), col="#4040ff80", lty=2)
lines(lowess(z.range, z.diff), lwd=2, col="Red")
}
m <- 7; n <- 7
x <- seq(-1, 1, length.out=m)
y <- seq(-1, 1, length.out=n)
#
# Plot differences against ranges.
#
require(raster)
image <- function(x, ...) plot(raster(x), asp=dim(x)[1]/dim(x)[2],
axes=FALSE, frame.plot=FALSE, ...)
par(mfrow=c(3, 5))
set.seed(17); error <- rnorm(m*n)
for (f in c(function(x,y) (x^2 + y^2)/2,
function(x,y) x * y^3,
function(x,y) cos(pi*x) * sin(2*pi*y))) {
f.body <- deparse(f)[2]
f.xy <- outer(x, y, f)
for (sigma in c(0, 1/10, 3/10)) {
g.xy <- f.xy + error * sigma
additive.test(g.xy, main=paste("sd =", round(sigma, 2)),
xlab="", ylab="", cex.main=9/10)
}
image(f.xy, main=f.body, cex.main=9/10)
image(f.xy + error * 3/10, main="SD = 0.3", cex.main=9/10)
}