Introduction

This vignette presents a simple stochastic process model in which the quantity of interest is the time needed to reach a given number of successes. The model is defined at the level of individual units and is driven by exponential and Bernoulli random variables.

The goal is to compute Sobol indices for this model when the distributional parameters are uncertain.

Process model

The elementary unit is implemented by sobol4r_one_unit. The individual level model sobol4r_process_indiv aggregates units up to a target number of successes.

Sobol4R:::one_unit(
  lambda1 = 1 / 60,
  lambda2 = 1 / 1012,
  lambda3 = 1 / 72,
  p1      = 0.18,
  p2      = 0.5
)
#> [1] 0.000000000 0.005804484 0.000000000

process_fun_indiv(
  X_indiv = c(
    lambda1 = 1 / 60,
    lambda2 = 1 / 1012,
    lambda3 = 1 / 72,
    p1      = 0.18,
    p2      = 0.5
  ),
  M = 50
)
#> [1] 7.957008

Design for distributional parameters

We build two designs X1 and X2 for the uncertain distributional parameters, which are interpreted row wise by the process model.

n <- 200

draw_params <- function(n) {
  data.frame(t(replicate(
    n,
    c(
      1 / runif(1, min = 20,  max = 100),
      1 / runif(1, min = 24,  max = 2000),
      1 / runif(1, min = 24,  max = 120),
      runif(1,  min = 0.05, max = 0.3),
      runif(1,  min = 0.3,  max = 0.7)
    )
  )))
}

X1 <- draw_params(n)
X2 <- draw_params(n)

gensol_proc <- sobol4r_design(
  X1    = X1,
  X2    = X2,
  order = 2,
  nboot = 100
)
gensol_proc_s2007 <- sensitivity::sobol2007(
  model=NULL,
  X1    = X1,
  X2    = X2,
  nboot = 100
)

Sobol indices based on a single trajectory

Yproc_1 <- process_fun_row_wise(gensol_proc$X, M = 50)

Yproc_s2007 <- process_fun_row_wise(gensol_proc_s2007$X, M = 50)

xproc_1 <- tell(gensol_proc, Yproc_1)
xproc_s2007 <- tell(gensol_proc_s2007, Yproc_s2007)

print(xproc_1)
#> 
#> Call:
#> sensitivity::soboljansen(model = NULL, X1 = X1, X2 = X2, nboot = nboot)
#> 
#> Model runs: 1400 
#> 
#> First order indices:
#>       original       bias std. error  min. c.i. max. c.i.
#> X1 -0.06035885 -0.1943190  0.5252157 -0.5512794  1.530225
#> X2 -0.68748459 -0.2748701  0.8112947 -1.4527005  1.734165
#> X3 -0.62974779 -0.2672645  0.7794533 -1.3206662  1.772060
#> X4 -0.12271592 -0.2529045  0.7187522 -0.7039559  2.237898
#> X5 -0.56784187 -0.2390499  0.7183536 -1.2652686  1.515888
#> 
#> Total indices:
#>      original        bias std. error   min. c.i.  max. c.i.
#> X1 0.49688442 0.034074188 0.12012408 0.177370310 0.62348132
#> X2 0.05721835 0.003060691 0.01557337 0.018536193 0.08086854
#> X3 0.04479584 0.004307687 0.01549492 0.002083132 0.06269704
#> X4 0.79925259 0.074954219 0.23305675 0.253568296 1.08792579
#> X5 0.16104597 0.017976787 0.05480637 0.010849987 0.22749485
Sobol4R::autoplot(xproc_1, ncol = 1)

print(xproc_s2007)
#> 
#> Call:
#> sensitivity::sobol2007(model = NULL, X1 = X1, X2 = X2, nboot = 100)
#> 
#> Model runs: 1400 
#> 
#> First order indices:
#>       original          bias std. error  min. c.i.  max. c.i.
#> X1  0.43693012  0.0959874281 0.31896707 -0.4259164 0.79794331
#> X2 -0.07951464 -0.0026190458 0.04689638 -0.1701162 0.02444539
#> X3 -0.04666871  0.0009732841 0.03526860 -0.1127169 0.02859585
#> X4  0.85982140  0.1081248383 0.50227481 -0.4351312 1.34498153
#> X5  0.14986606  0.0198570238 0.11160263 -0.1782947 0.28665192
#> 
#> Total indices:
#>     original          bias std. error   min. c.i. max. c.i.
#> X1 0.6214603 -0.0228249164 0.15835255  0.40201418 1.0570390
#> X2 0.1811083 -0.0002272604 0.04834517  0.06796494 0.2760796
#> X3 0.1499875 -0.0124287828 0.04986474  0.07810300 0.2869787
#> X4 0.3287415 -0.0448496063 0.21856534  0.05531533 0.8203062
#> X5 0.1922343 -0.0135348126 0.11601118 -0.01730357 0.4411842
Sobol4R::autoplot(xproc_s2007)

