Variance components and generalized Sobol' indices

Introduction

Sobol’ indices are certain sums of variance components in an ANOVA decomposition. They are used to understand the importance of various subsets of variables in global sensitivity analysis. Saltelli et al., (2008) give an extensive introduction to Sobol’ indices and variance based methods in general for investigating computer models. Linear combinations of Sobol’ indices are also used to measure the effective dimension of functions for quasi-Monte Carlo integration.

This article reviews Sobol’ indices for a statistical audience, relating them to well known ideas in experimental design, particularly crossed random effects. Moving from physical experiments to computer experiments brings important changes in both the costs and goals of the analysis. In physical experiments one may be interested in all components of variance, or at least all of the low order ones. In computer experiments interest centers instead on sums of variance components. While the ANOVA for computer experiments is essentially the same as that for physical ones, the experimental designs are different.

Of particular interest are what are known as ‘fixing methods’ for estimation of Sobol’ indices. These evaluate the function at two points. Those two points have identical random values for some of the input components (the ones that are ‘fixed’) but have independently sampled values for the other components. Sample variances and covariances of point pairs are then used to estimate the Sobol’ indices.

As a basic example, let $f$ be a deterministic function on $^{5}$. One kind of Sobol’ index estimate takes a form like

An example of the second kind of Sobol’ index is

The sampling design to estimate (2) is that same as that for (1), but the quantity estimated is now the sum of all variance components that involve any of the first three variables. The difference between (1) and (2) is that the latter includes interactions between the first three and last two variables while the former excludes them.

The great convenience of Sobol’s measures is that they can be directly estimated by integration, without explicitly estimating all of the necessary interaction effects, squaring them, integrating their squares and summing those integrated squared estimates. Sobol’ provides a kind of tomography: integrals of cross-products of $f$ reveal facts about the internal structure of $f$.

The goal of this paper is to exhibit the entire space of linear combinations of cross-moments of function evaluations with some variables fixed and others independently sampled. Such a linear combination is a generalized Sobol’ index, or GSI. Then, using this space of functions, we make a systematic search for estimators of interpretable quantities with desirable computational or statistical properties.

This systematic approach yields some new and useful estimators. Some have reduced cost compared to previously known ones. Some have reduced sampling variance. It also encompasses some earlier work. In particular, an efficient strategy to estimate all two factor interaction mean squares due to Saltelli, (2002) appears as a special case.

Section 2 introduces some notation and reviews the ANOVA of $^{d}$ and Sobol’ indices. A compact notation is necessary to avoid cumbersome expressions with many indices. Section 3 defines the generalized Sobol’ indices and gives an expression for their value. It also defines several special classes of GSI based on interpretability, computational efficiency, or statistical considerations. These are contrasts, squares, sums of squares and bilinear GSIs. Section 4 shows that many GSIs including the Sobol’ index (1) cannot be estimated by unbiased sums of squares. Section 5 considers estimation of a specific variance component for a proper subset containing $k$ of the variables. A direct approach requires $2^{k}$ function evaluations per Monte Carlo sample, while a bilinear estimate reduces the cost to $2^{\lfloor k/2\rfloor}+2^{k-\lfloor k/2\rfloor}$. That section also introduces a bilinear estimate for the superset importance measure defined in Section 2. Section 6 considers some GSIs for high dimensional problems. It includes a contrast GSI which estimates the mean square dimension of a function of $d$ variables using only $d+1$ function evaluations per Monte Carlo trial as well as some estimators of the mean dimension in the truncation sense. Section 7 presents a bias correction for GSIs that are not contrasts. Section 8 makes some comparisons among alternative methods and Section 9 has conclusions.

Background and notation

The analysis of variance originates with Fisher and Mackenzie, (1923). It partitions the variance of a quantity among all non-empty subsets of factors, defined on a finite Cartesian grid.

The ANOVA was generalized by Hoeffding, (1948) to functions in L^{2}^{d} for integer $d\geqslant 1$. That generalization extends the one for factorial experimental designs in a natural way, and can be applied to $L^{2}$ functions on any tensor product domain. For $d=\infty$, see Owen, (1998).

The ANOVA of L^{2}^{d} is also attributed to Sobol’, (1969). For historical interest, we note that Sobol’ used a different approach than Hoeffding. He represented $f$ by an expansion in a complete orthonormal basis (tensor products of Haar functions) and gathered together terms corresponding to each subset of variables. That is, where Hoeffding has an analysis, Sobol’ has a synthesis.

