# Interpretation of second-order sensitivity indices from 'sensitivity' package in R?

I am trying to understand the importance of input factors in constructing a composite indicator. I am using global sensitivity analyses in R's 'sensitivity' package. Using 'sobolSalt', there is a second-order calculation. I'm having trouble figure out how these second-order coefficients relate to the first- and total indices calculated by the same function. I would expect that:

Total Effect[k] = (First order)[k] + sum(second order interactions)[k] + sum(higher order interactions)[k]

where [k] is an individual input factor.

When I look at one specific example, and using the above equation, I find that:

0.8 = 0.72 + 0.16 + sum(Higher Order),

where the sum of the first and second order is already greater than the total effect. I've looked through the documentation but couldn't find an interpretation and also found the supplementary information in this article to be useful (in addition to Saltelli et al's 'Global Sensitivity Analysis: A Primer'):

Pons-Salort M, Serra-Cobo J, Jay F, López-Roig M, Lavenir R, et al. (2014) Insights into Persistence Mechanisms of a Zoonotic Virus in Bat Colonies Using a Multispecies Metapopulation Model. PLOS ONE 9(4): e95610. doi: 10.1371/journal.pone.0095610

Am I simply misunderstanding what the second-order coefficients are estimating with this function?

Here is a toy example of code I used to extract the pertinent bits of information:

# 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]
library(sensitivity)
library(boot)
set.seed(1234)
n <- 1000
X1 <- data.frame(matrix(runif(8 * n), nrow = n))
X2 <- data.frame(matrix(runif(8 * n), nrow = n))

x <- sobolSalt(model = sobol.fun, X1, X2, scheme="B", nboot = 100)
print(x)
plot(x, choice=1)

x$S # first order indices x$S2 # second order indices
x$T # total indices # Generating new variable names for second order indices for easy subsetting rownames(x$S)
n <- length(rownames(x$S)) second.varnames <- list() # Generating the names of the interactions in order of x$S2, but for easier subsetting
for (i in 1:n){
(var <- rownames(x$S)[i]) (others <- rownames(x$S)[(i+1):n])
(newvar <- paste(var,others))
second.varnames[[i]] <- newvar
}

# Concatenate, but don't include the last instance of second.varnames
second.order.names <- do.call("c", second.varnames[-length(second.varnames)])

first.order <- x$S second.order <- x$S2
total <- x\$T

rownames(second.order) <- second.order.names

(var <- rownames(first.order)[1]) # single variable
(var.Si <- first.order[var,1]) # first order index
(var.Ti <- total[var,1]) # total index
# get all of the second-order interactions
(var.int <- second.order[grep(paste(var, " ", sep=""),     rownames(second.order)),])

var.Si2 <- colSums(var.int)[1]

# I would expect: Total > (var.Si + var.Si2)


I think that you simply have a too important estimation error in this example. n = 1000 is not sufficient here to have a precise estimation error on second order indices

When increasing n to 100,000, I got

• First order index: 0.716
• Sum of second order indices: 0.0686
• Total index: 0.7873 > 0.7846 = 0.7162 + 0.0686

The sensitivity function seems to be actually computing second order indices from this example. It's just that it needs much more points to be estimated at the same precision as first order or total indices.