testHSIC: Tests of Independence based on the Hilbert-Schmidt...

In sensitivity: Global Sensitivity Analysis of Model Outputs and Importance Measures

Tests of Independence based on the Hilbert-Schmidt Independence Criterion (HSIC)

Description

testHSIC allows to test independence among all input-output pairs (X_i, Y) after a preliminary sensitivity analysis based on HSIC indices. testHSIC takes an object of class sensiHSIC (produced by a prior call to the function sensiHSIC that estimates HSIC indices) and it returns the estimated p-values after testing independence among all input-output pairs. For each input-output pair, having access to the p-value helps the user decide whether the null hypothesis H_0: "X_i and Y are independent" must be accepted or rejected. If all kernels selected in sensiHSIC are characteristic, H_0 can be rewritten "HSIC(X_i, Y)=0" and this paves the way to several test procedures.

Depending on the sample size, the selected HSIC-based sensitivity measure and the estimator type (either a U-statistic or a V-statistic), there are up to four different methods to test H_0. The asymptotic test is recommended when the sample size n is around a few hundreds (or more). When n is smaller, a permutation-based test must be considered instead. As a general rule, permutation-based tests can always be applied but a much heavier computational load is to be expected. However, if HSIC indices were initially estimated with V-statistics, the Gamma test is a parametric method that offers an enticing tradeoff.

If the total-order HSIC-ANOVA indices were estimated during the preliminary sensitivity analysis, their numerators T_i also allow to test independence among all input-output pairs (X_i, Y). For a given input variable X_i, let T_i denote the sum of HSIC(X_i, Y) and all HSIC-ANOVA interaction effects between X_i and other inputs. If all kernels selected in sensiHSIC are characteristic, the null hypothesis H_0: "X_i and Y are independent" can be rewritten H_0: "T_i=0". The resulting test procedures prove to be more powerful, which means the probability of truly detecting dependence among (X_i, Y) is increased.

The test statistic associated to T_i can be either a U-statistic or a V-statistic. Several kinds of test procedures can then be considered. They are all inspired from what is done for HSIC indices. A permutation-based test can always been applied but it becomes too expensive (from a computational viewpoint) when the sample size n is large. A much cheapier strategy consists in using the Gamma test but it is limited to the specific case where T_i is estimated with a V-statistic.

Usage

testHSIC(sensi, test.method = "Asymptotic", B = 3000,
         seq.options = list(criterion = "screening", alpha = 0.05,
                            Bstart = 200, Bfinal = 5000, Bbatch = 100, Bconv = 200,
                            graph = TRUE) )

## S3 method for class 'testHSIC'
print(x, ...)

## S3 method for class 'testHSIC'
plot(x, ylim = c(0, 1), err, ...)

Arguments

