sensiHSIC conducts a sensitivity analysis where the impact of an input variable is defined in terms of the distance between the input/output joint probability distribution and the product of their marginals when they are embedded in a Reproducing Kernel Hilbert Space (RKHS). This distance corresponds to the Hilbert-Schmidt Independence Criterion (HSIC) proposed by Gretton et al. (2005) and serves as a dependence measure between random variables, see Da Veiga (2015) for an illustration in the context of sensitivity analysis. The use of universal kernels ensures equivalence between HSIC nullity and independence of X and Y. In this case, a statistical test of independence with HSIC measure as statistic can be built. H0 is "X and Y are independent", against H1: X and Y are dependent. P-value can be computed either with asymptotic approximation (Gamma approximation), either with permutation method. See Meynaoui et al. (2019) for details. sensiHSIC can also be used to perform a Target Sensitivity Analysis (TSA). The basic idea of TSA is to measure the influence of the inputs on a critical domain of the model output, see Spagnol et al. (2019), Raguet and Marrel (2018) for details.

sensiHSIC is located in package sensitivity. Please install and load package sensitivity before use.

sensiHSIC(model = NULL, X, target = NULL, kernelX = "rbf", paramX = NA,
kernelY = "rbf", paramY = NA, nboot = 0, conf = 0.95,
estimator.type = "V-stat", test.method = "Asymptotic", B = 5000,
crit.option = list(stop.criterion = "screening", alpha = 0.05, Bstart = 100,
Bfinal = 5000, Bbatch = 100, lo = 100), ...)
## S3 method for class 'sensiHSIC'
tell(x, y = NULL, ...)
## S3 method for class 'sensiHSIC'
print(x, ...)
## S3 method for class 'sensiHSIC'
plot(x, ylim = c(0, 1), ...)
## S3 method for class 'sensiHSIC'
ggplot(x, ylim = c(0, 1), ...)

install.packages("sensitivity", repo="http://cran.r-project.org", dep=T)
library(sensitivity)
library(ggplot2)
library(boot)

# Test case : the non-monotonic Sobol g-function
# Only one kernel is provided with default hyperparameter value

n <- 100
X <- data.frame(matrix(runif(8 * n), nrow = n))

# HSIC-based GSA and asymptotic independence test
x <- sensiHSIC(model = sobol.fun, X, kernelX = "raquad", kernelY = "rbf",
nboot = 100, test.method = "Asymptotic")
print(x)

# HSIC-based GSA and permutation independence test
x <- sensiHSIC(model = sobol.fun, X, kernelX = "rbf", kernelY = "rbf",
estimator.type = "U-stat", test.method = "Permutation")
print(x)

# HSIC-based GSA and independence test with optimized number of permutations
x <- sensiHSIC(model = sobol.fun, X, kernelX = "raquad", kernelY = "rbf",
nboot = 100, test.method = "Seq_Permutation",
crit.option = list(stop.criterion = "ranking", alpha = 0.05, Bstart = 100, Bfinal = 3000,
Bbatch = 100, lo = 100))
print(x)

# HSIC-Target GSA in case of given model
x <- sensiHSIC(model = sobol.fun, X, target = list(c = 0.4, upper = TRUE, type = "indicTh",
param = 1), kernelX = "raquad", kernelY = "rbf",nboot = 100, test.method = "Permutation")

# HSIC-target GSA in case of given data
x <- sensiHSIC(model = NULL, X, target = list(c = 0.4, upper = TRUE, type = "indicTh", param = 1),
kernelX = "raquad", kernelY = "rbf",nboot = 100, test.method = "Permutation")
y <- sobol.fun(X)
tell(x,y)

# Test case : the Ishigami function

# A list of kernels is given with default hyperparameter value
n <- 100
X <- data.frame(matrix(-pi+2*pi*runif(3 * n), nrow = n))
x <- sensiHSIC(model = ishigami.fun, X, kernelX = c("rbf","matern3","dcov"),
kernelY = "rbf")
print(x)
ggplot(x)

# A combination of kernels is given and a dummy value is passed for
# the first hyperparameter
x <- sensiHSIC(model = ishigami.fun, X, kernelX = c("ssanova1","matern3","dcov"),
paramX = c(1,2,1), kernelY = "ssanova1")
print(x)
ggplot(x)

# Example in case of given data
n <- 100
X <- data.frame(matrix(-pi+2*pi*runif(3 * n), nrow = n))
Y <- ishigami.fun(X)

x <- sensiHSIC(model = NULL, X)
tell(x,Y)
print(x)
ggplot(x)

Return Values: sensiHSIC returns a list of class "sensiHSIC", containing all the input arguments detailed before, plus the following components:

Details: The HSIC sensitivity indices are obtained as a normalized version of the Hilbert-Schmidt independence criterion:

Si = HSIC(Xi,Y) / (√ HSIC(Xi,Xi) √ HSIC(Y,Y)),

see Da Veiga (2014) for details. When kernelX="dcov" and kernelY="dcov", the kernel is given by k(u,u')=1/2(||u||+||u'||-||u-u'||) and the sensitivity index is equal to the distance correlation introduced by Szekely et al. (2007) as was recently proven by Sejdinovic et al. (2013). The target sensitivity measures are defined via weight functions w which depend on c = threshold. Indicator function 1_c and smooth relaxations are available (according to target$type):

if type = "indicTh" --> w = 1_(Y>c),

if type = "logistic" --> w = 1/(1 + exp(-param*(Y-c)/abs(c))), cf. Spagnol et al. (2019)

if type = "exp1side" --> w = exp{-max(c - Y, 0)/(param.σ(Y)/5)}, cf. Raguet and Marrel (2018)

where σ(Y) is an estimation of the standard deviation of Y and param = 1 is a parameter tuning the smoothness.

See Also: kde, sensiFdiv