Sobol4R: Sobol indices for stochastic models

Introduction

This vignette illustrates the use of the Sobol4R package to compute Sobol indices for deterministic and stochastic versions of the classical Sobol g function.

The designs are generated with sensitivity::sobol and the models are provided by Sobol4R, in particular

  • sobol_g_function for the deterministic g function
  • sobol_g2_additive_noise for a version of the model with additive Gaussian noise
  • sobol_g2_with_covariate_noise for a version with covariate dependent noise

The stochastic models can be combined with the generic quantity of interest wrapper sobol4r_qoi.

Deterministic Sobol g function

Order 1 and Total via sensitivity::sobol()

n  <- 1e4
p  <- 8
X1 <- data.frame(matrix(runif(p * n), nrow = n))
X2 <- data.frame(matrix(runif(p * n), nrow = n))

sob_det <- sobol4r_design(X1 = X1, X2 = X2, order = 2, nboot = 50)

Y <- sobol_g_function(sob_det$X)
sensitivity::tell(sob_det, Y)
print(sob_det)
Sobol4R::autoplot(sob_det, ncol = 1)

Order 1 and Total via sensitivity::sobol2007()

n  <- 1e4
p  <- 8
X1 <- data.frame(matrix(runif(p * n), nrow = n))
X2 <- data.frame(matrix(runif(p * n), nrow = n))

sob_det_2007 <- sobol4r_design(
  X1    = X1,
  X2    = X2,
  nboot = 50,
  type  = "sobol2007"
)

Y <- sobol_g_function(sob_det_2007$X)
sensitivity::tell(sob_det_2007, Y)
print(sob_det_2007)
Sobol4R::autoplot(sob_det_2007)

Random effect on the output

We restrict the g function to the first two inputs and add a Gaussian noise term with zero mean and unit variance.

Order 1 and Total via sensitivity::sobol()

sob_noise_add <- sobol4r_design(X1 = X1[, 1:2], X2 = X2[, 1:2], order = 2, nboot = 50, type = "sobol")

Y <- sobol_g2_additive_noise(sob_noise_add$X)
sensitivity::tell(sob_noise_add, Y)
print(sob_noise_add)
Sobol4R::autoplot(sob_noise_add)

Order 1 and Total via sensitivity::sobol2007()

sob_noise_add <- sobol4r_design(X1 = X1[, 1:2], X2 = X2[, 1:2], nboot = 50, type = "sobol2007")

Y <- sobol_g2_additive_noise(sob_noise_add$X)
sensitivity::tell(sob_noise_add, Y)
print(sob_noise_add)
Sobol4R::autoplot(sob_noise_add)

Quantity of interest based on repeated runs

Instead of a single noisy run, we can focus on a quantity of interest, here the conditional mean of the output given the inputs. This is approximated by repeated calls to the stochastic model.

Order 1 and Total via sensitivity::sobol()

sob_noise_const_qoi <- sobol4r_qoi_indices(
  model = sobol_g2_additive_noise,
  X1 = X1[, 1:2], 
  X2 = X2[, 1:2], 
  order = 2, 
  nboot = 50, 
  type = "sobol")
print(sob_noise_const_qoi)
Sobol4R::autoplot(sob_noise_const_qoi)

Order 1 and Total via sensitivity::sobol2007()

sob_noise_const_qoi <- sobol4r_qoi_indices(
  model = sobol_g2_additive_noise,
  X1 = X1[, 1:2], 
  X2 = X2[, 1:2], 
  nboot = 50, 
  type = "sobol2007")
print(sob_noise_const_qoi)
Sobol4R::autoplot(sob_noise_const_qoi)

Covariate dependent noise

We now add a third input which controls the mean of the Gaussian noise term. The mean is equal to the third input, and the variance is fixed.

X1_cov <- data.frame(
  C1 = runif(n),
  C2 = runif(n),
  C3 = runif(n, min = 1, max = 100)
)
X2_cov <- data.frame(
  C1 = runif(n),
  C2 = runif(n),
  C3 = runif(n, min = 1, max = 100)
)

Order 1 and Total via sensitivity::sobol()

sob_cov_single <- sobol4r_design(X1 = X1_cov, X2 = X2_cov, order = 2, nboot = 50, type = "sobol")

Y <- sobol_g2_additive_noise(sob_cov_single$X)
sensitivity::tell(sob_cov_single, Y)
print(sob_cov_single)
Sobol4R::autoplot(sob_cov_single)

Order 1 and Total via sensitivity::sobol2007()

sob_cov_qoi <- sobol4r_qoi_indices(
  model = sobol_g2_with_covariate_noise,
  X1    = X1_cov,
  X2    = X2_cov,
  nboot = 50, 
  type = "sobol2007")
print(sob_cov_qoi)
Sobol4R::autoplot(sob_cov_qoi)

Conclusion

This vignette shows how Sobol4R can be used to study the impact of randomness in the model output on Sobol indices. More advanced examples, including models with random distributional parameters, are presented in a separate vignette.