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))
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.