We use $\boldsymbol{x}=(x_{1},x_{2},\dots,x_{d})$ to represent a typical point in $^{d}$. The set of indices is ${\cal D}=\{1,2,\dots,d\}$. We write $u\subset v$ to denote a proper subset, that is $u\subsetneq v$. For $u\subseteq{\cal D}$ we use $|u|$ to denote the cardinality of $u$, and either $-u$ or $u^{c}$ (depending on typographical clarity) to represent the complementary set ${\cal D}-u$. The expression $u+v$ means $u\cup v$ where $u$ and $v$ are understood to be disjoint.

The ANOVA decomposition represents $f(\boldsymbol{x})$ via

where the functions $f_{u}$ are defined recursively by

In statistical language, $u$ is a set of factors and $f_{u}$ is the corresponding effect. We get $f_{u}$ by subtracting sub-effects of $f_{u}$ from $f$ and then averaging the residual over $x_{j}$ for $j\not\in u$.

This section introduces the Sobol’ indices that we generalize, and mentions some of the methods for their estimation.

There are various ways that one might define the importance of a variable $x_{j}$. The importance of variable $j\in\{1,\dots,d\}$ is due in part to $\sigma^{2}_{\{j\}}$, but also due to $\sigma^{2}_{u}$ for other sets $u$ with $j\in u$. More generally, we may be interested in the importance of $\boldsymbol{x}_{u}$ for a subset $u$ of the variables.

Sobol’, (1993) introduced two measures of variable subset importance, which we denote

The upper index counts every ANOVA component that touches the set $u$ in any way. If $\underline{\tau}_{u}^{2}$ is large then the subset $u$ is clearly important. If $\overline{\tau}_{u}^{2}$ is small, then the subset $u$ is not important and Sobol’ et al., (2007) investigate the effects of fixing such $\boldsymbol{x}_{u}$ at some specific values. The second measure includes interactions between $\boldsymbol{x}_{u}$ and $\boldsymbol{x}_{-u}$ while the first measure does not.

These Sobol’ indices satisfy $\underline{\tau}_{u}^{2}\leqslant\overline{\tau}_{u}^{2}$ and $\underline{\tau}_{u}^{2}+\overline{\tau}_{-u}^{2}=\sigma^{2}$. Sobol’ usually normalizes these quantities, yielding global sensitivity indices $\underline{\tau}_{u}^{2}/\sigma^{2}$ and $\overline{\tau}_{u}^{2}/\sigma^{2}$. In this paper we work mostly with unnormalized versions.

Sobol’s original work was published in Sobol’, (1990) before being translated in Sobol’, (1993). Ishigami and Homma, (1990) independently considered computation of $\underline{\tau}^{2}_{\{j\}}$.

To estimate Sobol’ indices, one pairs the point $\boldsymbol{x}$ with a hybrid point $\boldsymbol{y}$ that shares some but not all of the components of $\boldsymbol{x}$. We denote the hybrid point by $\boldsymbol{y}=\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u}$ where $y_{j}=x_{j}$ for $j\in u$ and $y_{u}=z_{j}$ for $j\not\in u$.

From the ANOVA properties one can show directly that

Mauntz, (2002) and Kucherenko et al., (2011) use an estimator for $\underline{\tau}^{2}_{u}$ derived as a sample version of the identity

There are $2^{d}-1$ variance components $\sigma^{2}_{u}$ as well as $2^{d}-1$ Sobol’ indices $\underline{\tau}^{2}_{u}$ and $\overline{\tau}^{2}_{u}$ of each type. We can recover any desired $\sigma^{2}_{u}$ as a linear combination of $\underline{\tau}^{2}_{v}$. For example $\sigma^{2}_{\{1,2,3\}}=\underline{\tau}^{2}_{\{1,2,3\}}-\underline{\tau}^{2}_{\{1,2\}}-\underline{\tau}^{2}_{\{1,3\}}-\underline{\tau}^{2}_{\{2,3\}}+\underline{\tau}^{2}_{\{1\}}+\underline{\tau}^{2}_{\{2\}}+\underline{\tau}^{2}_{\{3\}}$. More generally, we have the Moebius-type relation

Because $f$ is defined on a unit cube and can be computed at any desired point, methods other than simple Monte Carlo can be applied. Quasi-Monte Carlo (QMC) sampling (see Niederreiter, (1992)) can be used instead of plain Monte Carlo, and Sobol’, (2001) reports that QMC is more effective. For functions $f$ that are very expensive, a Bayesian numerical analysis approach (Oakley and O’Hagan,, 2004) based on a Gaussian process model for $f$ is an attractive way to compute Sobol’ indices.

2 Related indices

Another measure of the importance of $\boldsymbol{x}_{u}$ is the superset importance measure

