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

Confidence intervals #7

Closed paulrbuckley closed 3 years ago

paulrbuckley commented 3 years ago

Hi,

More questions than an issue - I am firstly wondering whether it is possible using the 'plot' function to filter parameters for which to visualise?

Secondly, in the event I'd like to make my own visualisations using the 'out' list, I am wondering, say for total order sesitivity, how the confidence interval is calculated? In the $tSI object, I cannot find multiple values for some parameters at timepoints where the 'plot' function shows a CI.

Best,

Paul

nanhung commented 3 years ago

Hi Paul,

For the first question. Unfortunately, the current function does not provide the feature to filter the specific parameter to display. The only way you can do this is by using heat_check(). That will filter out the non-influential parameter based on your cut-off. The example can calculate the confidence interval I provided below.

# 1
## Example code for PK model
library(pksensi)

pbtk1cpt <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    dAgutlument = - kgutabs * Agutlument
    dAcompartment = kgutabs * Agutlument - ke * Acompartment
    dAmetabolized = ke * Acompartment
    Ccompartment = Acompartment / vdist * BW;
    list(c(dAgutlument, dAcompartment, dAmetabolized),
         "Ccompartment" = Ccompartment)
  })
}

params <- c("vdist", "ke", "kgutabs", "BW")
q <- c("qunif", "qunif", "qunif", "qnorm")
q.arg <- list(list(min = 0.55, max = 2.2),
              list(min = 0.1, max = 0.4),
              list(min = 1.0, max = 4.0),
              list(mean = 70.0, sd = 5.0))

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

outputs <- c("Ccompartment", "Ametabolized")
initState <- c(Agutlument = 10, Acompartment = 0, Ametabolized = 0)
t <- seq(from = 0.01, to = 24.01, by = 1)
out <- solve_fun(x, time = t, func = pbtk1cpt, initState = initState, outnames = outputs)

plot(out)

## Use heat_check, changing SI.cutoff to filter parameter
heat_check(out)
heat_check(out, SI.cutoff = 0.01)

# 2
## Select the parameter (vdist) and variable (Ccompartment) to plot 
## Use total order as example 
## the main and interaction are stored in mSI and iSI
tSI <- out$tSI[, "vdist", "Ccompartment"]
tCI <- out$tCI[, "vdist", "Ccompartment"]
tCI_lwr <- tSI - tCI
tCI_upr <- tSI + tCI

plot(t, tSI, type="l")
lines(t, tCI_lwr, lty=2)
lines(t, tCI_upr, lty=2)
paulrbuckley commented 3 years ago

Thanks so much.

Best

P