Sobol indices based on a quantity of interest

Instead of relying on a single trajectory, we can define a quantity of interest for each parameter vector, for instance the mean time to reach the target number of successes. This is implemented by sobol4r_process_qoi.

Yproc_2 <- process_fun_mean_to_M(gensol_proc$X, M = 50)

xproc_2 <- tell(gensol_proc, Yproc_2)

print(xproc_2)
#> 
#> Call:
#> sensitivity::soboljansen(model = NULL, X1 = X1, X2 = X2, nboot = nboot)
#> 
#> Model runs: 1400 
#> 
#> First order indices:
#>       original         bias std. error  min. c.i. max. c.i.
#> X1 -0.02101140 -0.002997617  0.3959696 -0.6270899 0.9668445
#> X2 -0.63854756 -0.013630715  0.5362616 -1.3849414 0.7294311
#> X3 -0.64280818 -0.012788184  0.5341491 -1.3804438 0.7296107
#> X4 -0.03407319 -0.024270442  0.4122531 -0.5787039 1.1992003
#> X5 -0.48849339 -0.028587104  0.5161051 -1.2116467 0.8876513
#> 
#> Total indices:
#>       original          bias  std. error   min. c.i.   max. c.i.
#> X1 0.525063851  2.379760e-02 0.116238613 0.160381030 0.682258078
#> X2 0.005804016 -6.200323e-05 0.000864052 0.004161653 0.007676320
#> X3 0.005227879 -5.736972e-05 0.001330573 0.002548517 0.007647047
#> X4 0.748342102  1.539982e-02 0.194466556 0.192483589 1.042712257
#> X5 0.134914602  1.893345e-03 0.026865972 0.071892230 0.172516716
Sobol4R::autoplot(xproc_2, ncol = 1)

Qoi Mean

res_sobol_mean <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = base::mean,
  nrep    = 1000, 
  order   = 2,
  nboot   = 20,
  M       = 50,
  type  = "sobol"
)
print(res_sobol_mean)
#> 
#> Call:
#> sensitivity::sobol(model = NULL, X1 = X1, X2 = X2, order = order,     nboot = nboot)
#> 
#> Model runs: 3200 
#> 
#> Sobol indices
#>         original         bias std. error   min. c.i. max. c.i.
#> X1     0.6822643  0.034951868  0.2455683  0.06838191 1.0551364
#> X2     0.3008562 -0.030890270  0.1886502 -0.11682758 0.6020245
#> X3     0.3001573 -0.030759329  0.1883693 -0.11688103 0.6008908
#> X4     0.4364389 -0.011430001  0.1439652  0.19268134 0.6272249
#> X5     0.3134053 -0.033600451  0.1647855  0.01660074 0.5838436
#> X1*X2 -0.3027257  0.030599139  0.1894415 -0.60331095 0.1170975
#> X1*X3 -0.2996425  0.031135020  0.1876340 -0.59835804 0.1121831
#> X1*X4 -0.2036307  0.010896313  0.2599571 -0.55088006 0.4082261
#> X1*X5 -0.2662471  0.007527398  0.1998650 -0.58671563 0.1048185
#> X2*X3 -0.2995978  0.030972645  0.1888865 -0.60147915 0.1176453
#> X2*X4 -0.3010065  0.031437567  0.1885819 -0.60193172 0.1174187
#> X2*X5 -0.2990000  0.030879817  0.1882249 -0.59955878 0.1193742
#> X3*X4 -0.3022471  0.031121540  0.1882804 -0.60228288 0.1145267
#> X3*X5 -0.3002707  0.030517425  0.1878131 -0.59997343 0.1155203
#> X4*X5 -0.2780720  0.022580157  0.2088999 -0.59855564 0.1780694
Sobol4R::autoplot(res_sobol_mean, ncol = 1)

res_sobol2007_mean <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = base::mean,
  nrep    = 1000, 
  nboot   = 20,
  M       = 50,
  type  = "sobol2007"
)
print(res_sobol2007_mean)
#> 
#> Call:
#> sensitivity::sobol2007(model = NULL, X1 = X1, X2 = X2, nboot = nboot)
#> 
#> Model runs: 1400 
#> 
#> First order indices:
#>         original         bias  std. error    min. c.i.   max. c.i.
#> X1  5.934945e-01 3.907145e-02 0.332062917 -0.228291092 0.931882225
#> X2 -7.446782e-04 1.503936e-04 0.001239672 -0.003133059 0.002356756
#> X3 -3.304737e-05 4.999201e-05 0.001038250 -0.002016916 0.002641938
#> X4  8.513953e-01 1.439132e-01 0.562921530 -1.083195460 1.307876120
#> X5  1.650172e-01 8.834170e-03 0.108949504 -0.107114907 0.304261753
#> 
#> Total indices:
#>        original          bias  std. error    min. c.i.     max. c.i.
#> X1  0.566461949  0.0122810503 0.082013001  0.426079700  0.7402123156
#> X2 -0.002497144  0.0001186913 0.001934865 -0.008189405  0.0006489632
#> X3 -0.005057377  0.0001839905 0.001668025 -0.008628433 -0.0027416368
#> X4  0.521584070 -0.0106300678 0.119975374  0.308361178  0.7568794388
#> X5  0.126198631 -0.0431350400 0.107725867  0.006095012  0.3463368820
Sobol4R::autoplot(res_sobol2007_mean)

