[R] Global sensitivity indices using sensitivity package: sobol, sobol2002

[Adiaba]

Dear R users,
I have been trying to estimate global sensitivity indices such as the sobol
1st and 2nd order indices. I managed to obtain the PRCC.

The example presented in the sensitivity package on sobol2002 seems to work
well for linear models: for example: calculate y for given x values. 

However, when trying to apply this technique to dynamic models (SIR type),
the error messages just keep spinning from one angle to another. Even the
decoupled approach still gives errors. These may be likely due to
misunderstanding of the indices or the application of the approach to
dynamic models. Here is an example: 

rm(list=ls())#clears objects
##Triangular distribution used for generating samples
rtri <- function(n=1,min=0,max=1,ml=0.5) {
  if((ml<min)||(ml>max)) {
    stop("ml outside of range [min max]")
  }
    u <- runif(n)
  mode <- (ml-min)/(max-min) # "mode" defined in range [0 1] (rescaling will
be done last)

  s1 <- which(u<=mode)
  s2 <- which(u>mode)

  u[s1] <- sqrt(mode*u[s1])
  u[s2] <- 1-sqrt((1-mode)*(1-u[s2]))

  min+(max-min)*u
  }



###Loading required packages
#install.packages("odesolve")
library(odesolve)
#install.packages("sensitivity",dependencies=T)
library(sensitivity)




###function that generates responses
simfunzcom<-function(xnews1,So=4250, I=250, R=0, No=4500){
SIRdob1<- function(t, x, parms){
with(as.list(c(parms,x)),{
dS <- (xnews1[, 2]-xnews1[, 4]*N)*(S+R+I*xnews1[, 6]*(1-xnews1[,
9]))-xnews1[, 3]*S + xnews1[,7]*R-(xnews1[, 1]*I*S)/N
dI <- (xnews1[, 1]*I*S)/N + I*(xnews1[, 2]-xnews1[, 3]*N)*xnews1[,
9]*xnews1[, 6]-(xnews1[, 8]+xnews1[, 3]+xnews1[, 5])*I
dR <- xnews1[, 5]*I -(xnews1[, 3] + xnews1[,7])*R
dN<-(xnews1[, 2]-xnews1[, 4]*N)*N-xnews1[, 8]*I-(xnews1[, 2]-xnews1[,
4]*N)*(1-xnews1[, 6])*I
der <- c(dS, dI,dR,dN)
list(der)  # the output must be returned
}) # end of 'with'
}  # end of function definition
yout<-matrix(0,100,100)
parms1<-xnews1
dt<- seq(1,50,1)
inits1 <- c(S=xnews1[, 10], I=xnews1[, 11], R=xnews1[, 12],N=xnews1[, 13])
for(j in 1:100){
simulation2 <- as.data.frame(lsoda(inits=inits1, times=dt, funct=SIRdob1,
parms=parms1[j, ]))
attach(simulation2)
 yout[j]<-as.numeric(simulation2[50 , 3])
}
yout
}


####Input data 1
n<-10
lambda<-rtri(n,min=4,max=7,ml=5.8)
a<-rtri(n,min=0.51,max=0.87,ml=0.64)
b<-rtri(n,min=0.0001,max=0.01,ml=0.001)
phi<-rtri(n,min=0.000006,max=0.05,ml=0.00004)
v<-rtri(n,min=0.3,max=0.8,ml=0.5)
rho <-rtri(n,min=0.1,max=0.9,ml=0.5)
delta <-rtri(n,min=0.01,max=0.65,ml=0.2)
alpha<-rtri(n,min=0.0001,max=0.1,ml=0.005)
e<-rtri(n,min=0.4,max=1.3,ml=0.9)
xnews11<-data.frame(lambda, a, b, phi, v, rho, delta, alpha,
e,So=4250,Io=250,Ro=0,No=4500)




###input data 2
n<-10
lambda <- rtri(n,min=4,max=7,ml=5.8)
a<- rtri(n,min=0.51,max=0.87,ml=0.64)
b <- rtri(n,min=0.0001,max=0.01,ml=0.001)
phi <- rtri(n,min=0.000006,max=0.05,ml=0.00004)
v<- rtri(n,min=0.3,max=0.8,ml=0.5)
rho <- rtri(n,min=0.1,max=0.9,ml=0.5)
delta <- rtri(n,min=0.01,max=0.65,ml=0.2)
alpha<- rtri(n,min=0.0001,max=0.1,ml=0.005)
e<- rtri(n,min=0.4,max=1.3,ml=0.9)
xnews2<-data.frame(lambda, a,b, phi, v, rho, delta, alpha,
e,So=4250,Io=250,Ro=0,No=4500)

###sibol2002
finsibcom <- sobol2002(model = simfunzcom, X1 = xnews1, X2 = xnews2, nboot =
100)



####decoupled approach
sa <- sobol(model = simfunzcom, X1 = xnews11, X2 = xnews2,nboot = 100)
xnews3<-data.frame(lambda=sa$X$lambda, a=sa$X$a,b=sa$X$b, phi=sa$X$phi,
v=sa$X$v, rho=sa$X$rho, delta=sa$X$delta, alpha=sa$X$alpha, e=sa$X$e)
resp<-simfunzcom(xnews2)
tell(sa,resp)



Please kindly direct this analysis or suggest articles that have
successfully applied this technique. The emails of the authors of
sensitivity package don't seem to work. All emails came back..

All assistance and directions is appreciated.




--
View this message in context: http://r.789695.n4.nabble.com/Global-sensitivity-indices-using-sensitivity-package-sobol-sobol2002-tp4310570p4310570.html
Sent from the R help mailing list archive at Nabble.com.