# Determining if a function is additive

I have data collected on a grid of values $f(x,y)$. I have a hypothesis that the data is "additively separable", i.e.

$$f(x,y) = h_1(x) + h_2(y) + g(x,y)$$

where $g(x,y)$ is small and only represent the noise in my measurements. Is there a way to test for this? If I had infinite sampling resolution I could measure

$$\frac{\partial^2 f(x,y) }{\partial{x}\partial{y}} = \frac{\partial^2 h_1(x) }{\partial{x}\partial{y}} + \frac{\partial^2 h_2(y) }{\partial{x}\partial{y}} + \frac{\partial^2 g(x,y) }{\partial{x}\partial{y}} = g_{xy}(x,y)$$

Then I could test if for some threshold value $\epsilon$, $|g_{xy}(x,y)| > \epsilon$. However, my data is real, with noise and limited samples. What statistical test could I use to check for additive separability?

• I doubt these are the right tags, but I wasn't sure what to use. Please feel free to correct them. – Hooked Jan 31 '14 at 19:14
• You currently have no noise/variation about your model, which makes me suspect this isn't a statistical problem. With noise, and if you were trying to test if $g=0$ was zero a statistical test might be relatively easy (with some assumptions about the $h$'s). To do much more you're going to need to be precise about what "small" means and likely put some limits on $g$ (such as some kind of assumption about smoothness at least). If you can form estimates with $g$ both restricted to be small and not restricted, you might be able to do some kind of resampling-based test -- perhaps. – Glen_b Jan 31 '14 at 21:43
• @Glen_b it may be that I'm not using the proper phrasings here, I'll take any help I can get. My data is assumed to have noise, and imprecision in measurement which I have some knowledge about via repeated measurements of the same point. Is it enough of a condition that $h,g$ are smooth and continuous (which I suspect is the case for me)? Shouldn't the test tell me what "small" is given a particular confidence rating? – Hooked Jan 31 '14 at 21:56
• If there's noise, why does your model have an equality there without a noise term? Are the $f$ values not observed? The part about small is not something the test tells us. Estimation might tell us something about the size of $g$, but only you can declare it to be sufficiently small by your own criteria. How are we (or some test) to know what you mean when you say 'small'? – Glen_b Jan 31 '14 at 22:01
• @Glen_b One solution is to compute all possible values of $f(x_1,y_1)+f(x_2,y_2)-f(x_1,y_2)-f(x_2,y_1).$ Under the null they should be random with zero mean (although correlated in a predictable way), so we can apply a t-test of the null hypothesis provided we compute the standard error appropriately for the correlation. Notice that estimation of $h_1$ and $h_2$ is unnecessary! – whuber Jan 31 '14 at 22:02

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:

1. 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.)

2. 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.

3. 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

• let $a_k$ be the range of the four values of $f$ found in its block and

• let $b_k$ be the square root of the absolute value of $D_{12}f$ for this block.

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

1. Tend to be near the horizontal axis and

2. 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); n <- dim(z)
ii <- combn(1:m, 2); jj <- combn(1:n, 2)
sapply(1:dim(jj), function(j) {
sapply(1:dim(ii), function(i) fun(z[ii[,i], jj[,j]]))
})
}
#
# Compare ranges to differences.
#
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)/dim(x),
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)
f.xy <- outer(x, y, f)
for (sigma in c(0, 1/10, 3/10)) {
g.xy <- f.xy + error * sigma