used by Hooker, (2004) to quantify the effect of dropping all interactions containing the set $u$ of variables from a black box function.

Sums of ANOVA components are also used in quasi-Monte Carlo sampling. QMC is, in general, more effective on integrands $f$ that are dominated by their low order ANOVA components. Two versions of $f$ that are equivalent in Monte Carlo sampling may behave quite differently in QMC. For example the basis used to sample Brownian paths has been seen to affect the accuracy of QMC integrals (Caflisch et al.,, 1997; Acworth et al.,, 1997; Imai and Tan,, 2002).

The function $f$ has effective dimension $s$ in the superposition sense (Caflisch et al.,, 1997), if $\sum_{|u|\leqslant s}\sigma^{2}_{u}\geqslant(1-\epsilon)\sigma^{2}$. Typically $\epsilon=0.01$ is used as a default. Similarly, $f$ has effective dimension $s$ in the truncation sense (Caflisch et al.,, 1997), if $\sum_{u\subseteq\{1,2,\dots,s\}}\sigma^{2}_{u}=\underline{\tau}^{2}_{\{1,2,\dots,s\}}\geqslant(1-\epsilon)\sigma^{2}$.

It is much easier to estimate the mean dimension (superposition sense) defined as $\sum_{u}|u|\sigma^{2}_{u}/\sigma^{2}$ than the effective dimension, because the mean dimension is a linear combination of variance components. The mean dimension also offers better resolution than effective dimension. For instance, two functions having identical effective dimension $2$ might have different mean dimensions, say $1.03$ and $1.05$. Similarly one can estimate a mean square dimension $\sum_{u}|u|^{2}\sigma^{2}_{u}/\sigma^{2}$:

This follows from Theorem 2 of Liu and Owen, (2006). ∎

Generalized Sobol’ indices

Here we consider a general family of quadratic indices similar to those of Sobol’. The general form of these indices is

for coefficients $\Omega_{uv}$. If we think of $f(\boldsymbol{x})$ as being the standard evaluation, then the generalized Sobol’ indices (9) are linear combinations of all possible second order moments of $f$ based on fixing two subsets, $u$ and $v$, of input variables.

The Sobol’ matrix is symmetric. It also satisfies $\Theta_{uv}=\Theta_{u^{c}v^{c}}$. Theorem 2 gives the general form of a Sobol’ matrix entry.

Equation (9) describes a $2^{2d}$ dimensional family of linear combinations of pairwise function products. There are only $2^{d}-1$ ANOVA components to estimate. Accordingly we are interested in special cases of (9) with desirable properties.

A generalized Sobol’ index is a square if it takes the form

Squares and sums of squares have the advantage that they are non-negative and hence avoid the problems associated with negative sample variance components. A square GSI can be written compactly as ${\text{tr}}(\lambda\lambda^{\mathsf{T}}\Theta)=\lambda^{\mathsf{T}}\Theta\lambda$ where $\lambda$ is a vector of $2^{d}$ coefficients. If $\lambda_{u}$ is sparse (mostly zeros) then a square index is inexpensive to compute.

A generalized Sobol’ index is bilinear if it takes the form

Bilinear estimates have the advantage of being rapidly computable. If there are $\|\lambda\|_{0}$ nonzero elements in $\lambda$ and $\|\gamma\|_{0}$ nonzero elements in $\gamma$ then the integrand in a bilinear generalized Sobol’ index can be computed with at most $\|\gamma\|_{0}+\|\lambda\|_{0}$ function calls and sometimes fewer (see Section 3.3) even though it combines values from $\|\gamma\|_{0}\times\|\lambda\|_{0}$ function pairs. We can write the bilinear GSI as ${\text{tr}}(\lambda\gamma^{\mathsf{T}}\Theta)=\gamma^{\mathsf{T}}\Theta\lambda$. The sum of a small number of bilinear GSIs is a low rank GSI.

A GSI is simple if it is written as a linear combination of entries in just one row or just one column of $\Theta$, such as

It is convenient if the chosen row or column corresponds to $u$ or $v$ equal to $\varnothing$ or ${\cal D}$. Any linear combination $\sum_{u}\delta_{u}(\mu^{2}+\underline{\tau}^{2}_{u})$ of variance components and $\mu^{2}$ can be written as a simple GSI taking $\lambda_{u}=\delta_{-u}$. There are computational advantages to some non-simple representations.

2 Sample GSIs

We can derive a matrix expression for the estimator by introducing the vectors

The vectors $F_{i}$ have covariance $\Theta-\mu^{2}$. Then

3 Cost per (𝒙,𝒛)𝒙𝒛(\boldsymbol{x},\boldsymbol{z}) pair