sensi	An object of class "sensiHSIC" which is produced by a prior call to the function sensiHSIC. In particular, sensi must contain objects named "KX" (3D-array filled with all input Gram matrices), "KY" (output Gram matrix), "HSICXY" (estimated HSIC indices) and "estimator.type" (either "U-stat" or "V-stat"). If the intended test procedure is based on total-order HSIC-ANOVA indices, sensi must also contain an object named "TO.num" (numerators of the estimated total-order HSIC-ANOVA indices). If the intended test procedure is based on C-HSIC indices, sensi must also contain objects named "cond" (list of options describing the conditioning event) and "weights" (conditioning weights after normalization).
test.method	A string specifying the numerical procedure used to estimate the p-values. For a test of independence based on HSIC indices, available procedures include "Asymptotic" (asymptotic test), "Permutation" (permutation-based test), "Seq_Permutation" (sequential permutation-based test) and "Gamma" (Gamma test). If sensi contains V-statistics, the asymptotic test (resp. the Gamma test) is recommended for large (resp. small) sample sizes. Otherwise, permutation-based tests can be used as well. If sensi contains U-statistics, the Gamma test must not be employed. The asymptotic test is recommended for large sample sizes. Otherwise, permutation-based tests can be used as well. For a test of independence based on total-order HSIC-ANOVA indices, available procedures include "Tot_Permutation" (permutation-based test), "Tot_Seq_Permutation" (sequential permutation-based test) and "Tot_Gamma" (Gamma test). If sensi contains V-statistics, the Gamma test is recommended. Otherwise, permutation-based tests can be used as well. If sensi contains U-statistics, only permutation-based tests can be used.
B	Number of random permutations carried out on the output samples before the non-parametric estimation of p-values. Only relevant if test.method is "Permutation" or "Tot_Permutation".
seq.options	A list of options guiding the sequential procedure. Only relevant if test.method is "Seq_Permutation" or "Tot_Seq_Permutation". criterion is a string specifying the stopping criterion. Available criteria include "screening" (permutations stop as soons as the estimated p-values have sufficiently converged so that they can be compared to the reference threshold alpha), "ranking" (permutations stop as soon as the estimated p-values have sufficiently converged so that they can be ranked) and "both" (permutations stop as soon as the two previous criteria are fulfilled). alpha is a scalar value (between 0 and 1) specifying the type I error (probability of wrongly accepting H_0). Only relevant if criterion is "screening" or "both". Bstart is the initial number of random permutations before the first criterion check. Bfinal is the maximum number of random permutations. Bbatch is the number of permutations at each new iteration of the sequential procedure. Bconv is the number of permutations that is used to determine whether convergence has already occured or not. For criterion="screening", convergence is assumed to be reached if the positions of the estimated p-values with respect to alpha no longer evolve after the Bconv latest permutations. For criterion="ranking", convergence is assumed to be reached if the rankings of the estimated p-values no longer evolve after the Bconv latest permutations. graph is a boolean indicating whether the estimated p-values have to be plotted against the number of permutations.
x	An object of class "testHSIC" storing the parameters and results of independence testing.
ylim	A vector of two values specifying the y-coordinate plotting limits.
err	A scalar value (between 0 and 1) specifying the reference type I error. This value is used to plot a vertical line.
...	Additional options.

Details

Why and how to properly choose kernels?

For a given input-output pair of variables, the Hilbert-Schmidt independence criterion (HSIC) is a dissimilarity measure between the joint bivariate distribution and the product of marginal distributions. Dissimilarity between those two distributions is measured through the squared norm of the distance between their respective embeddings in a reproducing kernel Hilbert space (RKHS) that directly depends on the selected input kernel K_{X_i} and the selected output kernel K_Y.

It must always be kept in mind that this criterion allows to detect independence within the pair (X_i, Y) provided that the two kernels are characteristic.

• If both kernels are characteristic, H_0: "X_i and Y are independent" is equivalent to H_0: "HSIC(X_i, Y)=0" and any estimator of HSIC(X_i, Y) emerges as a relevant test statistic.

• If they are not, testing H_0: "HSIC(X_i, Y)=0" is no longer sufficient for testing H_0: "X_i and Y are independent".

The reader is referred to Fukumizu et al. (2004) for the mathematical definition of a characteristic kernel and to Sriperumbur et al. (2010) for an overview of the major related results.

Moreover, in order to achieve the HSIC-ANOVA decomposition, specific kernels called ANOVA kernels must be used for the input variables. See Da Veiga (2021) for further details on ANOVA kernels, especially examples of ANOVA kernels and several methods to transform a basic kernel into an ANOVA kernel. Using ANOVA kernels is a necessary condition in order to define rigorously first-order and total-order HSIC-ANOVA indices.

In a nutshell, two situations must be distinguished:

1. If the intended test of independence is built upon HSIC(X_i, Y), all kernels must be characteristic.

2. If the intended test of independence is built upon T_i, all kernels must be characteristic and ANOVA kernels must be assigned to all input variables.

Responsability for kernel selection is left to the user while calling the function sensiHSIC. Let us simply recall which kernels verify the above properties:

• Characteristic kernel on any discrete set: "categ" (categorical kernel).

