search

nlmixr2 / nlmixr2est

nlmixr2 estimation routines
https://nlmixr2.github.io/nlmixr2est/
GNU General Public License v2.0
7 stars 3 forks source link

Error - Drug/metabolite modeing after zero oder absorption. #394

Closed Lee21A closed 1 year ago

Lee21A commented 1 year ago

I tried to fit Parent and metabolite simultaneously after 0 order absorption of parent dug based on 1-compartment model.

km: conversion rate from parent to metabolite cl : clearance for parent clm : clearance for metabolite v : volume of distribution dur0: duration time of 0 order absorption

Step 1. I made a model like this.

library(nlmixr2)
one.0abs.PM <- function() {
  ini({
    tkm <- log(c(0,1)) 
    tcl <- log(c(0,2)) 
    tclm <- log(c(0,10))  
    tv <- log(c(0,10))  
    tdur0 <- log(c(0,1))  
    eta.cl ~ 1
    eta.v ~ 1
    eta.km ~ 1 
    eta.clm ~ 1
    eta.dur0 ~ 1
    prop.err <- 1 
    prop.err2 <- 1 
  })

  model({
    km <- exp(tkm + eta.km)
    cl <- exp(tcl + eta.cl)
    clm <- exp(tclm + eta.clm)
    v <- exp(tv + eta.v)
    dur0    <- exp(tdur0+eta.dur0)
    d/dt(center) = - cl / v * center - km * center
    dur(center)=dur0
    d/dt(meta) = km * center - clm/ v * meta
    center(0)= 0
    meta(0)= 0
    cp = center / v
    cm = meta / v
    cp ~ prop(prop.err)   
    cm ~  prop(prop.err2) 
  })
}
nlmixr2(one.0abs.PM)

The result was fine, and information for multiple endpoint was

── Multiple Endpoint Model ($multipleEndpoint): ──  
  variable               cmt               dvid*
1  cp ~ … cmt='cp' or cmt=3 dvid='cp' or dvid=1
2  cm ~ … cmt='cm' or cmt=4 dvid='cm' or dvid=2
  * If dvids are outside this range, all dvids are re-numered sequentially, ie 1,7, 10 becomes 1,2,3 etc

Step 2. I made a data file considering the upper information.

image

I used numbers for CMT and DVID based on the upper multiple endpoint model.

Step 3. I tried to fit. 

dat <- read.csv("0ABS_PM_Report3.csv",header = T,encoding = "UTF-8", na.strings=c("."))
fit <- nlmixr2(one.0abs.PM, dat, list(print=0), est="saem")

Step 4. The fitting result was weird, and there was a wrong assignment for DV data. (for example 3.8) image For every ID, DV (such as 3.8) of time 0.5 assigned as DVID=1 is “cm”, not “cp”. This weird phenomenon applied for only time 0.5.

Step 5. I changed the model like this and the CMT info, But it was the same result.

cp ~ prop(prop.err)   |center  --> changed the number of CMT 3 to 1
cm ~  prop(prop.err2) |meta -->  changed the number of CMT 4 to 2

Step 5. I tried to change the info of DVID of CMT to “cp”, “cm” or “center”, “meta”. But the result was the same or I encountered an error like this.

'DVID'/'CMT' translation:
'CMT': 4
Error : 'dvid'->'cmt' or 'cmt' on observation record on a undefined compartment (use 'cmt()' 'dvid()') id: 1 row: 0
Error: 'dvid'->'cmt' or 'cmt' on observation record on a undefined compartment (use 'cmt()' 'dvid()') id: 1 row: 0

Could somebody help me? I attached csv files I used. Thanks in advance.

Lee 0ABS_PM_Report3.csv

mattfidler commented 1 year ago

Hi @Lee21A

It seems the first step, the translation of the model to the input dataset, was correct. You can see it with fit$dataSav:

library(nlmixr2est)
#> Warning: package 'nlmixr2est' was built under R version 4.3.1
#> Loading required package: nlmixr2data
#> Warning: package 'nlmixr2data' was built under R version 4.3.1
one.0abs.PM <- function() {
  ini({
    tkm <- log(c(0,1)) 
    tcl <- log(c(0,2)) 
    tclm <- log(c(0,10))  
    tv <- log(c(0,10))  
    tdur0 <- log(c(0,1))  
    eta.cl ~ 1
    eta.v ~ 1
    eta.km ~ 1 
    eta.clm ~ 1
    eta.dur0 ~ 1
    prop.err <- 1 
    prop.err2 <- 1 
  })

  model({
    km <- exp(tkm + eta.km)
    cl <- exp(tcl + eta.cl)
    clm <- exp(tclm + eta.clm)
    v <- exp(tv + eta.v)
    dur0    <- exp(tdur0+eta.dur0)
    d/dt(center) = - cl / v * center - km * center
    dur(center)=dur0
    d/dt(meta) = km * center - clm/ v * meta
    center(0)= 0
    meta(0)= 0
    cp = center / v
    cm = meta / v
    cp ~ prop(prop.err)   
    cm ~  prop(prop.err2) 
  })
}

dat <- read.csv("c:/tmp/0ABS_PM_Report3.csv",header = T,encoding = "UTF-8", na.strings=c("."))
fit <- nlmixr2(one.0abs.PM, dat, list(print=0), est="saem")
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> Calculating covariance matrix
#> Error : Error calculating covariance via linearization
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> → optimizing duplicate expressions in saem predOnly model 0...
#> → finding duplicate expressions in saem predOnly model 1...
#> → optimizing duplicate expressions in saem predOnly model 1...
#> → finding duplicate expressions in saem predOnly model 2...
#> → optimizing duplicate expressions in saem predOnly model 2...
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 4208
#> → compress phiM in nlmixr2 object, save 63992
#> → compress parHist in nlmixr2 object, save 18048
#> → compress saem0 in nlmixr2 object, save 35848

# Fit at time =  0.5
print(fit %>% as.data.frame() %>% dplyr::filter(TIME==0.5))
#>   ID TIME CMT  DV     PRED         RES    IPRED        IRES      IWRES
#> 1  1  0.5  cm 3.8 1.052377  2.74762342 1.052458  2.74754182 15.7994069
#> 2  1  0.5  cm 0.9 1.052377 -0.15237658 1.052458 -0.15245818 -0.8766923
#> 3  2  0.5  cm 3.4 1.052377  2.34762342 1.052145  2.34785485 13.5050722
#> 4  2  0.5  cm 0.7 1.052377 -0.35237658 1.052145 -0.35214515 -2.0255705
#> 5  3  0.5  cm 2.8 1.052377  1.74762342 1.052563  1.74743673 10.0474203
#> 6  3  0.5  cm 1.1 1.052377  0.04762342 1.052563  0.04743673  0.2727519
#>         eta.cl         eta.v        eta.km       eta.clm      eta.dur0       cp
#> 1  0.002990595  2.090554e-07  3.036159e-04  4.454576e-07 -1.765297e-05 4.800266
#> 2  0.002990595  2.090554e-07  3.036159e-04  4.454576e-07 -1.765297e-05 4.800266
#> 3  0.002855822  1.039170e-06 -3.130864e-05 -1.773866e-06  4.417569e-07 4.800563
#> 4  0.002855822  1.039170e-06 -3.130864e-05 -1.773866e-06  4.417569e-07 4.800563
#> 5 -0.005186057 -2.373407e-08 -2.022258e-04 -2.352790e-08 -4.039316e-06 4.804493
#> 6 -0.005186057 -2.373407e-08 -2.022258e-04 -2.352790e-08 -4.039316e-06 4.804493
#>         cm   center     meta        km        cl       clm         v      dur0
#> 1 1.052458 4.766634 1.045084 0.9084866 0.4343098 0.8241340 0.9929936 0.7634786
#> 2 1.052458 4.766634 1.045084 0.9084866 0.4343098 0.8241340 0.9929936 0.7634786
#> 3 1.052145 4.766933 1.044774 0.9081824 0.4342513 0.8241321 0.9929944 0.7634924
#> 4 1.052145 4.766933 1.044774 0.9081824 0.4342513 0.8241321 0.9929944 0.7634924
#> 5 1.052563 4.770830 1.045188 0.9080272 0.4307731 0.8241336 0.9929934 0.7634890
#> 6 1.052563 4.770830 1.045188 0.9080272 0.4307731 0.8241336 0.9929934 0.7634890
#>   tad dosenum
#> 1 0.5       1
#> 2 0.5       1
#> 3 0.5       1
#> 4 0.5       1
#> 5 0.5       1
#> 6 0.5       1

