search

nlmixr2 / rxode2random

https://nlmixr2.github.io/rxode2random/
1 stars 0 forks source link

Simulation with IOV -- not accepting lotri matrix for omega #25

Closed SmithNM closed 3 months ago

SmithNM commented 3 months ago

Hello,

I'm trying to put together a few toy examples for training students, and I ran into an issue with implementing IOV where using a modified/simplified version of the Chapter 12 example is producing an error stating:

using C compiler: 'gcc.exe (GCC) 13.2.0'

Error: 'omega' needs to be a matrix or a numeric vector that can be converted to a matrix

The section of relevant code is as follows:

library(rxode2)
# library(ABxPKPD)
# library(ggplot2)
library(dplyr)
# library(tidyr)
# library(tibble)
library(lotri)

model.code = "

TABS = TV_TABS * exp(eta.TABS + iov.TABS)
TR_Fbio = TV_TR_Fbio + eta.TR_Fbio + iov.TR_Fbio

CL = TV_CL * exp(eta.CL)
V1 = TV_V1 * exp(eta.V1)
V2 = TV_V2 * exp(eta.V2)
CLD = TV_CLD * exp(eta.CLD)

KA = log(2) / (TABS/60)
FBIO = 1 / (exp(-TR_Fbio) + 1)

DC1 = AMT1/V1
DC2 = AMT2/V2

d/dt(AMTa) =       -KA * AMTa
d/dt(AMT1) = FBIO * KA * AMTa - CLD * DC1 + CLD * DC2 - CL * DC1
d/dt(AMT2) =                  + CLD * DC1 - CLD * DC2

d/dt(AUC) = DC1

"

Live.Model = rxode2(model = model.code)

n.sub.simulate = 10
n.trial.simulate = 1

theta = c("TV_TABS" = 45,
          "TV_TR_Fbio" = logit(x = 0.85),
          "TV_CL" = 10,
          "TV_V1" = 10,
          "TV_V2" = 65,
          "TV_CLD" = 25)

# tmp <- matrix(rnorm(6^2), 6, 6)
# tMat <- tcrossprod(tmp, tmp) / (6 ^ 2)
# dimnames(tMat) <- list(names(theta), names(theta))
# 
# tMat

omega = lotri(lotri(eta.TABS~0.25,
                    eta.TR_Fbio~0.20,
                    eta.CL~0.30,
                    eta.V1~0.30,
                    eta.V2~0.45,
                    eta.CLD~0.15) | id(nu=n.sub.simulate),
              lotri(iov.TABS~0.15,
                    iov.TR_Fbio~0.15) | occ(nu=n.sub.simulate*2))

omega

# omega.names = c(rownames(omega$id), rownames(omega$occ))
# 
# omega.means = rep(0,length=length(omega.names))
# names(omega.means) = omega.names
# 
# param = c(theta,omega.means)

dosing = et(amt=1000,
            addl=6,
            ii=24,
            evid=1,
            cmt="AMTa",
            time=0) %>%
  et(amt=1000,
     addl=6,
     ii=24,
     evid=4,
     cmt="AMTa",
     time = 336) %>%
  et(seq(0,168,0.5)) %>%
  et(seq(336,672,0.5)) %>%
  et(id=seq(1,n.sub.simulate))

dosing = mutate(dosing, occ = 1) %>%
  mutate(occ = ifelse(time>=336,2,occ))

Simulated.Data = rxSolve(object = Live.Model,
                         theta,
                         omega=omega,
                         ev=dosing,
                         nDisplayProgress=100L)

Any help on identifying what I might be missing would be greatly appreciated!

Thanks a ton!

mattfidler commented 3 months ago

I get what you get, though the example still runs fine.

Let me see what I can find out.

mattfidler commented 3 months ago

Thanks for reporting; There was a bug in the omega simulation code found in the rxode2random package.

If you update from R-universe, you can get the development version of rxode2:

install.packages(c("dparser", "rxode2ll", "rxode2parse",
                   "rxode2random", "rxode2et", "rxode2"),
                 repos=c(nlmixr2="https://nlmixr2.r-universe.dev",
                         CRAN="https://cloud.r-project.org"))

Then the code works for me:

library(rxode2)
#> rxode2 2.1.3 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
# library(ABxPKPD)
# library(ggplot2)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
# library(tidyr)
# library(tibble)
library(lotri)

model.code = "

TABS = TV_TABS * exp(eta.TABS + iov.TABS)
TR_Fbio = TV_TR_Fbio + eta.TR_Fbio + iov.TR_Fbio

CL = TV_CL * exp(eta.CL)
V1 = TV_V1 * exp(eta.V1)
V2 = TV_V2 * exp(eta.V2)
CLD = TV_CLD * exp(eta.CLD)

KA = log(2) / (TABS/60)
FBIO = 1 / (exp(-TR_Fbio) + 1)

DC1 = AMT1/V1
DC2 = AMT2/V2

d/dt(AMTa) =       -KA * AMTa
d/dt(AMT1) = FBIO * KA * AMTa - CLD * DC1 + CLD * DC2 - CL * DC1
d/dt(AMT2) =                  + CLD * DC1 - CLD * DC2

d/dt(AUC) = DC1

"

Live.Model = rxode2(model = model.code)
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

n.sub.simulate = 10
n.trial.simulate = 1

theta = c("TV_TABS" = 45,
          "TV_TR_Fbio" = logit(x = 0.85),
          "TV_CL" = 10,
          "TV_V1" = 10,
          "TV_V2" = 65,
          "TV_CLD" = 25)