We suppose that the cost of computing a sample GSI is dominated by the number of function evaluations required. If the GSI requires $C(\Omega)$ (defined below) distinct function evaluations per pair $(\boldsymbol{x}_{i},\boldsymbol{z}_{i})$ for $i=1,\dots,n$, then the cost of the sample GSI is proportional to $nC(\Omega)$.

If the row $\Omega_{uv}$ for given $u$ and all values of $v$ is not entirely zero then we need the value $f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})$. Let

indicate whether $f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})$ is needed as the ‘left side’ of a product $f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})f(\boldsymbol{x}_{v}{:}\boldsymbol{z}_{-v})$. The number of function evaluations required for the GSI ${\text{tr}}(\Omega^{\mathsf{T}}\Theta)$ is:

We count the number of rows of $\Omega$ for which $f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})$ is needed, add the number of columns and then subtract the number of double counted sets $u$.

Squares and sums of squares

A square or sum of squares yields a nonnegative estimate. An unbiased and nonnegative estimate is especially valuable. When the true GSI is zero, an unbiased nonnegative estimate will always return exactly zero as Fruth et al., (2012) remark. The Sobol’ index $\overline{\tau}^{2}_{u}$ is of square form, but $\underline{\tau}^{2}_{u}$ is not. Theorem 1 leads to a sum of squares for $\sum_{u}|u|\sigma^{2}_{u}$.

Liu and Owen, (2006) express the superset importance measure as a square:

Fruth et al., (2012) find that a sample version of (11) is the best among four estimators of $\Upsilon^{2}_{u}$.

In classical crossed mixed effects models (Montgomery,, 1998) every ANOVA expected mean square has a contribution from the highest order variance component, typically containing measurement error. A similar phenomenon applies for Sobol’ indices. In particular, no square GSI or sum of squares will yield $\underline{\tau}^{2}_{u}$ for $|u|<d$, as the next proposition shows.

It is enough to show this for $R=1$ with $\lambda_{1,u}=\lambda_{u}$. First,

Next, the only elements of $\Theta$ containing $\sigma^{2}_{{\cal D}}$ are the diagonal ones, equal to $\sigma^{2}$. Therefore the coefficient of $\sigma_{{\cal D}}^{2}$ in (12) is $\sum_{u}\lambda^{2}_{u}$. ∎

For a square or sum of squares to be free of $\sigma^{2}_{{\cal D}}$ it is necessary to have $\sum_{r=1}^{R}\sum_{u}\lambda_{r,u}^{2}=0$. That in turn requires all the $\lambda_{r,u}$ to vanish, leading to the degenerate case $\Omega=0$. As a result, we cannot get an unbiased sum of squares for any GSI that does not include $\sigma^{2}_{{\cal D}}$. In particular, $\underline{\tau}^{2}_{u}$ cannot have an unbiased sum of squares estimate for $|u|<d$.

Specific variance components

For any $w\subseteq{\cal D}$ the variance component for $w$ is given in (5) as an alternating sum of $2^{|w|}$ lower Sobol’ indices. It can thus be estimated by a simple GSI,

where $\lambda_{v}=(-1)^{|w-v|}$. The cost of this simple GSI is $C=2^{|w|}+1_{|w|<d}$. If $w={\cal D}$, then $f(\boldsymbol{x})$ appears twice, but otherwise it is only used once. The GSI can also be estimated by some bilinear GSIs using fewer function evaluations as we show here.

We begin by noting that for $u,v\subseteq w$,

To illustrate, suppose that $w=\{1,2,3\}$. Let $u,v\subseteq w$. Then we can work out a $2\times 8$ submatrix of the Sobol’ matrix using

where we omit braces and commas from the set notation.

It follows now that we can use a non-simple bilinear GSI, $\lambda^{\mathsf{T}}\Theta\gamma$, where $\lambda$ and $\gamma$ are given by

with $u$ and $v-w^{c}$ given along the top labels in (17) while the remaining $2^{d}-8$ elements of $\lambda$ and $\gamma$ are all zero. Specifically, the expected value of

is $\sigma^{2}_{\{1,2,3\}}$. While equation (13) for $w=\{1,2,3\}$ requires $9$ function evaluations per $(\boldsymbol{x},\boldsymbol{z})$ pair, equation (18) only requires $6$ function evaluations. For $|w|<d$, the $u=\varnothing$ and $v=\varnothing$ evaluations are different due to the presence of $w^{c}$, so no evaluations are common to both the $\lambda$ and $\gamma$ expressions. There are also two variants of (17) that single out variables $2$ and $3$ respectively, analogously to the way that (18) treats variable $1$.

