[Python] How to tackle weird sobol indices plot in R

I am currently developping a model for transmission of health care acquire infection using desolve. In this model, HCWs are considered as vectors. I have implemented a Sobol's global sensitivity analysis of the model, however some of the parameters appear to not have values (or values are larger than one) (see plot 2 - sobol odd looking attached)

Could anybody provide help on how to perform the analysis in case i did anything wrong.

Below is a reproducible example.

library(deSolve)
library(tidyverse)
library(sensitivity)

model <- function(time, state, parameters){

with(as.list(c( state, parameters)),{

Sp <- state[1]
Cp <- state[2]
Ip <- state[3]
Un <- state[4]
Cn <- state[5]

# Parameters of patient flow and disease progression
Lambda <- parameters["Lambda"] #Admission rate
lambda_c <- parameters["lambda_c"] #Fraction of colonized
lambda_i <- parameters["lambda_i"] #Fraction of infected
x <- parameters["x"] #Fraction of Sp going to Ip

#Disease transmission
a <- parameters["a"] #Total number of contacts
b_np <- parameters["b_np"] #Transmission probability HCW-Patients
b_pn <- parameters["b_pn"] #Transmission probability patient-HCW
theta <- parameters["theta"] #Precaution compliance
tau_s <-parameters["tau_s"] #Discharge rate in Sp

kappa <- parameters["kappa"] #Transition rate from Cp
m_1 <-parameters["m_1"] #Probability of becoming infectious
#phi <- parameters["phi"] #Patient to HCW ratio - to be defined separately as Np/Nn
psi <- parameters["psi"] #Treatment rate
iota <- parameters["iota"] #Treatment success probability
mu_i <- parameters["mu_i"] #Death rate in Ip

#Decontamination
decon_min <- parameters["decon_min"]
decon_max <- parameters["decon_max"]

#Derived quantities
tau_c <- kappa * (1 - m_1) #Discharge rate in Cp including death from other causes
gamma_c <- kappa * m_1 #Rate of progression to Cp -> Ip

Nn <- Cn + Un
Np <- Sp + Cp + Ip

phi <- Np / Nn

beta_np <- a * b_np * (1 - theta) * Cn/Nn
beta_pn <- a * b_pn * (1 - theta) * phi * (Cp + Ip)/Np

Psi <- psi * iota #Treatment discharge rate

Lambda <- tau_s * Sp + tau_c * Cp + (Psi + mu_i) * Ip
#ODE patients

dSp <- (1 - lambda_c - lambda_i) * Lambda - beta_np * Sp - tau_s * Sp
dCp <- lambda_c * Lambda + (1 - x) * beta_np * Sp - (gamma_c + tau_c) * Cp
dIp <- lambda_i * Lambda + x * beta_np * Sp + gamma_c * Cp - (Psi + mu_i) * Ip

#ODE HCWs
dUn <- - beta_pn * Un + ((1 - theta) * decon_min + theta * decon_max) * Cn
dCn <- beta_pn * Un - ((1 - theta) * decon_min + theta * decon_max) * Cn

return(list(c(dSp, dCp, dIp, dUn, dCn)))

})
}

parameters <- c(

# Parameters of patient flow and disease progression
lambda_c = 0.4, #Fraction of colonized
lambda_i = 0.001, #Fraction of infected
x = 0.05, #Fraction of Sp going to Ip
a = 8, #Total number of contacts
b_np = 0.01, #Transmission probability HCW-Patients
b_pn = 0.15, #Transmission probability patient-HCW
theta = 0.4, #Precaution compliance
tau_s = 0.2, #Discharge rate in Sp

kappa = 0.1428, #Transition rate from Cp
m_1 = 0.3, #Probability of becoming infectious
psi = 0.0714, #Treatment rate
iota = 0.65, #Treatment success probability
mu_i = 0.01, #Death rate in Ip

decon_min = 2,
decon_max = 4
)

state <- c(

Sp = 590,
Cp = 5,
Ip = 5,
Un = 148,
Cn = 2
)

times <- seq(0, 100, by = 1)

model_results <- ode(
y = state,
times = times,
func = model,
parms = parameters
)



And below is the code for the sobols sensitivity analysis I performed

require(sensitivity)

simulate_hospital_model <- function(parameters) {

# Solve the system
result <- ode(y = state, times = times, func = model, parms = parameters)

# Extract the max number of infected patients
max(result[, "Ip"])
}

set.seed(190793)
N <- 1000 # Number of Sobol samples

# Define parameter ranges (assuming uniform distribution between min/max values)
param_ranges <- data.frame(
lambda_c = runif(N, 0.2, 0.6), # Fraction of colonized
lambda_i = runif(N, 0.0005, 0.005), # Fraction of infected
x = runif(N, 0.01, 0.1), # Fraction of Sp going to Ip
a = runif(N, 5, 15), # Total number of contacts
b_np = runif(N, 0.005, 0.02),# Transmission probability HCW-Patients
b_pn = runif(N, 0.1, 0.3), # Transmission probability patient-HCW
theta = runif(N, 0.2, 0.85), # Precaution compliance
tau_s = runif(N, 0.1, 0.5), # Discharge rate in Sp

kappa = runif(N, 0.1, 0.2), # Transition rate from Cp
m_1 = runif(N, 0.1, 0.5), # Probability of becoming infectious
psi = runif(N, 0.05, 0.2), # Treatment rate
iota = runif(N, 0.3, 0.8), # Treatment success probability
mu_i = runif(N, 0.005, 0.03), # Death rate in Ip

decon_min = runif(N, 1, 3), # Minimum decontamination rate
decon_max = runif(N, 3, 10) # Maximum decontamination rate
)
# Generate two independent Sobol sequences
X1 <- param_ranges
X2 <- data.frame(matrix(runif(N * ncol(X1), min = 0, max = 1), ncol = ncol(X1)))
colnames(X2) <- colnames(X1)

sobol_result <- soboljansen(
model = function(X) apply(X, 1, simulate_hospital_model),
X1 = X1,
X2 = X2,
#scheme = "A",
nboot = 100
)

# Display results
print(sobol_result)

#require(here)
#png(file =here("results", "visualizations", "sobol_totals.png"))
plot(sobol_result, ylim = c(0,1))


As I said, I would expect to have at least first and total effect for each variable (see attached plot 1 ---> sobol good looking, this is the result of a previous model), which in not the case (see attached plot 2 ---> sobol odd looking), could someone guide me through this, is there anything wrong anywhere in the code?

I would appreciate any guidance and alternative solutions.

Thanks.

Continue reading...