sobolSalt function from sensitivity package

[jhupiterz]

Hi there,

I am trying to perform a variance-based sensitivity analysis on a modified Bayesian statistical model:

posterior = (prior x likelihood)/(prior x likelihood +(false_positive x (1-prior)))

Each term of the equation is a probability (i.e. comprised between 0 and 1) and I make the assumption they follow a uniform distribution. I would like to understand the influence of each parameter on the output of the model so I wish to compute the Sobol index (1st order and total index). Here's what I've written so far:

library(ggplot2)
library(sensitivity)
library(boot)

#-------------------------------------------------------------------------------------------------------#
# number of parameters combinations
n <- 100

# Matrix containing the n parameters combinations, dim = n * number of parameters in the model (here n*3)
X1 <- data.frame(matrix(runif(3 * n), nrow = n))
X2 <- data.frame(matrix(runif(3 * n), nrow = n))


# Bayesian inference when performing 1 test, elements in M[,1] = prior, M[,2] = likelihhood, and M[,3] = fasle_positive
SingleTest_Bayesian <- function(M) {
  posterior <- (M[,2]*M[,1])/((M[,2]*M[,1])+(M[,3]*(1-M[,1])))
  return(posterior)
}

# 1st order and total Sobol index calculation using sobolSalt
x <- sobolSalt(model = SingleTest_Bayesian, X1 = X1, X2 = X2, scheme = 'A', nboot = 0)
print(x)

So far this bit of code works but what I'm having trouble with is when I use an iterative Bayesian inference (i.e. for multiple tests in a row). For each parameter combination I want my function to output the final posterior probability and then compute the Sobol index on these final posterior probabilities.

# Iterative Bayesian inference with number_test = number of iterations 

MultipleTest_Bayesian <- function(M, number_test) {
  # initializing matrix containing all posterior values
  posterior <- data.frame(matrix(nrow=j, ncol=number_bio))   
  for (j in 1:n) {
    # initializing the first posterior for each row of M
    posterior[j,1] <- (M[j,2]*M[j,1])/((M[j,2]*M[j,1])+(M[j,3]*(1-M[j,1])))
    for (k in 2:number_bio) {
      # calculating the following posterior values 
      posterior[j,k] <- (M[j,2]*posterior[j,k-1])/((M[j,2]*posterior[j,k-1])+(M[j,3]*(1-posterior[j,k-1])))
    }
  }
  return(posterior[,number_test]) # returns the final postetior values 
}

# 1st order and total Sobol index calculation using sobolSalt
y <- sobolSalt(model = MultipleTest_Bayesian, X1 = X1, X2 = X2, number_test = 4, nboot = 0)
print(y)

And I always get the same error:

y <- sobolSalt(model = Catling_inference, X1 = X1, X2 = X2, number_bio = 4, nboot = 0)
Error in rep(0, p) : argument 'times' incorrect

In each case (single and multiple test) the functions return a numeric vector of length n so I don't understand why the SingleTest_Bayesian works and the MultipleTest_Bayesian doesn't... I checked the R documentation but the examples provided don't help me much. It must come from the fact I don't fully understand the sobolSalt function so if anyone could help me it would be much appreciated. I hope my post is clear enough. Thanks a lot