# input data (seen by saem) at time = 0.5
print(fit$dataSav %>% dplyr::filter(TIME==0.5))
#> Warning: Undefined 'dvid' integer values in data: 1, 2
#>   ID TIME EVID AMT II  DV CMT nlmixrRowNums
#> 1  1  0.5    0  NA  0 3.8   3             2
#> 2  1  0.5    0  NA  0 0.9   4            15
#> 3  2  0.5    0  NA  0 3.4   3            28
#> 4  2  0.5    0  NA  0 0.7   4            41
#> 5  3  0.5    0  NA  0 2.8   3            54
#> 6  3  0.5    0  NA  0 1.1   4            67

Created on 2023-09-11 with reprex v2.0.2

It is likely in the table output step. I will transfer this to nlmixr2est to see if I can investigate this a bit further.

mattfidler commented 1 year ago

This seems to happen in some edge case where the values for different end-points are tied and they are not synced appropriately. https://github.com/nlmixr2/rxode2/issues/581 has a more precise reproducible example that shows this. While all the times are tied for some reason it only applies to the 0.5 time point.

In the mean-time, by simply offsetting the TIME by a very small amount for one of the compartments will make the fit perform a bit better, and the time-point at approximately 0.5 shows the right compartment:

library(nlmixr2est)
#> Warning: package 'nlmixr2est' was built under R version 4.3.1
#> Loading required package: nlmixr2data
#> Warning: package 'nlmixr2data' was built under R version 4.3.1

one.0abs.PM <- function() {
  ini({
    tkm <- log(c(0,1)) 
    tcl <- log(c(0,2)) 
    tclm <- log(c(0,10))  
    tv <- log(c(0,10))  
    tdur0 <- log(c(0,1))  
    eta.cl ~ 1
    eta.v ~ 1
    eta.km ~ 1 
    eta.clm ~ 1
    eta.dur0 ~ 1
    prop.err <- 1 
    prop.err2 <- 1 
  })

  model({
    km <- exp(tkm + eta.km)
    cl <- exp(tcl + eta.cl)
    clm <- exp(tclm + eta.clm)
    v <- exp(tv + eta.v)
    dur0    <- exp(tdur0+eta.dur0)
    d/dt(center) = - cl / v * center - km * center
    dur(center)=dur0
    d/dt(meta) = km * center - clm/ v * meta
    center(0)= 0
    meta(0)= 0
    cp = center / v
    cm = meta / v
    cp ~ prop(prop.err)   
    cm ~  prop(prop.err2) 
  })
}

dat <- read.csv("c:/tmp/0ABS_PM_Report3.csv",header = T,encoding = "UTF-8", na.strings=c("."))
dat$TIME[dat$CMT==4] <- dat$TIME[dat$CMT==4]+0.0001

fit <- nlmixr2(one.0abs.PM, dat, list(print=0), est="saem")
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in saem model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> Calculating covariance matrix
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of saem model...
#> ✔ done
#> → finding duplicate expressions in saem predOnly model 0...
#> → optimizing duplicate expressions in saem predOnly model 0...
#> → finding duplicate expressions in saem predOnly model 1...
#> → optimizing duplicate expressions in saem predOnly model 1...
#> → finding duplicate expressions in saem predOnly model 2...
#> → optimizing duplicate expressions in saem predOnly model 2...
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 4128
#> → compress phiM in nlmixr2 object, save 77096
#> → compress parHist in nlmixr2 object, save 20704
#> → compress saem0 in nlmixr2 object, save 38696

