Help please:different estimation results in using time data or addl data.

[nmarusan]

Hi nlmixr team, I was trying to fit a 1-compartment model using the two datasets. I have tried:

1. ADDL and II-based dataset
2. only time-based dataset

The first dataset seems to work fine, but somehow, the second dataset cannot be estimated well. data2 is a conversion of data1. Why would the estimation results be different if I did not use ADDl and II, but only time?

I have struggled with this problem forever. Would be so appreciated if I could get some help, please. I have attached the dataset as well.

#### Setup 
library(tidyverse)
library(nlmixr)
library(RxODE)

#### Data import
data1 <- read_csv("data1.csv")
data2 <- read_csv("data2.csv")

#### Model definition
one_add <- function(){
    # initial estimates
    ini({
        # Fixed effects
        tCl <- log(c(0, 3))  # log Cl (L/hr)
        tV  <- log(c(0, 10)) # log V (L)
        # Random effects
        eta.Cl ~ 0.1 
        eta.V ~ 0.1 
        # Residual error 
        add.err <- 12  # additive error
    })
    # Model
    model({
        # Model parameters
        Cl <- exp(tCl + eta.Cl)
        V  <- exp(tV + eta.V)
        # Elimination rate constant 
        Kel = Cl / V
        # ODEs
        d/dt(centr)  = -Kel * centr;
        # Concentration
        C = centr / V;
        # Additive error (estimated by the parameter add.err)
        C ~ add(add.err)
    })
}

#### Estimation
This is the first dataset (ADDL, II);

first_fit <- 
nlmixr(
    object=one_add, 
    data=data1,  
    est="foce",
    control=foceiControl(print=0,outerOpt="bobyqa"),
    table=list(cwres=TRUE, npde=TRUE)
)

This is the second dataset (TIME);

second_fit <- 
nlmixr(
    object=one_add, 
    data=data2,  
    est="foce",
    control=foceiControl(print=0,outerOpt="bobyqa"),
    table=list(cwres=TRUE, npde=TRUE)
)

#### Environment
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.5 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] RxODE_1.1.5         GGally_2.1.2.9000   plotly_4.10.0      
 [4] huxtable_5.4.0      broom.mixed_0.2.9.4 latex2exp_0.9.4    
 [7] skimr_2.1.4         nlmixr_2.0.7        forcats_0.5.1      
[10] stringr_1.4.0       dplyr_1.0.9         purrr_0.3.4        
[13] readr_2.1.2         tidyr_1.2.0         tibble_3.1.6       
[16] ggplot2_3.3.6       tidyverse_1.3.1     httr_1.4.2         
[19] R.utils_2.11.0      R.oo_1.24.0         R.methodsS3_1.8.1  
[22] googledrive_2.0.0  

loaded via a namespace (and not attached):
 [1] minqa_1.2.4         colorspace_2.0-3    ellipsis_0.3.2     
 [4] lbfgsb3c_2020-3.2   IRdisplay_1.1       base64enc_0.1-3    
 [7] fs_1.5.2            rstudioapi_0.13     listenv_0.8.0      
[10] furrr_0.3.0         bit64_4.0.5         fansi_1.0.3        
[13] lubridate_1.8.0     xml2_1.3.3          codetools_0.2-18   
[16] splines_4.2.0       cachem_1.0.6        knitr_1.39         
[19] IRkernel_1.3        jsonlite_1.8.0      broom_0.8.0        
[22] dbplyr_2.1.1        compiler_4.2.0      backports_1.4.1    
[25] Matrix_1.4-1        assertthat_0.2.1    fastmap_1.1.0      
[28] lazyeval_0.2.2      gargle_1.2.0        cli_3.3.0          
[31] htmltools_0.5.2     tools_4.2.0         gtable_0.3.0       
[34] glue_1.6.2          Rcpp_1.0.8.3        lotri_0.4.1        
[37] cellranger_1.1.0    dparser_1.3.1-5     vctrs_0.4.1        
[40] n1qn1_6.0.1-10      nlme_3.1-157        xfun_0.30          
[43] globals_0.15.0      rvest_1.0.2         lifecycle_1.0.1    
[46] sys_3.4             symengine_0.2.1     future_1.25.0      
[49] scales_1.2.0        vroom_1.5.7         hms_1.1.1          
[52] rex_1.2.1           parallel_4.2.0      RColorBrewer_1.1-3 
[55] qs_0.25.3           curl_4.3.2          memoise_2.0.1      
[58] reshape_0.8.9       stringi_1.7.6       checkmate_2.1.0    
[61] repr_1.1.4          rlang_1.0.2         pkgconfig_2.0.3    
[64] PreciseSums_0.5     evaluate_0.15       lattice_0.20-45    
[67] htmlwidgets_1.5.4   bit_4.0.4           tidyselect_1.1.2   
[70] parallelly_1.31.1   plyr_1.8.7          magrittr_2.0.3     
[73] R6_2.5.1            generics_0.1.2      pbdZMQ_0.3-7       
[76] DBI_1.1.2           pillar_1.7.0        haven_2.5.0        
[79] withr_2.5.0         modelr_0.1.8        crayon_1.5.1       
[82] uuid_1.1-0          utf8_1.2.2          RApiSerialize_0.1.0
[85] tzdb_0.3.0          grid_4.2.0          readxl_1.4.0       
[88] data.table_1.14.2   reprex_2.0.1        digest_0.6.29      
[91] brew_1.0-7          RcppParallel_5.1.5  munsell_0.5.0      
[94] stringfish_0.15.7   viridisLite_0.4.0

dataset.zip

[mattfidler]