In general, bilinear GSIs let us estimate $\sigma^{2}_{w}$ using $2^{k}+2^{|w|-k}$ function evaluations per $(\boldsymbol{x},\boldsymbol{z})$ pair for integer $1\leqslant k<|w|$ instead of the $2^{|w|}$ evaluations that a simple GSI requires.

Let $w$ be a nonempty subset of ${\cal D}$ for $d\geqslant 1$. Let f\in L^{2}^{d}. Choose $w_{1}\subseteq w$ and put $w_{2}=w-w_{1}$. Then

after a change of variable from $u_{j}$ to $w_{j}-u_{j}$ for $j=1,2$. We may write the above as

Consider the set $v\subseteq{\cal D}$. The coefficient of $\sigma^{2}_{v}$ in (20) is if $v\cap w^{c}\neq\varnothing$. Otherwise, we may write $v=v_{1}+v_{2}$ where $v_{j}\subseteq w_{j}$, $j=1,2$. Then the coefficient of $\sigma^{2}_{v}$ in (20) is

These alternating sums over $u_{j}$ with $v_{j}\subseteq u_{j}\subseteq w_{j}$ equal $1$ if $v_{j}=w_{j}$ but otherwise they are zero. Therefore the coefficient of $\sigma^{2}_{v}$ in (20) is $1$ if $v=w$ and is otherwise. ∎

We can use Theorem 3 to get a bilinear (but not square) estimator of $\sigma^{2}_{{\cal D}}=\Upsilon^{2}_{{\cal D}}$. A similar argument to that in Theorem 3 yields a bilinear estimator of superset importance $\Upsilon^{2}_{w}$ for a general set $w$.

Let $w$ be a nonempty subset of ${\cal D}$ for $d\geqslant 1$. Let f\in L^{2}^{d}. Choose $w_{1}\subseteq w$ and put $w_{2}=w-w_{1}$. Then

Now write $v=(v\cap w^{c})+v_{1}+v_{2}$ with $v_{j}\subseteq w_{j}$, $j=1,2$. The coefficient of $\sigma^{2}_{v}$ is

which vanishes unless $w_{1}=v_{1}$ and otherwise equals $1$. Therefore the coefficient of $\sigma^{2}_{v}$ is $1$ if $v\supseteq w$ and is otherwise. ∎

The cost of the estimator (21) is $C=2^{|w_{1}|}+2^{|w_{2}|}-1$, because the evaluation $f(\boldsymbol{x}_{w^{c}}{:}\boldsymbol{z}_{w})$ can be used for both $u_{1}=\varnothing$ and $u_{2}=\varnothing$.

GSIs with O​(d)𝑂𝑑O(d) function evaluations per pair

Some problems, like computing mean dimension, can be solved with $O(d)$ different integrals instead of the $O(2^{d})$ required to estimate all ANOVA components. In this section we enumerate what can be estimated by certain GSIs based on only $O(d)$ carefully chosen function evaluations per $(\boldsymbol{x}_{i},\boldsymbol{z}_{i})$ pair.

and hence the accessible elements of the Sobol’ matrix are:

Using (23) we can construct estimates of $\sum_{u}|u|\sigma^{2}_{u}=\sum_{j}\overline{\tau}^{2}_{\{j\}}$ and $\sum_{|u|=1}\sigma^{2}_{u}=\sum_{j}\underline{\tau}^{2}_{\{j\}}$ at cost $C=d+1$. Simple GSI estimates are available using $u=\varnothing$ and either $v=\{j\}$ or $v=-\{j\}$ for $j=1,\dots,d$. More interestingly, it is possible to compute all $d(d-1)/2$ indices $\overline{\tau}^{2}_{\{j,k\}}$ along with all $\underline{\tau}^{2}_{\{j\}}$ and $\overline{\tau}^{2}_{\{j\}}$ for $j=1,\dots,d$, at total cost $C=d+2$ as was first shown by Saltelli, (2002, Theorem 1). Given $C=2d+2$ evaluations one can also compute all of the $\underline{\tau}^{2}_{\{j,k\}}$ indices by pairing up $u=\{j\}$ and $v=-\{k\}$ (Saltelli,, 2002, Theorem 2).

For the remainder of this section we present some contrast estimators. The estimate

is both a contrast and a sum of squares. It has expected value $\sum_{u}|u|\sigma^{2}_{u}$ and cost $C=d+1$.

Next, to estimate $\sum_{u}1_{|u|=1}\sigma^{2}_{u}$ by a contrast using $d+2$ function evaluations per $(\boldsymbol{x}_{i},\boldsymbol{z}_{i})$ pair, let

