Sensitivity analysis for ODE with list as parameters - result gives standard deviation of 0

Question

Note: Initial problem was "Sensitivity analysis for ODE with parameters that include lists", as the sensRange-function gave an error due to the lists passed in the parameters. The question evolved as the list-parameters were fixed but a different problem where the sensitivity analysis showed strange results with a standard deviation of 0 for all runs.

I have a model simulating the concentrations of a chemical over time in R using the deSolve package. The parameters I use are the chemical properties (volume of distribution, halflife etc) as well as weight over time. Weight is given in a list which iterates over each time step in the ODE due to weight increase over time, and the chemical properties are given as a numeric.

I would like to perform a sensitivity analysis for this model, and only test the chemical properties. I have not done a sensitivity analysis before, but I am trying to follow examples using sensRange(). However, it doesn't seem like the sensRange() function allows that one of the parameters is given as a list. I get the error:

Error in yRef[, ivar] : invalid subscript type 'list'

My code for the model and global sensitivity analysis is set up like this:

library(FME)
library(deSolve)

c.weight <- c(3.5,  4,  5,  5,  6, 7,  7,  7,  8,  8,  8,  9,  9,  9,  9, 9, 10, 10, 10, 11, 11, 11, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 29, 29, 29, 30, 30, 30, 31, 31, 31, 32, 32, 32, 33, 33, 33, 34, 34, 34, 35, 35, 35, 36, 36, 36, 37, 37, 37, 38, 38, 38, 39, 39, 39, 40, 40, 40, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 51, 51, 52, 54, 55)

params.sens <- c(vd = 0.2,
                 pm = 0.05,
                 halflife = 5,
                 dose = 0.001,
                 c.weight = c.weight)

solve.sens <- function(pars) {
  sens.model <- function(times, state, parameters) {
    with(as.list(c(state, parameters)), {
      
      if (times <= 5) {
        volume <- (((0.2 * times) * c.weight[times+1]) * 30)
        transferred <- pm * volume
      } else if (times > 5 & times <= 14) {
        volume <- ((0.1 * c.weight[times+1]) * 30)
        transferred <- pm * volume
      } else {
        transferred <- 0
      }
      
      intake.c <- (dose * c.weight[times+1])
      elimination.c <- concentration * vd * c.weight[times+1] * log(2) / (halflife * 12)
      
      #total
      concentration <- (intake.c + transferred - elimination.c) / (vd * c.weight[times+1])
      
      list(c(concentration))
      
    })
  }
  
  state <- c(concentration = 0.5)
  months <- seq(0, 144, 1)
  return(as.data.frame(ode(y = state, times = months, func = sens.model, parms = params.sens)))
  
}

out <- solve.sens(params.sens)

parRanges <- data.frame(min = c(0.02, 0.1, 2.1), max = c(0.09, 0.2, 5.5))
rownames(parRanges) <- c("pm", "vd", "halflife")

sens <- sensRange(func = solve.sens, parms = params.sens, dist = "latin", sensvar = "concentration", parRange = parRanges, num = 50)

head(summary(sens))
summ.sens <- summary(sens)
plot(summ.sens, xlab = "months", ylab = "concentration")

I don't know how to go forward, does anyone have any tips or see where my mistake is??

Edit: followed the Bacterial growth model from Soetart and Herman, 2009, to correct my =/<- errors and added the parameter values from the comments into the model. Now it runs with no error however the summary shows identical values for all (mean, min, max, and all quantiles) so I am assuming it is not running correctly

               x      Mean Sd       Min       Max       q05       q25       q50       q75       q95
concentration0 0  0.500000  0  0.500000  0.500000  0.500000  0.500000  0.500000  0.500000  0.500000
concentration1 1  1.246348  0  1.246348  1.246348  1.246348  1.246348  1.246348  1.246348  1.246348
concentration2 2  3.475493  0  3.475493  3.475493  3.475493  3.475493  3.475493  3.475493  3.475493
concentration3 3  7.170403  0  7.170403  7.170403  7.170403  7.170403  7.170403  7.170403  7.170403
concentration4 4 12.314242  0 12.314242 12.314242 12.314242 12.314242 12.314242 12.314242 12.314242
concentration5 5 18.890365  0 18.890365 18.890365 18.890365 18.890365 18.890365 18.890365 18.890365

Answer 1

Thanks to fixing the assignment mistakes, the script runs now trough and the behavior can be reproduced.

Now we can see that the parameters passed to solve.sens were not passed down to the ode function.

A fix is to replace parms.sens in the ode call with pars, that was passed to the calling function:

return(as.data.frame(ode(y = state, times = months, func = sens.model, parms = pars)))

Then

plot(summ.sens, xlab = "months", ylab = "concentration")

results in:

Answer 2

The original post confused = and <- to create a named vector, so I recommended the following code snippet:

params.sens <- c(vd = p.vd,
              pm = p.pm,
              halflife = p.hl,
              dose = concentrations$dose,
              c.weight = variables$c.weight)

After the edit made by the original poster, this answer became mostly obsolete and was deleted. However, it turned out in a later post, that the distinction between <-(assignment) and = (parameter matching; here creation of a named vector) was not yet completely clear.

Here an example that shows the difference. To avoid confusion, run it in a fresh R session or delete the workspace before:

#rm(list=ls()) # uncomment this to clear the work space
x <- c(a <- 2, b <- 3)
y <- c(d = 2, e = 3)

x
y
ls()

where we see that y is a named vector, whereas x has no named elements. Instead a and b were "on the fly" created as variables in the user work space:

> x
[1] 2 3
> y
d e 
2 3 
> ls()
[1] "a" "b" "x" "y"