print(fit)
#> ── nlmixr² SAEM OBJF by FOCEi approximation ──
#> 
#>  Gaussian/Laplacian Likelihoods: AIC() or $objf etc. 
#>  FOCEi CWRES & Likelihoods: addCwres() 
#> 
#> ── Time (sec $time): ──
#> 
#>         setup saem table compress other
#> elapsed 0.002 6.94  0.03     0.03 4.088
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>              Est.    SE %RSE Back-transformed(95%CI)  BSV(CV%) Shrink(SD)%
#> tkm       -0.0368 0.224  611      0.964 (0.621, 1.5)   0.00970      94.6% 
#> tcl         -1.34 0.967 72.2    0.262 (0.0393, 1.74)     0.124      95.9% 
#> tclm       0.0805 0.236  294      1.08 (0.682, 1.72)   0.00370      94.9% 
#> tv          0.163 0.188  115       1.18 (0.814, 1.7) 0.0000599      97.1% 
#> tdur0     -0.0984 0.144  147      0.906 (0.683, 1.2)    0.0335      98.4% 
#> prop.err    0.594                              0.594                      
#> prop.err2   0.332                              0.332                      
#>  
#>   Covariance Type ($covMethod): linFim
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Information about run found ($runInfo):
#>    • Undefined 'dvid' integer values in data: 1, 2 
#>    • package 'rxode2' was built under R version 4.3.1 
#>   Censoring ($censInformation): No censoring
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 75 × 25
#>   ID      TIME CMT      DV         PRED      RES   IPRED     IRES  IWRES  eta.cl
#>   <fct>  <dbl> <fct> <dbl>        <dbl>    <dbl>   <dbl>    <dbl>  <dbl>   <dbl>
#> 1 1     0.0001 cm      0   0.0000000451 -4.51e-8 4.51e-8 -4.51e-8 -3.01  5.08e-5
#> 2 1     0.5    cp      3.8 3.54          2.65e-1 3.54e+0  2.65e-1  0.126 5.08e-5
#> 3 1     0.500  cm      0.9 0.802         9.82e-2 8.02e-1  9.82e-2  0.369 5.08e-5
#> # ℹ 72 more rows
#> # ℹ 15 more variables: eta.v <dbl>, eta.km <dbl>, eta.clm <dbl>,
#> #   eta.dur0 <dbl>, cp <dbl>, cm <dbl>, center <dbl>, meta <dbl>, km <dbl>,
#> #   cl <dbl>, clm <dbl>, v <dbl>, dur0 <dbl>, tad <dbl>, dosenum <dbl>

print(fit %>% as.data.frame() %>% dplyr::mutate(TIME= round(TIME,3)) %>%
        dplyr::filter(TIME==0.5) %>%
        dplyr::select(ID, TIME, CMT, IPRED, DV))
#>   ID TIME CMT     IPRED  DV
#> 1  1  0.5  cp 3.5352614 3.8
#> 2  1  0.5  cm 0.8017849 0.9
#> 3  2  0.5  cp 3.5352739 3.4
#> 4  2  0.5  cm 0.8017826 0.7
#> 5  3  0.5  cp 3.5352973 2.8
#> 6  3  0.5  cm 0.8017842 1.1

Created on 2023-09-11 with reprex v2.0.2

While it isn't an ideal fix, it allows you to use the CRAN version of nlmixr2 and still get reasonable results.

I am hoping to have a CRAN release of the nlmixr2 family of packages soon once a few things get fixed (including this particular issue now).

Thank you for raising this issue so we can fix this case of 2 tied endpoints not synced appropriately.

mattfidler commented 1 year ago

Since this is no longer the issue (and I have a different issue tagged for this problem). I will close this for now and tag you in the rxode2 issue so you can track progress.

mattfidler commented 1 year ago

@Lee21A this has been updated in development.

For the ODE system you are using it should work just fine.

This means with the new r-universe build you should be able to run without issue.

This would update everything needed; You will likely need to wait for an hour or so before the build is complete.

install.packages(c("dparser", "nlmixr2data", "lotri", "rxode2ll",
                   "rxode2parse", "rxode2random", "rxode2et",
                   "rxode2", "nlmixr2est", "nlmixr2extra", "nlmixr2plot",
                   "nlmixr2"),
                 repos = c('https://nlmixr2.r-universe.dev',
                           'https://cloud.r-project.org'))
Lee21A commented 1 year ago

I really appreciate your quick answer. It works!

mattfidler commented 1 year ago

Great!