Then the contrast $\sum_{u}\lambda_{u}f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})f(\boldsymbol{z})$ has expected value

The total of second order interactions can be estimated with a contrast at cost $C=2d+2$. Taking

Thus $\lambda^{\mathsf{T}}\widehat{\Theta}\gamma/2$ estimates $\sum_{|u|=2}\sigma^{2}_{u}$, at cost $C=2d+2$.

2 Consecutive index GSIs

A second way to reduce function evaluations to $O(d)$ is to consider only subsets $u$ and $v$ of the form $\{1,2,\dots,j\}$ and $\{j+1,\dots,d\}$. We write these as $(0,j]$ and $(j,d]$ respectively. If $f(\boldsymbol{x})$ is the result of a process evolving in discrete time then $(0,j]$ represents the effects of inputs up to time $j$ and $(j,d]$ represents those after time $j$. A small value of $\overline{\tau}^{2}_{(0,j]}$ then means that the first $j$ inputs are nearly forgotten while a large value for $\underline{\tau}^{2}_{(0,j]}$ means the initial conditions have a lasting effect.

and hence the accessible elements of the Sobol’ matrix are $\mu^{2}+\underline{\tau}^{2}_{u}$ for the sets $u$ in the array above.

The same strategies used on singletons and their complements can be applied to consecutive indices. They yield interesting quantities related to mean dimension in the truncation sense. To describe them, we write $\lfloor u\rfloor=\min\{j\mid j\in u\}$ and $\lceil u\rceil=\max\{j\mid j\in u\}$ for the least and greatest indices in the non-empty set $u$.

Let f\in L^{2}^{d} have variance components $\sigma^{2}_{u}$. Then

Since these are contrasts, we may suppose that $\mu=0$. Then, using (25)

Next, $\Theta_{\varnothing,{\cal D}}=0$, and

Combining these yields the first result. The second is similar. ∎

Using Proposition 2 we can obtain an estimate of $\sum_{u}\lceil u\rceil\sigma^{2}_{u}/\sigma^{2}$, the mean dimension of $f$ in the truncation sense. We also obtain a contrast

which measures the extent to which indices at distant time lags contribute important interactions.

We can also construct GSIs based on pairs of segments. For example,

Bias corrected GSIs

When we are interested in estimating a linear combination of variance components, then the corresponding GSI is a contrast. Sometimes estimating a contrast requires an additional function evaluation per $(\boldsymbol{x}_{i},\boldsymbol{z}_{i})$ pair. For instance the unbiased estimator (4) of $\underline{\tau}^{2}_{u}$ requires three function evaluations per pair compared to the two required by the biased estimator of Janon et al., (2012).

Proposition 3 supplies a bias-corrected version of Janon et al.’s (2011) estimator of $\underline{\tau}^{2}_{u}$ using only two function evaluations per $(\boldsymbol{x}_{i},\boldsymbol{z}_{i})$ pair.

is an unbiased estimate of $\sum_{u}\sum_{v}\Omega_{uv}(\Theta_{uv}-\mu^{2})$.

This follows by applying Proposition 3 term by term. ∎

The computational burden for the unbiased estimator in Proposition 4 is not much greater than that for the possibly biased estimator ${\text{tr}}(\Omega^{\mathsf{T}}\widehat{\Theta})$. It requires no additional function evaluations. The quantities $\hat{\mu}_{u}$ and $s_{u}^{2}$ need only be computed for sets $u\subseteq{\cal D}$ for which $\Omega_{uv}$ or $\Omega_{vu}$ is nonzero for some $v$. If $\Omega$ is a sum of bilinear estimators then the $-\sum_{u}\sum_{v}\Omega_{uv}\hat{\mu}_{u}\hat{\mu}_{v}/2$ cross terms also have that property.

The bias correction in estimator (26) complicates calculation of confidence intervals for ${\text{tr}}(\Omega^{\mathsf{T}}\Theta)$. Jackknife or bootstrap methods will work but confidence intervals for contrasts are much simpler because the estimators are simple averages.

Comparisons

There is a $2^{2d}$–dimensional space of GSIs but only a $2^{d}-1$–dimensional space of linear combinations of variance components to estimate. As a result there is more than one way to estimate a desired linear combination of variance components.

As a case in point the Sobol’ index $\underline{\tau}^{2}_{u}$ can be estimated by either the original method or by the contrast (4). Janon et al., (2012) prove that their estimate of $\hat{\mu}$ improves on the simpler one and establish asymptotic efficiency for their estimator within a class of methods based on exchangeability, but that class does not include the contrast. Similarly, inspecting the Sobol’ matrix yields at least four ways to estimate the variance component $\sigma^{2}_{\{1,2,3\}}$, and superset importance can be estimated via a square or a bilinear term.