• Characteristic kernels on \mathbb{R}: "laplace" (exponential kernel), "matern3" (Matern 3/2 kernel), "matern5" (Matern 5/2 kernel) and "rbf" (Gaussian kernel). They remain characteristic when they are restricted to a compact domain D within \mathbb{R}.

• Characteristic kernels on [0,1]: "sobolev1" (Sobolev kernel with smoothness parameter r=1) and "sobolev2" (Sobolev kernel with smoothness parameter r=2). In addition, the four abovementioned kernels (which are characteristic on \mathbb{R}) remain characteristic on [0,1]. Their ANOVA variants are also characteristic on [0,1].

• Characteristic ANOVA kernel for any discrete distribution: "categ_anova" (transformed categorical kernel).

• Characteristic ANOVA kernels for the uniform distribution on [0,1]: "sobolev1" (Sobolev kernel with smoothness parameter r=1), "sobolev2" (Sobolev kernel with smoothness parameter r=2), "laplace_anova" (transformed exponential kernel), "matern3_anova" (transformed Matern 3/2 kernel), "matern5_anova" (transformed Matern 5/2 kernel) and "rbf_anova" (transformed Gaussian kernel).

Which test method is most appropriate?

First of all, the user must decide whether H_0 has to be tested with HSIC(X_i,Y) or with T_i. It mainly depends on what is stored in the object sensi of class "sensiHSIC". When both test statistics are available, it is strongly recommended to choose T_i because the associated test procedures proved to be more powerful in most cases. Unfortunately, in many situations (especially when the input variables X_i are not mutually independent), T_i cannot be estimated and HSIC(X_i,Y) is the only available option.

Tests of independence based on HSIC indices

The test statistic for the pair (X_i, Y) is either the U-statistic or the V-statistic associated to HSIC(X_i, Y).

If a V-statistic was used in sensiHSIC, four different test methods can be considered.

• The asymptotic test can be used if the sample size n is large enough (at least a hundred of samples). The asymptotic distribution of the test statistic is approximated by a Gamma distribution whose parameters are estimated with the method of moments. See Gretton et al. (2007) for more details about how to estimate the first two moments of the asymptotic Gamma distribution.

• The permutation-based test is more expensive in terms of computational cost but it can be used whatever the sample size n is. The initial output samples (stored in the object of class sensiHSIC) are randomly permuted B times and the test statistic is recomputed as many times. This allows to simulate B observations of the test statistic under H_0 and to estimate the p-value in a non-parametric way. See Meynaoui (2019) for more details on how to correctly estimate the p-value in order to preserve the expected level of the test.

• The sequential permutation-based test is a goal-oriented variant of the previous test. The main idea is to reduce the computational cost by stopping permutations as soon as the estimation of the p-value has sufficiently converged so that it can be compared to a reference threshold or be given a final ranking. See El Amri and Marrel (2022) for more details on how to implement this sequential approach for the three stopping criteria (namely "ranking", "screening" or "both").

• The Gamma test is a parametric alternative to permutation-based tests when n is not large enough to resort to the asymptotic test. The permutation-based test reveals the test statistic under H_0 follows a unimodal distribution having significant positive skewness. Thus, it seems quite natural to estimate the p-value with a Gamma distribution, especially in view of the fact that the asymptotic distribution is properly approximated by this parametric family. See El Amri and Marrel (2021) for more details on how to estimate the parameters of the Gamma distribution with the method of moments. In particular, the first two moments under H_0 of the test statistic are computed thanks to the formulas that were initially provided in Kazi-Aoual et al. (1995).

If a U-statistic was used in sensiHSIC, the estimated value of HSIC(X_i,Y) may be negative.