# tmp <- matrix(rnorm(6^2), 6, 6)
# tMat <- tcrossprod(tmp, tmp) / (6 ^ 2)
# dimnames(tMat) <- list(names(theta), names(theta))
# 
# tMat

omega = lotri(lotri(eta.TABS~0.25,
                    eta.TR_Fbio~0.20,
                    eta.CL~0.30,
                    eta.V1~0.30,
                    eta.V2~0.45,
                    eta.CLD~0.15) | id(nu=n.sub.simulate),
                    lotri(iov.TABS~0.15,
                          iov.TR_Fbio~0.15) | occ(nu=n.sub.simulate*2))

omega
#> $id
#>             eta.TABS eta.TR_Fbio eta.CL eta.V1 eta.V2 eta.CLD
#> eta.TABS        0.25         0.0    0.0    0.0   0.00    0.00
#> eta.TR_Fbio     0.00         0.2    0.0    0.0   0.00    0.00
#> eta.CL          0.00         0.0    0.3    0.0   0.00    0.00
#> eta.V1          0.00         0.0    0.0    0.3   0.00    0.00
#> eta.V2          0.00         0.0    0.0    0.0   0.45    0.00
#> eta.CLD         0.00         0.0    0.0    0.0   0.00    0.15
#> 
#> $occ
#>             iov.TABS iov.TR_Fbio
#> iov.TABS        0.15        0.00
#> iov.TR_Fbio     0.00        0.15
#> 
#> Properties: nu

# omega.names = c(rownames(omega$id), rownames(omega$occ))
# 
# omega.means = rep(0,length=length(omega.names))
# names(omega.means) = omega.names
# 
# param = c(theta,omega.means)

dosing = et(amt=1000,
            addl=6,
            ii=24,
            evid=1,
            cmt="AMTa",
            time=0) %>%
et(amt=1000,
   addl=6,
   ii=24,
   evid=4,
   cmt="AMTa",
   time = 336) %>%
et(seq(0,168,0.5)) %>%
et(seq(336,672,0.5)) %>%
et(id=seq(1,n.sub.simulate))

dosing = mutate(dosing, occ = 1) %>%
mutate(occ = ifelse(time>=336,2,occ))

Simulated.Data = rxSolve(object = Live.Model,
                                theta,
                                omega=omega,
                                ev=dosing,
                                nDisplayProgress=100L)
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’

Simulated.Data
#> ── Solved rxode2 object ──
#> ── Parameters (x$params): ──
#> # A tibble: 10 × 17
#>    id    `ETA[1]` `ETA[3]` `ETA[2]` `ETA[4]` TV_TABS eta.TABS TV_TR_Fbio
#>    <fct>    <dbl>    <dbl>    <dbl>    <dbl>   <dbl>    <dbl>      <dbl>
#>  1 1        0.848   0.469  -0.200      0.183      45   0.570        1.73
#>  2 2        0.474   0.191   0.322      0.620      45  -1.34         1.73
#>  3 3       -0.242   0.449  -0.414     -0.550      45  -0.0430       1.73
#>  4 4       -0.177   0.829   0.00349   -0.317      45  -0.467        1.73
#>  5 5        0.463   0.0244  0.463      0.215      45   0.0791       1.73
#>  6 6       -0.153   0.371  -0.213      0.546      45   0.0283       1.73
#>  7 7        0.994  -0.0657  0.0911     0.324      45   1.11         1.73
#>  8 8       -0.179  -0.289  -0.0638     0.502      45   0.0222       1.73
#>  9 9        0.335   0.176  -0.00588   -0.813      45   0.443        1.73
#> 10 10      -0.206  -0.209   0.0725     0.283      45  -0.0612       1.73
#> # ℹ 9 more variables: eta.TR_Fbio <dbl>, TV_CL <dbl>, eta.CL <dbl>,
#> #   TV_V1 <dbl>, eta.V1 <dbl>, TV_V2 <dbl>, eta.V2 <dbl>, TV_CLD <dbl>,
#> #   eta.CLD <dbl>
#> ── Initial Conditions (x$inits): ──
#> AMTa AMT1 AMT2  AUC 
#>    0    0    0    0 
#> ── First part of data (object): ──
#> # A tibble: 10,100 × 19
#>      id  time iov.TABS iov.TR_Fbio  TABS TR_Fbio    CL    V1    V2   CLD    KA
#>   <int> <dbl>    <dbl>       <dbl> <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1     1   0      0.848      -0.200  186.    2.81  1.20 0.856  247.  11.8 0.224
#> 2     1   0.5    0.848      -0.200  186.    2.81  1.20 0.856  247.  11.8 0.224
#> 3     1   1      0.848      -0.200  186.    2.81  1.20 0.856  247.  11.8 0.224
#> 4     1   1.5    0.848      -0.200  186.    2.81  1.20 0.856  247.  11.8 0.224
#> 5     1   2      0.848      -0.200  186.    2.81  1.20 0.856  247.  11.8 0.224
#> 6     1   2.5    0.848      -0.200  186.    2.81  1.20 0.856  247.  11.8 0.224
#> # ℹ 10,094 more rows
#> # ℹ 8 more variables: FBIO <dbl>, DC1 <dbl>, DC2 <dbl>, AMTa <dbl>, AMT1 <dbl>,
#> #   AMT2 <dbl>, AUC <dbl>, occ <fct>

Created on 2024-06-06 with reprex v2.1.0

SmithNM commented 3 months ago

Brilliant! Thanks so much for the help and guidance, Matt!