Qoi Median

res_sobol_median <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = stats::median,
  nrep    = 1000, 
  order   = 2,
  nboot   = 20,
  M       = 50,
  type    = "sobol"
)
print(res_sobol_median)
#> 
#> Call:
#> sensitivity::sobol(model = NULL, X1 = X1, X2 = X2, order = order,     nboot = nboot)
#> 
#> Model runs: 3200 
#> 
#> Sobol indices
#>         original        bias std. error   min. c.i. max. c.i.
#> X1     0.6775352  0.06127344  0.2732840  0.15326335 1.1708096
#> X2     0.2967152  0.08375007  0.2504626 -0.27627983 0.6944652
#> X3     0.2969483  0.08421302  0.2506960 -0.27783804 0.6937520
#> X4     0.4306871  0.03354405  0.1850199  0.03238601 0.6973292
#> X5     0.3123924  0.07460259  0.1931403 -0.05570025 0.6054862
#> X1*X2 -0.2993177 -0.08358247  0.2506151 -0.69754945 0.2756497
#> X1*X3 -0.2982007 -0.08407894  0.2498971 -0.69276534 0.2737260
#> X1*X4 -0.1964990 -0.07645700  0.2852624 -0.68748562 0.4710141
#> X1*X5 -0.2670671 -0.08131557  0.2612472 -0.71638963 0.2667171
#> X2*X3 -0.2965184 -0.08358587  0.2501122 -0.69217565 0.2783826
#> X2*X4 -0.2946446 -0.08444970  0.2503829 -0.69208927 0.2764071
#> X2*X5 -0.2981067 -0.08439250  0.2505819 -0.69562433 0.2758800
#> X3*X4 -0.2979195 -0.08505411  0.2525015 -0.69773657 0.2807934
#> X3*X5 -0.2993837 -0.08543498  0.2522567 -0.69626785 0.2836709
#> X4*X5 -0.2711645 -0.08861232  0.2663254 -0.70065204 0.2664554
Sobol4R::autoplot(res_sobol_median, ncol = 1)

res_sobol2007_median <- sobol4r_qoi_indices(
  model = process_fun_row_wise,
  X1      = X1,
  X2      = X2,
  qoi_fun = stats::median,
  nrep    = 1000, 
  nboot   = 20,
  M       = 50,
  type    = "sobol2007"
)
print(res_sobol2007_median)
#> 
#> Call:
#> sensitivity::sobol2007(model = NULL, X1 = X1, X2 = X2, nboot = nboot)
#> 
#> Model runs: 1400 
#> 
#> First order indices:
#>       original          bias  std. error    min. c.i.   max. c.i.
#> X1 0.590658606  0.0207370716 0.283222196 -0.204575989 0.923893887
#> X2 0.001081989  0.0001480158 0.001651758 -0.001485676 0.003946955
#> X3 0.001357549  0.0001640472 0.001586851 -0.003624942 0.002963604
#> X4 0.846054073 -0.0559033206 0.335420421  0.048460008 1.293452109
#> X5 0.165117072 -0.0027593259 0.104057674 -0.092891752 0.320207117
#> 
#> Total indices:
#>        original          bias  std. error    min. c.i.   max. c.i.
#> X1  0.570136419 -0.0136792398 0.073011144  0.483343764 0.731758331
#> X2  0.002311887 -0.0004408162 0.002922670 -0.004255009 0.007037544
#> X3 -0.001660063 -0.0014236614 0.003068034 -0.007959045 0.004917187
#> X4  0.525838025  0.0230391382 0.100506898  0.325803289 0.743804803
#> X5  0.130326087  0.0134582060 0.103782660 -0.019154379 0.363442151
Sobol4R::autoplot(res_sobol2007_median)

Summary

This vignette illustrates how Sobol4R can be used to compute Sobol indices for models in which both the mechanisms and their distributional parameters are stochastic. The same workflow can be applied to more complex simulators and to various quantities of interest (QoI).

Sobol4R: Sobol indices for a stochastic process model

Frédéric Bertrand

2025-11-27