• The asymptotic test can no longer be conducted with a Gamma distribution (whose support is limited to [0,+\infty[). It is replaced by a Pearson III distribution (which is a left-shifted Gamma distribution).

• The permutation-based tests can be applied directly.

• The Gamma test has no longer any theoretical justification.

Tests of independence based on total-order HSIC-ANOVA indices

The test statistic for the pair (X_i, Y) is either the U-statistic or the V-statistic associated to T_i.

If a V-statistic was used in sensiHSIC, three different test methods can be considered.

• The permutation-based test is a simple adaptation of what is done for HSIC(X_i,Y). A slightly different permutation scheme is repeated B times and the test statistic is recomputed as many times. This allows to simulate B observations of the test statistic under H_0 and to estimate the p-value in a non-parametric way. This test may be very expensive if n is large.

• The sequential permutation-based test is the goal-oriented variant of the previous test. It is a simple adaptation of what is done for HSIC(X_i,Y). It allows to reduce as much as possible the computational load induced while recomputing the test statistic for a series of random permutations.

• The Gamma test is a parametric procedure that may be used when the computational load of permutation-based tests becomes too heavy (which happens when the sample size n is large). It is an extension of what is done for HSIC(X_i,Y). In order to compute the first two moments of the test statistic under H_0, the analytical formulas proposed in Kazi-Aoual et al. (1995) were altered so as to comply with less accomodating assumptions.

If a U-statistic was used in sensiHSIC, the estimated value of T_i may be negative. Only permutation-based tests can be considered.

What about target and conditional HSIC indices?

In Marrel and Chabridon (2021), HSIC indices were adapted to target sensitivity analysis (thus becoming T-HSIC indices) and to conditional sensitivity analysis (thus becoming C-HSIC indices). Tests of independence can still be useful after estimating T-HSIC indices or C-HSIC indices.

• For T-HSIC indices, the null hypothesis is H_0: "X_i and w(Y) are independent" where w is the weight function selected in target and passed to the function sensiHSIC. Everything works just as for basic HSIC indices (apart from the fact that w is applied on the original output variable Y). Available test methods include "Asymptotic", "Permutation", "Seq_Permutation" and "Gamma" (for V-statistics only).

• For C-HSIC indices, the null hypothesis is H_0: "X_i and Y are independent if the event described in cond occurs". In this specific context, testing conditional independence is only relevant if the weight function is an indicator function. For this reason, if conditional independence has to be tested, the user must select type="indicTh" in cond while calling the function sensiHSIC. Let us recall that only V-statistic estimators can be used for C-HSIC indices. As a result, available test methods include "Asymptotic", "Permutation", "Seq_Permutation" and "Gamma".

Target sensitivity analysis can also be considered with total-order HSIC-ANOVA indices (when the context allows to define them, which requires among other things mutually independent input variables). If their numerators were estimated during the prior call to sensiHSIC, the null hypothesis H_0: "X_i and w(Y) are independent" can be tested with T_i. Associated test methods include "Tot_Permutation", "Tot_Seq_Permutation" and "Tot_Gamma" (for V-statistics only).

Conditional sensitivity analysis cannot be considered with total-order HSIC-ANOVA indices because the conditioning operation is feared to introduce statistical dependence among input variables, which forbids using HSIC-ANOVA indices.

Value

testHSIC returns a list of class "testHSIC". It contains test.method, B (for the permutation-based test), seq.options (for the sequential permutation-based test) and the following objects:

call	The matched call.
pval	The estimated p-values after testing independence for all input-output pairs.
prop	A vector of two strings. The first string indicates if the chosen test method is asymptotic or non-asymptotic. The second string indicates if the chosen test method is parametric or non-parametric.
family	Only if test.method is "Asymptotic", "Gamma" or "Tot_Gamma". A string indicating the parametric family used to estimate p-values.
param	Only if test.method is "Asymptotic", "Gamma" or "Tot_Gamma". A 2-column (resp. 3-column) matrix containing the parameters of the Gamma (resp. Pearson III) distributions used to estimate p-values.
Hperm	Only if test.method is "Permutation" or "Tot_Permutation". A B-column matrix containing simulated values of the test statistics after randomly permuting the output samples. Each column in Hperm corresponds to one random permutation.
paths	Only if test.method is "Seq_Permutation" or "Tot_Seq_Permutation". A matrix containing all estimated p-values over the sequential test procedure. The i-th row provides all estimates of the i-th p-value as the number of permutations increases. If one row ends with a sequence of missing values NA, it means permutations were stopped earlier for this input variable. This situation can only happen if criterion="screening" in the list seq.options.

Author(s)

Sebastien Da Veiga, Amandine Marrel, Anouar Meynaoui, Reda El Amri and Gabriel Sarazin.

References

Da Veiga, S. (2021), Kernel-based ANOVA decomposition and Shapley effects: application to global sensitivity analysis, arXiv preprint arXiv:2101.05487.

El Amri, M. R. and Marrel, A. (2022), Optimized HSIC-based tests for sensitivity analysis: application to thermalhydraulic simulation of accidental scenario on nuclear reactor, Quality and Reliability Engineering International, 38(3), 1386-1403.

El Amri, M. R. and Marrel, A. (2021), More powerful HSIC-based independence tests, extension to space-filling designs and functional data.

Fukumizu, K., Bach, F. R. and Jordan, M. I. (2004), Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces, Journal of Machine Learning Research, 5(Jan), 73-99.

Gretton, A., Fukumizu, K., Teo, C., Song, L., Scholkopf, B. and Smola, A. (2007), A kernel statistical test of independence, Advances in Neural Information Processing Systems, 20.

Kazi-Aoual, F., Hitier, S., Sabatier, R. and Lebreton, J. D. (1995), Refined approximations to permutation tests for multivariate inference, Computational Statistics & Data Analysis, 20(6), 643-656.

Marrel, A. and Chabridon, V. (2021), Statistical developments for target and conditional sensitivity analysis: application on safety studies for nuclear reactor, Reliability Engineering & System Safety, 214, 107711.

Meynaoui, A. (2019), New developments around dependence measures for sensitivity analysis: application to severe accident studies for generation IV reactors (Doctoral dissertation, INSA de Toulouse).

Sriperumbudur, B., Fukumizu, K. and Lanckriet, G. (2010), On the relation between universality, characteristic kernels and RKHS embedding of measures, Proceedings of the 13th International Conference on Artificial Intelligence and Statistics (pp. 773-780). JMLR Workshop and Conference Proceedings.

See Also

sensiHSIC, weightTSA

Examples

  
# Test case: the Matyas function.

n <- 300  # nb of samples
p <- 3    # nb of input variables (including 1 dummy variable)

########################################
### PRELIMINARY SENSITIVITY ANALYSIS ###
########################################

X <- matrix(runif(n*p), n, p)
sensi <- sensiHSIC(model=matyas.fun, X, 
                   kernelX="sobolev1", anova=list(obj="both", is.uniform=TRUE))
print(sensi)
plot(sensi)
title("GSA for the Matyas function")

#############################
### TESTS OF INDEPENDENCE ###
#############################

### HSIC indices ###

test.asymp <- testHSIC(sensi)

test.perm <- testHSIC(sensi, test.method="Permutation")

test.seq.screening <- testHSIC(sensi, test.method="Seq_Permutation")

test.seq.ranking <- testHSIC(sensi, test.method="Seq_Permutation", 
                             seq.options=list(criterion="ranking"))

test.seq.both <- testHSIC(sensi, test.method="Seq_Permutation", 
                          seq.options=list(criterion="both"))

test.gamma <- testHSIC(sensi, test.method="Gamma")

# comparison of p-values

res <- rbind( t(as.matrix(test.asymp$pval)), t(as.matrix(test.perm$pval)), 
              t(as.matrix(test.seq.screening$pval)), t(as.matrix(test.seq.ranking$pval)),
              t(as.matrix(test.seq.both$pval)), t(as.matrix(test.gamma$pval)) )

rownames(res) <- c("asymp", "perm", "seq_perm_screening", 
                   "seq_perm_ranking", "seq_perm_both", "gamma")
res

### total-order HSIC-ANOVA indices ###

test.tot.perm <- testHSIC(sensi, test.method="Tot_Permutation")

test.tot.seq.screening <- testHSIC(sensi, test.method="Tot_Seq_Permutation")

test.tot.seq.ranking <- testHSIC(sensi, test.method="Tot_Seq_Permutation", 
                             seq.options=list(criterion="ranking"))

test.tot.seq.both <- testHSIC(sensi, test.method="Tot_Seq_Permutation", 
                          seq.options=list(criterion="both"))

test.tot.gamma <- testHSIC(sensi, test.method="Tot_Gamma")

res <- rbind( t(as.matrix(test.tot.perm$pval)), 
              t(as.matrix(test.tot.seq.screening$pval)), 
              t(as.matrix(test.tot.seq.ranking$pval)),
              t(as.matrix(test.tot.seq.both$pval)), 
              t(as.matrix(test.tot.gamma$pval)) )

rownames(res) <- c("tot_perm", "tot_seq_perm_screening", 
                   "tot_seq_perm_ranking", "tot_seq_perm_both", "tot_gamma")
res

#####################
### VISUALIZATION ###
#####################

# simulated values of HSIC indices under H0 (random permutations)
Hperm <- t(unname(test.perm$Hperm))

# simulated values of total-order HSIC-ANOVA indices under H0 (random permutations)
tot.Hperm <- t(unname(test.tot.perm$Hperm))

for(i in 1:p){
  
  ttl <- paste0("First-order and total-order HSIC-ANOVA indices for X", i)
  
  ######################################
  ### FIRST-ORDER HSIC-ANOVA INDICES ###
  ######################################
  
  # histogram of the test statistic under H0 (random permutations)
  
  hist(Hperm[,i], probability=TRUE,
       nclass=70, main=ttl, xlab="", ylab="", col="cyan")
  
  xx <- seq(0, max(tot.Hperm[,i]), length.out=200)
  
  # asymptotic Gamma distribution
  
  shape.asymp <- test.asymp$param[i, "shape"]
  scale.asymp <- test.asymp$param[i, "scale"]
  
  dens.asymp <- dgamma(xx, shape=shape.asymp, scale=scale.asymp)
  
  lines(xx, dens.asymp, lwd=2, col="darkorchid")
  
  # finite-sample Gamma distribution
  
  shape.perm <- test.gamma$param[i, "shape"]
  scale.perm <- test.gamma$param[i, "scale"]
  
  dens.perm <- dgamma(xx, shape=shape.perm, scale=scale.perm)
  
  lines(xx, dens.perm, lwd=2, col="blue")
  
  ######################################
  ### TOTAL-ORDER HSIC-ANOVA INDICES ###
  ######################################
  
  # histogram of the test statistic under H0 (random permutations)
  
  hist(tot.Hperm[,i], probability=TRUE, add=TRUE,
       nclass=70, xlab="", ylab="", col="deeppink")
  
  # finite-sample Gamma distribution
  
  shape.tot.perm <- test.tot.gamma$param[i, "shape"]
  scale.tot.perm <- test.tot.gamma$param[i, "scale"]
  
  dens.tot.perm <- dgamma(xx, shape=shape.tot.perm, scale=scale.tot.perm)
  
  lines(xx, dens.tot.perm, lwd=2, col="red")
  
  ### legend ###
  
  txt.1 <- paste0("Histogram of HSIC(X", i, ",Y)")
  txt.11 <- "Asymptotic Gamma distribution"
  txt.12 <- "Finite-sample Gamma distribution"
  
  txt.2 <- paste0("Histogram of T", i, " = HSIC(X,Y) - HSIC(X", 
                  paste(setdiff(1:p, i), collapse=""), ",Y)")
  txt.21 <- "Finite-sample Gamma distribution"
  
  all.cap <- c(txt.1, txt.11, txt.12, txt.2, txt.21)
  all.col <- c("cyan", "darkorchid", "blue", "deeppink", "red")
  all.lwd <- c(7, 2, 2, 7, 2)
  
  legend("topright", legend=all.cap, col=all.col, lwd=all.lwd, y.intersp=1.3)

}
 

sensitivity documentation built on Nov. 5, 2025, 6:55 p.m.