qosa {sensitivity}

Quantile-oriented sensitivity analysis

Description

qosa implements the estimation of first-order quantile-oriented sensitivity indices as defined in Fort et al. (2016) with a kernel-based estimator of conditonal probability density functions closely related to the one proposed by Maume-Deschamps and Niang (2018). qosa also supports a kernel-based estimation of Sobol first-order indices (i.e. Nadaraya-Watson).

Usage

qosa(model = NULL, X1, X2 = NULL, type = "quantile", alpha = 0.1, split.sample = 2/3, 
nsample = 1e4, nboot = 0, conf = 0.95, ...)
## S3 method for class 'qosa'
tell(x, y = NULL, ...)
## S3 method for class 'qosa'
print(x, ...)
## S3 method for class 'qosa'
plot(x, ylim = c(0, 1), ...)
## S3 method for class 'qosa'
ggplot(data, mapping = aes(), ylim = c(0, 1), ..., environment
                 = parent.frame())

Arguments

model

a function, or a model with a predict method, defining the model to analyze.

X1

a random sample of the inputs used for the estimation of conditional probability density functions. If X2 is NULL, X1 is split in two samples, with the first split.sample proportion of observations assigned to X1 and the rest to X2.

X2

a random sample of the inputs used to evaluate the conditional probability density functions. If NULL, it is constructed with the last (1-split.sample) proportion of observations from X1, see above.

type

a string specifying which first-order sensitivity indices must be estimated: quantile-oriented indices (type="quantile") or Sobol' indices (type="mean").

alpha

if type="quantile" the quantile level.

split.sample

if X2=NULL the proportion of observations from X1 assigned to the estimation of conditional probability density functions.

nsample

the number of samples from the conditional probability density functions used to estimate the conditional quantiles (if type="quantile") or the conditional means (if type="mean").

nboot

the number of bootstrap replicates.

conf

the confidence level for confidence intervals.

x

a list of class "sobolrank" storing the state of the sensitivity study (parameters, data, estimates).

data

a list of class "sobolrank" storing the state of the sensitivity study (parameters, data, estimates).

y

a vector of model responses.

ylim

y-coordinate plotting limits.

mapping

Default list of aesthetic mappings to use for plot. If not specified, must be supplied in each layer added to the plot.

environment

[Deprecated] Used prior to tidy evaluation.

...

any other arguments for model which are passed unchanged each time it is called.

Details

Quantile-oriented sensitivty indices were defined as a special case of sensitivity indices based on contrast functions in Fort et al. (2016). The estimator used by qosa follows closely the one proposed by Maume-Deschamps & Niang (2018). The only difference is that Maume-Deschamps and Niang (2018) use the following kernel-based estimate of the conditional cumulative distribution function:

\hat{F}(y\Vert X=x) = \frac{ \sum_{i=1}^n K_{h_x}(x - X_i) \bold{1}\{Y_i< y\}}{\sum_{i=1}^n K_{h_x}(x - X_i)}

whereas we use

\hat{F}(y\vert X=x) = \frac{ \sum_{i=1}^n K_{h_x}(x - X_i) \int_{-\infty}^y K_{h_y}(t - Y_i)dt} {\sum_{i=1}^n K_{h_x}(x - X_i)},

meaning that \bold{1}\{Y_i< y\} is replaced by \int_{-\infty}^y K_{h_y}(t - Y_i)dt = \Phi (\frac{y-Y_i}{h_y}) where \Phi is the cumulative distribution function of the standard normal distribution (since kernel K is Gaussian). The two definitions thus coincide when h_y \rightarrow 0. Our formula arises from a kernel density estimator of the joint pdf with a diagonal bandwidth. In a future version, it will be genralized to a general bandwidth matrix for improved performance.

Value

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

call

the matched call.

X

a data.frame containing the design of experiments.

X1

a data.frame containing the design of experiments used for the estimation of conditional probability density functions.

X

a data.frame containing the design of experiments used for the evaluation of conditional probability density functions.

y

a vector of model responses.

S

the estimations of the Sobol' sensitivity indices.

Author(s)

Sebastien Da Veiga

References

Fort, J. C., Klein, T., and Rachdi, N. (2016). New sensitivity analysis subordinated to a contrast. Communications in Statistics-Theory and Methods, 45(15), 4349-4364.

Maume-Deschamps, V., and Niang, I. (2018). Estimation of quantile oriented sensitivity indices. Statistics & Probability Letters, 134, 122-127.

Examples

 
library(ks)
library(ggplot2)
library(boot)

# Test case : difference of two exponential distributions (Fort et al. (2016))
# We use two samples with different sizes
n1 <- 5000
X1 <- data.frame(matrix(rexp(2 * n1,1), nrow = n1))
n2 <- 1000
X2 <- data.frame(matrix(rexp(2 * n2,1), nrow = n2))
Y1 <- X1[,1] - X1[,2]
Y2 <- X2[,1] - X2[,2]
x <- qosa(model = NULL, X1, X2, type = "quantile", alpha = 0.1)
tell(x,c(Y1,Y2))
print(x)
ggplot(x)

# Test case : difference of two exponential distributions (Fort et al. (2016))
# We use only one sample
n <- 1000 # put n=10000 for more consistency
X <- data.frame(matrix(rexp(2 * n,1), nrow = n))
Y <- X[,1] - X[,2]
x <- qosa(model = NULL, X1 = X, type = "quantile", alpha = 0.7)
tell(x,Y)
print(x)
ggplot(x)

# Test case : the Ishigami function
# We estimate first-order Sobol' indices (by specifying 'mean')
# Next lines are put in comment because too long fro CRAN tests
#n <- 5000 
#nboot <- 50 
#X <- data.frame(matrix(-pi+2*pi*runif(3 * n), nrow = n))
#x <- qosa(model = ishigami.fun, X1 = X, type = "mean", nboot = nboot)
#print(x)
#ggplot(x)