Here we consider some theoretical aspects of the comparison, but they do not lead to unambiguous choices. Next we consider a small set of empirical investigations.

Ideally we would like to choose $\Omega$ to minimize the variance of the sample GSI. But, the variance of a GSI depends on fourth moments of ANOVA contributions which are ordinarily unknown and harder to estimate than the variance components themselves.

The same issue comes up in the estimation of variance components, where MINQE (minimum norm quadratic estimation) estimators were proposed in a series of papers by C. R. Rao in the 1970s. For a comprehensive treatment see Rao and Kleffe, (1988) who present MINQUE and MINQIE versions using unbiasedness or invariance as constraints. The idea in MINQUE estimation is to minimize a convenient quadratic norm as a proxy for the variance of the estimator.

The GSI context involves variance components for crossed random effects models with interactions of all orders. Even the two way crossed random effects model with an interaction is complicated enough that no closed form estimator appears to be known for that case. See Kleffe, (1980).

We can however generalize the idea behind MINQE estimators to the GSI setting. Writing

Using the proxy for variance suggests choosing the estimator which minimizes $C(\Omega)\times V(\Omega)$. The contrast estimator (4) of $\underline{\tau}^{2}_{u}$ has $C(\Omega)\times V(\Omega)=3\times 2=6$ while the original Sobol’ estimator has $C(\Omega)\times V(\Omega)=2\times 1=2$. The estimators (13) and (18) for $\sigma^{2}_{\{1,2,3\}}$ both have $V(\Omega)=8$. The former has cost $C(\Omega)=9$, while the latter costs $C(\Omega)=6$. As a result, the proxy arguments support the original Sobol’ estimator and the alternative estimator (18) for $\sigma^{2}_{\{1,2,3\}}$.

2 Test cases

To compare some estimators we use test functions of product form:

The third condition ensures that all GSIs have finite variance, while the first two allow us to write the variance components of $f$ as

We will compare Monte Carlo estimates and so smoothness or otherwise of $g_{j}(\cdot)$ play no role. Only $\mu_{j}$, $\tau_{j}$ and the third and fourth moments of $g$ play a role. Monte Carlo estimation is suitable when $f$ is inexpensive to evaluate, like surrogate functions in computer experiments. For our examples we take $g_{j}(x)=\sqrt{12}(x-1/2)$ for all $j$.

For an example function of non-product form, we take the minimum,

for this function. Taking $u={\cal D}$, gives $\sigma^{2}=d(d+1)^{-2}(d+2)^{-1}$.

We considered both simple and bilinear estimators of $\sigma^{2}_{\{1,2,3\}}$ in Section 5. The simple estimator requires $9$ function evaluations per $(\boldsymbol{x},\boldsymbol{z})$ pair, while three different bilinear ones each require only $6$.

For a function of product form, all four of these estimators yield the same answer for any specific set of $(\boldsymbol{x}_{i},\boldsymbol{z}_{i})$ pairs. As a result the bilinear formulas dominate the simple one for product functions.

For the minimum function, with $d=5$ we find that by symmetry,

Because we are interested in comparing the variance of estimators of a variance, a larger sample is warranted than if we were simply estimating a variance component. Based on $1{,}000{,}000$ function evaluations we find the estimated means and standard errors are given in Table 1. We see that the bilinear estimators give about half the standard error of the simple estimator, corresponding to about $(1.05/.571)^{2}\times 9/6\doteq 5.1$ times the statistical efficiency.

We consider two estimators of $\underline{\tau}^{2}_{u}$. The estimator (4) is a bilinear contrast, averaging $f(\boldsymbol{x})(f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})-f(\boldsymbol{z}))$. The estimator (3) using the estimator of $\hat{\mu}$ from Janon et al., (2012) is a modification of Sobol’s original simple estimator based on averaging $f(\boldsymbol{x})f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})$. The bias correction of Section 7 makes an asymptotically negligible difference, so we do not consider it here.

The contrast estimator requires $3$ function evaluations per $(\boldsymbol{x},\boldsymbol{z})$ pair, while Sobol’s only requires $2$. Both estimators make an adjustment to compensate for the bias $\mu^{2}$. Estimator (3) subtracts an estimate $\hat{\mu}^{2}$ based on combining all $2n$ function evaluations, the square of the most natural way to estimate $\mu$ from the available data. Estimator (4) subtracts $(1/n^{2})\sum_{i}\sum_{i^{\prime}}f(\boldsymbol{x}_{i})f(\boldsymbol{x}_{i,u}{:}\boldsymbol{z}_{i,-u})$, which may be advantageous when the difference $f(\boldsymbol{x}_{u}{:}\boldsymbol{z}_{-u})-f(\boldsymbol{z})$ involves considerable cancellation, as it might if $\boldsymbol{x}_{u}$ is unimportant. Thus we might expect (4) to be better when $\underline{\tau}^{2}_{u}$ is small. We compare the estimators on a product function, looking at $\underline{\tau}^{2}_{u}$ for three subsets $u$ of size $2$ and varying importance.

