sobolrep: Sobol' indices estimation based on replicated orthogonal arrays

Description

sobolrep generalizes the estimation of the Sobol' sensitivity indices introduced by Tissot & Prieur (2015) using two replicated orthogonal arrays. This function estimates either

all first-order and second-order indices at a total cost of \(2 \times N\) model evaluations,
or all first-order, second-order and total-effect indices at a total cost of \(N \times (d+2)\) model evaluations,

where \(N=q^{2}\) and \(q \geq d-1\) is a prime number corresponding to the number of levels of each orthogonal array.

Usage

sobolrep(model = NULL, factors, N, tail=TRUE, 
			conf=0.95, nboot=0, nbrep=1, total=FALSE, ...)
# S3 method for sobolrep
tell(x, y = NULL, ...)
# S3 method for sobolrep
print(x, ...)
# S3 method for sobolrep
plot(x, ylim = c(0,1), choice, ...)

Value

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

call: the matched call.
X: a data.frame containing the design of experiments (row concatenation of the two replicated designs).
y: the response used.
RP: the matrix of permutations.
V: the model variance.
S: a data.frame containing estimations of the first-order Sobol' indices plus confidence intervals if specified.
S2: a data.frame containing estimations of the second-order Sobol' indices plus confidence intervals if specified.
T: a data.frame containing estimations of the total-effect indices plus confidence intervals if specified.

Arguments

model: a function, or a model with a predict method, defining the model to analyze.
factors: an integer giving the number of factors, or a vector of character strings giving their names.
N: an integer giving the size of each replicated design (for a total of \(2 \times N\) model evaluations).
tail: a boolean specifying the method used to choose the number of levels of the orthogonal array (see "Warning messages").
conf: the confidence level for confidence intervals.
nboot: the number of bootstrap replicates.
nbrep: the number of times the estimation procedure is repeated (see "Details").
total: a boolean specifying whether or not total effect indices are estimated.
x: a list of class "sobolrep" storing the state of the sensitivity study (parameters, data, estimates).
y: the model response.
ylim: y-coordinate plotting limits.
choice: an integer specifying which indices to plot: 1 for first-order indices, 2 for second-order indices, 3 for total-effect indices.
...: any other arguments for model which are passed unchanged each time it is called.

Warning messages

"The value entered for N is not the square of a prime number. It has been replaced by: ": the number of levels q of each orthogonal array must be a prime number. If N is not a square of a prime number, then this warning message indicates that it was replaced depending on the value of tail. If tail=TRUE (resp. tail=FALSE) the new value of N is equal to the square of the prime number preceding (resp. following) the square root of N.
"The value entered for N is not satisfying the constraint \(N \geq (d-1)^2\). It has been replaced by: ": the following constraint must be satisfied \(N \geq (d-1)^{2}\) where \(d\) is the number of factors. This warning message indicates that N was replaced by the square of the prime number following (or equals to) \(d-1\).

Details

sobolrep automatically assigns a uniform distribution on [0,1] to each input. Transformations of distributions (between U[0,1] and the wanted distribution) have to be performed before the call to tell() (see "Examples").

nbrep specifies the number of times the estimation procedure is repeated. Each repetition makes use of the orthogonal array structure to obtain a new set of Sobol' indices. It is important to note that no additional model evaluations are performed (the cost of the procedure remains the same).

References

A.S. Hedayat, N.J.A. Sloane and J. Stufken, 1999, Orthogonal Arrays: Theory and Applications, Springer Series in Statistics.

J-Y. Tissot and C. Prieur, 2015, A randomized orthogonal orray-based procedure for the estimation of first- and second-order Sobol' indices, J. Statist. Comput. Simulation, 85:1358-1381.

Examples

Run this code

# Test case: the non-monotonic Sobol g-function

# The method of sobol requires 2 samples
# (there are 8 factors, all following the uniform distribution on [0,1])

x <- sobolrep(model = sobol.fun, factors = 8, N = 1000, nboot=100, nbrep=1, total=FALSE)
print(x)
plot(x,choice=1)
plot(x,choice=2)

# Test case: dealing with non-uniform distributions

x <- sobolrep(model = NULL, factors = 3, N = 1000, nboot=0, nbrep=1, total=FALSE)

# X1 follows a log-normal distribution:
x$X[,1] <- qlnorm(x$X[,1])

# X2 follows a standard normal distribution:
x$X[,2] <- qnorm(x$X[,2])

# X3 follows a gamma distribution:
x$X[,3] <- qgamma(x$X[,3],shape=0.5)

# toy example
toy <- function(x){rowSums(x)}
y <- toy(x$X)
tell(x, y)
print(x)
plot(x,choice=1)
plot(x,choice=2)

Run the code above in your browser using DataLab