Simulate a parameter that is not normally distributed (but rather a mixture from populations)

[mattfidler]

For rxode2 simulation, I want to sample a parameter value that is not normally distributed.

For instance, I tried the following.

kaPop <<- c(rnorm(n=0.2*1e6, mean = 1, sd = 0.1),
            rnorm(n=0.4*1e6, mean = 1.5, sd = 0.2),
            rnorm(n=0.4*1e6, mean = 2, sd = 0.3))

mod <- RxODE({
  # Each time, draw a uniformly distributed variable between 0 and 1, and select
  # the ka from a population
  p_val = runif(1,0,1)
  ka = quantile(kaPop, probs = p_val)
  
  A_dep(0)  =  1e4;

  d/dt(A)   = -ka*A;
  cp <- A/V * (1 + prop.sd) + add.sd 
})

It complains:

rxode2 model syntax error:
================================================================================
:001:     p_val = runif(1, 0, 1)

rxode2 syntax error:
:002:     ka = quantile(kaPop, probs = p_val)
                                     ^
:002:     ka = quantile(kaPop, probs = p_val)

rxode2 syntax error:
:003:     A_dep(0) = 10000
          ^
:003:     A_dep(0) = 10000

rxode2 syntax error:
:004:     d/dt(A) = -ka * A
            ^
:004:     d/dt(A) = -ka * A

rxode2 syntax error:
:005:     cp <- A/V * (1 + prop.sd) + add.sd
                     ^
================================================================================

[mattfidler]

You can simulate non-normal distributions outside the rxode2 solve and then supply them to the solving rxSolve

See

library(rxode2)
#> Warning: package 'rxode2' was built under R version 4.3.2
#> rxode2 2.1.1 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

set.seed(10)

kaPop <- c(rnorm(n=0.2*1e6, mean = 1, sd = 0.1),
            rnorm(n=0.4*1e6, mean = 1.5, sd = 0.2),
            rnorm(n=0.4*1e6, mean = 2, sd = 0.3))

pars <- do.call(rbind,
        lapply(seq_along(1:100), function(i) {
          u <- runif(1, 0, 1)
          data.frame(ka=quantile(kaPop, probs=u), V=10,
                     prop.sd=0.1, add.sd=0.1)  
        }))

mod <- function(){
  model({
    d/dt(A)   = -ka*A;
    A(0)  =  1e4;
    cp <- A/V * (1+ rnorm(0, prop.sd)) + rnorm(0, add.sd)
  })
}



s <- rxSolve(mod,params=pars, events=et(1:10, id=1:100))

print(s)
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#> # A tibble: 100 × 5
#>    id       ka     V prop.sd add.sd
#>    <fct> <dbl> <dbl>   <dbl>  <dbl>
#>  1 1      2.25    10     0.1    0.1
#>  2 2      1.70    10     0.1    0.1
#>  3 3      1.64    10     0.1    0.1
#>  4 4      1.89    10     0.1    0.1
#>  5 5      1.60    10     0.1    0.1
#>  6 6      1.52    10     0.1    0.1
#>  7 7      2.23    10     0.1    0.1
#>  8 8      1.24    10     0.1    0.1
#>  9 9      1.03    10     0.1    0.1
#> 10 10     1.81    10     0.1    0.1
#> # ℹ 90 more rows
#> ── Initial Conditions ($inits): ──
#>     A 
#> 10000 
#> ── First part of data (object): ──
#> # A tibble: 1,000 × 4
#>      id  time       cp         A
#>   <int> <dbl>    <dbl>     <dbl>
#> 1     1     1 104.     1058.    
#> 2     1     2  11.2     112.    
#> 3     1     3   1.12     11.8   
#> 4     1     4   0.185     1.25  
#> 5     1     5   0.0493    0.132 
#> 6     1     6   0.0296    0.0140
#> # ℹ 994 more rows

plot(s, cp, log="y")
#> Warning in self$trans$transform(x): NaNs produced
#> Warning: Transformation introduced infinite values in continuous y-axis
#> Warning: Removed 65 rows containing missing values (`geom_line()`).

^{Created on 2024-01-16 with reprex v2.0.2}

[mattfidler]

In my actual case, there are other parameters with random effects. Do I have to supply those parameters in the same pars data frame, or could I still specify the omega matrix and supply it to rxSolve function?

Hi ,

Yes you can, but not with the standard interface. You need to make the data.frame outside and then supply the parameters directly to rxode2’s rxSolve.

Here is the documentation for that:

https://nlmixr2.github.io/rxode2/articles/rxode2-pipeline.html?q=rxSimThetaOmega#simulating-from-a-data-frame-of-parameters