For the product function with $d=6$, $\tau=(1,1,1/2,1/2,1/4,1/4)$, all $\mu_{j}=1$ for $j=1,\dots,6$, and $g_{j}(x_{j})=\sqrt{12}(x_{j}-1/2)$, we may compute $\underline{\tau}^{2}_{\{1,2\}}=3\doteq 0.50\sigma^{2}$, $\underline{\tau}^{2}_{\{3,4\}}\doteq 1.56\doteq 0.093\sigma^{2}$, and $\underline{\tau}^{2}_{\{5,6\}}\doteq 0.13\doteq 0.021\sigma^{2}$.

Results from $R=10{,}000$ trials with $n=10{,}000$ $(\boldsymbol{x}_{i},\boldsymbol{z}_{i})$ pairs each, are shown in Table 2. The efficiency of the contrast estimator compared to the simple one ranges from about $0.5$ to about $2.5$ in this example, depending on the size of the effect being estimated, with the contrast being better for the small quantity $\underline{\tau}^{2}_{\{5,6\}}$. Sobol’ et al., (2007) also report superiority of the contrast estimator on a small $\underline{\tau}^{2}_{u}$.

Neither estimator is always more efficient than the other, hence no proxy based solely on $\Omega$ can reliably predict which of these is better for a specific problem.

The bias correction from Section 7 makes little difference here because for $n=10{,}000$ there is very little bias to correct. It does make a difference when $n=100$ (data not shown) but at such small sample sizes the standard deviation of $\widehat{}\underline{\tau}^{2}_{u}$ can be comparable to or larger than $\underline{\tau}^{2}_{u}$ itself for this function.

Here we compare two estimates of $\Upsilon^{2}_{\{1,2,3,4\}}$, the square (11) and the bilinear estimator (21) from Theorem 4. For a product function, $\Upsilon^{2}_{w}=\prod_{j\in w}\tau^{2}_{j}\prod_{j\not\in w}(\mu_{j}^{2}+\tau^{2}_{j}).$ Squares have an advantage estimating small GSIs so we consider one small and one large (for a four way interaction) $\Upsilon^{2}$.

For $d=8$, $\tau=c(4,4,3,3,2,2,1,1)/4$ and all $\mu_{j}=1$ we find that $\Upsilon^{2}_{\{1,2,3,4\}}\doteq 0.558\doteq 0.0334\sigma^{2}$ and $\Upsilon^{2}_{\{5,6,7,8\}}\doteq 0.00238\doteq 0.000147\sigma^{2}$. The bilinear estimate (21) based on $w_{1}=\{1,2\}$ and $w_{2}=\{3,4\}$ for $\Upsilon_{\{1,2,3,4\}}$ (respectively $w_{1}=\{5,6\}$ and $w_{2}=\{7,8\}$ for $\Upsilon_{\{5,6,7,8\}}$) requires $C=7$ function evaluations, while the square (11) requires $C=16$. From Table 3 we see that the square has an advantage that more than compensates for using a larger number of function evaluations and the advantage is overwhelming for the smaller effect.

The outlook for the bilinear estimator of $\Upsilon^{2}_{w}$ is pessimistic. Its cost advantage grows with $|w|$; for $|w|=20$ it has cost $1023$ compared to $2^{20}$ for the square. But $\Upsilon^{2}_{w}$ for such a large $w$ will often be so small that the variance advantage from using a square will be extreme.

Conclusions

We have generalized Sobol’ indices to estimators of arbitrary linear combinations of variance components. Sometimes there are multiple ways to estimate a generalized Sobol’ index with important efficiency differences. Square GSIs where available are very effective. When no square or sum of squares is available a bilinear or low rank GSI can at least save some function evaluations. Contrasts are simpler to work than other GSIs, because they avoid bias corrections.

Acknowledgments

This work was supported by the U.S. National Science Foundation under grant DMS-0906056. I thank Alexandra Chouldechova for translating Sobol’s description of the analysis of variance. Thanks to Sergei Kucherenko for discussions on Sobol’ indices. I also thank the researchers of the GDR MASCOT NUM for an invitation to their 2012 meeting which lead to the research presented here.

References