How to set up a Monte Carlo simulation in an ODE to estimate uncertainty?-CodePudding

I am aiming to do a Monte Carlo simulation of my ODE to see the variance of possible outcomes and graph the solutions to visualize the results for the entire time-sequence. However, I am struggling to set up the simulation and can't seem to find solutions for my specific problem. Previously I've used the package pksensi, but I'm using R version 3.4.4 and not able to update on this computer - so pksensi is not available. I have also seen examples using the FME-package using the modCRL-function. I've followed and run this example: https://rdrr.io/cran/FME/man/modCRL.html, however I struggle to understand how this is set up and the output. I've also tried seeing simpler models for potential solutions like this link Implement a Monte Carlo Simulation Method to Estimate an Integral in R, but I still can't find a solution.

I have a code showing what I've tried - I get errors and understand that this is not set up correctly, but hopefully my logic is understandable. I have never set up a Monte Carlo simulation before (besides pksensi), so any help as to how to set it up would be extremely appreciated.

library(deSolve)

months <- seq(0, 144, 1)
all.df <- data.frame(months, 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))

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)
    
    concentration <- (intake.c   transferred - elimination.c) / (vd * c.weight[times 1])
    
    list(c(concentration))
    
  })
}

#Following you can see my attempt to perform a Monte Carlo simulation with specific parameters using rnorm()
#The state can vary as well as the different parameters in params
state <- c(concentration = 0.5 * rnorm(1000, mean = 0.7, sd = 0.2)) #parameter to be tested
params <- list(vd = rnorm(1000, mean = 0.2, sd = 0.02), #parameter to be tested
               pm = rnorm(1000, mean = 0.05, sd = 0.03), #parameter to be tested
               halflife = rnorm(1000, mean = 4, sd = 2), #parameter to be tested
               dose = 0.001,
               c.weight = all.df$c.weight)
out <- as.data.frame(ode(y = state, times = months, func = model, parms = params))

CodePudding user response：

Following the approach from an answer given a few days ago, we may write a function that encapsulates random number generation and model simulation. Then we can call this in a loop or with replicate or lapply. The result is then a list of simulations that can be further analyzed. We can start with a plot. Here, the first argument of the plot function is an output of ode and the second a list of ode outputs. The first can be a deterministic reference simulation, here we use just the first case:

simulation <- function() {
  
  ## use only a single concentration, i.e. set rnorm(n=1, ....)
  state <- c(concentration = 0.5 * rnorm(1, mean = 0.7, sd = 0.2))
  params <- list(vd = rnorm(1, mean = 0.2, sd = 0.02),
                 pm = rnorm(1, mean = 0.05, sd = 0.03),
                 halflife = rnorm(1, mean = 4, sd = 2),
                 dose = 0.001,
                 c.weight = all.df$c.weight)
  ode(y = state, times = months, func = model, parms = params)
}

## here we use 1:10 for testing, but you can also use 1:1000 if you want
outlist <- lapply(1:10, function(i) simulation())

## plot results
plot(outlist[[1]], outlist)

Further ideas by request in form of a new SO question.