nanhung / pksensi

An R package for applying global sensitivity analysis in physiologically based kinetic modeling
https://nanhung.github.io/pksensi/
GNU Lesser General Public License v3.0
5 stars 3 forks source link

Error in quantile.default(newX[, i], ...) #8

Closed breisfeld closed 2 years ago

breisfeld commented 2 years ago

Hi Nan-Hung,

I am trying to run pksensi on a pbpk model that I've adapted from one created by Frederic.

The following is my job file:

source("utils/mcsim_utils.R")

library(pksensi)

set.seed(31415)

#----------------------------------------------------------------
# Sensitivity Analysis
#----------------------------------------------------------------

model <- "dmf.model"

q <- "qunif"

condition <- c("InhMag = 10", "InhDur = 7", "InhStart = 2", "InhPeriod = 48",
               "PCAir = 1", "log_PCFat = 1", "Frac = 0.8", "kGut = 0.3", 
               "Kp_sc_vs = 0.2", "fSA_exposed = 0.7", "Ke = 1", 
               "Michaelis = 0", "Vmax = 10", "Km = 1", "CLH = 10")

params <- c("fub", "Ke", "CLH")  
outputs <- c("CVen")

fub <- 0.5 
Ke <- 1
CLH <- 1 

q.arg <- list(list(min = fub * 0.8, max = fub * 1.2),
              list(min = Ke * 0.8, max = Ke * 1.2),
              list(min = CLH * 0.8, max = CLH * 1.2),
              list(min = 0.8, max = 1)) 
x <- rfast99(params, n = 400, q = q, q.arg = q.arg, replicate = 10)
times <- seq(0, 48, 0.5)

out_sens <- solve_mcsim(x, mName = model,
                        params = params, time = times,
                        vars = outputs, condition = condition)

Running the file in RStudio gives the following:

MCSim v6.2.0

Copyright (c) 1993-2020 Free Software Foundation, Inc.

MCSim comes with ABSOLUTELY NO WARRANTY;
This is free software, and you are welcome to redistribute it
under certain conditions; see the GNU General Public License.

* Using `modeling/dmf/dmf.model' model in file "dmf.model.c" created by D:/Projects/MCSim_under_R/MCSim/mod.exe v6.2.0

Reading experiment 1.
Doing analysis - 0 SetPoints runs... 1 experiment  each
0 runs specified for SetPoint(). Reading entire file.

Wrote results to "simmc.out"
Done.

Error in quantile.default(newX[, i], ...) : 
  missing values and NaN's not allowed if 'na.rm' is FALSE

Attached is the model file.

I didn't see anything obvious in the code for run_mcsim, but I am not an R expert. There is no traceback, so I don't know where to begin.

Thanks.

-Brad

dmf.model.txt

nanhung commented 2 years ago

Hi Brad,

I tested the code and found that the issue is because of the setting of the time point. The given time points are times <- seq(0, 48, 0.5). But, the concentration outputs are 0 at the beginning. In addition, there is no obvious uncertainty output. The issue can be solved by resetting the time points like times <- seq(10, 48, 0.5). I also noticed that only fub in the testing parameters have impacts on CVen, you can have more testing by using the code I provided as follow,

Rplot1

Rplot2

devtools::install_github("nanhung/RMCSim")
library(RMCSim)
library(pksensi)

install_mcsim()
model <- "dmf.model" 
makemcsim(model = model, dir = "MCSim") # The model file is placed in the subfolder named "MCSim" in the working directory

q <- "qunif"
condition <- c("InhMag = 10", "InhDur = 7", "InhStart = 2", "InhPeriod = 48",
               "PCAir = 1", "log_PCFat = 1", "Frac = 0.8", "kGut = 0.3", 
               "Kp_sc_vs = 0.2", "fSA_exposed = 0.7", "Ke = 1", 
               "Michaelis = 0", "Vmax = 10", "Km = 1", "CLH = 10")

params <- c("fub", "Ke", "CLH")  
outputs <- c("CVen")
dist <- rep("Uniform", 3)

fub <- 0.5 
Ke <- 1
CLH <- 1 

q.arg <- list(list(min = fub * 0.8, max = fub * 1.2),
              list(min = Ke * 0.8, max = Ke * 1.2),
              list(min = CLH * 0.8, max = CLH * 1.2)) 
pk.times <- seq(0, 48, 0.5)

# PK plot -----------------------------------------------------------------

set.seed(1111)
out <- solve_mcsim(mName = model, params = params, vars = outputs,
                   monte_carlo = 1000, dist = dist, q.arg = q.arg, 
                   time = pk.times, condition = condition)

# Check PK plot 
pksim(out, vars = "CVen", xlab = "Time", ylab = "Conc.", legend = FALSE)

# Sensitivity analysis ----------------------------------------------------

set.seed(1234)
x <- rfast99(params = params, n = 512, q = q, q.arg = q.arg, replicate = 5) 

# Reproduce error message
times1 <- seq(0, 48, 0.5) 
out1 <- solve_mcsim(x, mName = model, params = params, time = times1, vars = outputs, condition = condition)
## Error in quantile.default(newX[, i], ...) : 
## missing values and NaN's not allowed if 'na.rm' is FALSE

# Change time points. Started from 10 hr
times2 <- seq(10, 48, 0.5) # Set from 10 to 48 hr since the initial time points have no uncertainty output
out2 <- solve_mcsim(x, mName = model, params = params, time = times2, vars = outputs, condition = condition)
check(out2)
plot(out2) # Only fub have impact on CVen

# Additional check --------------------------------------------------------
# Set the lower fub (0.5 -> 0.1)
fub <- 0.1 
q.arg <- list(list(min = fub * 0.8, max = fub * 1.2),
              list(min = Ke * 0.8, max = Ke * 1.2),
              list(min = CLH * 0.8, max = CLH * 1.2)) 
q.arg # fub: U(0.08, 0.12)

set.seed(1111)
out3 <- solve_mcsim(mName = model, params = params, vars = outputs,
                   monte_carlo = 1000, dist = dist, q.arg = q.arg, 
                   time = times1, condition = condition)
pksim(out3)

# Set Ke (1 -> 10)
Ke <- 10 
q.arg <- list(list(min = fub * 0.8, max = fub * 1.2),
              list(min = Ke * 0.8, max = Ke * 1.2),
              list(min = CLH * 0.8, max = CLH * 1.2)) 
q.arg # Ke: U(8, 12)

set.seed(1111)
out4 <- solve_mcsim(mName = model, params = params, vars = outputs,
                   monte_carlo = 1000, dist = dist, q.arg = q.arg, 
                   time = times1, condition = condition)
pksim(out4)
# Does not have change on the output

# Set CLH (1 -> 100)
CLH <- 100 
q.arg <- list(list(min = fub * 0.8, max = fub * 1.2),
              list(min = Ke * 0.8, max = Ke * 1.2),
              list(min = CLH * 0.8, max = CLH * 1.2)) 
q.arg # CLH: U(80, 120)

set.seed(1111)
out5 <- solve_mcsim(mName = model, params = params, vars = outputs,
                   monte_carlo = 1000, dist = dist, q.arg = q.arg, 
                   time = times1, condition = condition)
pksim(out5)
# Does not have change on the output
breisfeld commented 2 years ago

Thanks, Nan-Hung!

This is a starter model that will be enhanced and should be more 'interesting' once some realistic data and parameters are included.

nanhung commented 2 years ago

Got it! Hoping this package can be helpful:)