According to rxode2 these two datasets are not equivalent

> library(rxode2)

one_add <- function(){
  # initial estimates
  ini({
    # Fixed effects
    tCl <- log(c(0, 3))  # log Cl (L/hr)
    tV  <- log(c(0, 10)) # log V (L)
    # Random effects
    eta.Cl ~ 0.1 
    eta.V ~ 0.1 
    # Residual error 
    add.err <- 12  # additive error
  })
  # Model
  model({
    # Model parameters
    Cl <- exp(tCl + eta.Cl)
    V  <- exp(tV + eta.V)
    # Elimination rate constant 
    Kel = Cl / V
    # ODEs
    d/dt(centr)  = -Kel * centr;
    # Concentration
    C = centr / V;
    # Additive error (estimated by the parameter add.err)
    C ~ add(add.err)
  })
}

f <- one_add()

d1 <- read.csv("c:/tmp/data1.csv")

d2 <- read.csv("c:/tmp/data2.csv")

rx <-f$simulationModel

d1 <- etTrans(d1, rx)

d2 <- etTrans(d2, rx)

library(waldo)

waldo::compare(d1, d2)
> > . + > > > > > > > > > > > > > > > > > `attr(class(old), '.rxode2.lst')$ndose`: 1200
`attr(class(new), '.rxode2.lst')$ndose`:  600

`attr(class(old), '.rxode2.lst')$nobs`: 104
`attr(class(new), '.rxode2.lst')$nobs`: 156

`attr(class(old), '.rxode2.lst')$allBolus`: FALSE
`attr(class(new), '.rxode2.lst')$allBolus`: TRUE 

`attr(class(old), '.rxode2.lst')$mxCmt`: 1
`attr(class(new), '.rxode2.lst')$mxCmt`: 2

attr(attr(class(old), '.rxode2.lst')$keepL, 'row.names')[606:1008] vs attr(attr(class(new), '.rxode2.lst')$keepL, 'row.names')[606:608]
  606
  607
  608
- 609
- 610
- 611
- 612
- 613
- 614
- 615
and 393 more ...

      attr(old, 'row.names') | attr(new, 'row.names')                 
[606] 606                    | 606                    [606]           
[607] 607                    | 607                    [607]           
[608] 608                    | 608                    [608]           
[609] 609                    -                                        
[610] 610                    -                                        
[611] 611                    -                                        
[612] 612                    -                                        
[613] 613                    -                                        
[614] 614                    -                                        
[615] 615                    -                                        
  ... ...                      ...                    and 393 more ...

old vs new
              ID TIME  EVID    AMT II     DV
- old[1, ]    1   0.0 10101  127.5  6     NA
+ new[1, ]    1   0.0   101  255.0  0     NA
- old[2, ]    1   2.0 10101 -127.5  0     NA
+ new[2, ]    1   0.0     0     NA  0   0.00
- old[3, ]    1   6.0 10101  127.5  0     NA
+ new[3, ]    1   6.0   101  255.0  0     NA
- old[4, ]    1   8.0 10101 -127.5  0     NA
+ new[4, ]    1  12.0   101  255.0  0     NA
- old[5, ]    1  12.0 10101  127.5  0     NA
+ new[5, ]    1  18.0   101  255.0  0     NA
- old[6, ]    1  14.0 10101 -127.5  0     NA
+ new[6, ]    1  24.0   101  255.0  0     NA
- old[7, ]    1  18.0 10101  127.5  0     NA
+ new[7, ]    1  30.0   101  255.0  0     NA
- old[8, ]    1  20.0 10101 -127.5  0     NA
+ new[8, ]    1  36.0   101  255.0  0     NA
- old[9, ]    1  24.0 10101  127.5  0     NA
+ new[9, ]    1  42.0   101  255.0  0     NA
- old[10, ]   1  26.0 10101 -127.5  0     NA
+ new[10, ]   1  48.0   101  255.0  0     NA
and 998 more ...

`old$ID[10:26]`: "1 " "1 " "1 " "1 " "1 " "1 " "1 " "1 " "1 " "1 " and 7 more...
`new$ID[10:15]`: "1 " "1 " "1 "                                    ...          

`old$ID[33:50]`: "2 " "2 " "2 " "2 " "2 " "2 " "2 " "2 " "4 " "4 " and 8 more...
`new$ID[22:27]`: "2 " "2 " "2 "                                    ...          

`old$ID[57:67]`: "4 " "4 " "4 " "4 " "6 " "6 " "6 " "6 " "6 " "6 " and 1 more...
`new$ID[34:39]`: "4 " "4 " "4 "                          "6 " "6 " ...          

And 87 more differences ...
> 

[mattfidler]

This is what is see:

`old$TIME`: 0 2 6 8 12 14 18 20 24 26 and 10 more...
`new$TIME`: 0 0 6   12    18    24               ...

     old$EVID | new$EVID                
 [1] 10101    - 101      [1]            
 [2] 10101    - 0        [2]            
 [3] 10101    - 101      [3]            
 [4] 10101    - 101      [4]            
 [5] 10101    - 101      [5]            
 [6] 10101    - 101      [6]            
 [7] 10101    - 101      [7]            
 [8] 10101    - 101      [8]            
 [9] 10101    - 101      [9]            
[10] 10101    - 101      [10]           
 ... ...        ...      and 10 more ...

Which means d2 is assuming bolus dosing and d1 assumes infusions.

If you are using RxODE or rxode2 with RATE you cannot use the classic evid information

[nmarusan]

Thank Dr. Fidler for responding to my personal question. I understand how the dataset is loaded through nlmixr.