Multiple compartment model cannot detect CMT/DVID correctly?

[billdenney]

When trying to run a model, I got the error below:

  multiple compartment models with expressions need to be conditioned by `|`
ie log(cp) ~ add(err) | cmt
The following endpoints need to be corrected: cp

This seems unusual because the same residual error model had worked in a previous version of the model (that was less complex in the math for one of the compartments, but was otherwise the same). In the dataset, CMT and DVID are both used. CMT includes values like "SC", "CENTRAL", and "TG"; while DVID is NA_character_ for CMT="SC" and "cp" for CMT="CENTRAL" and "TG" for CMT="TG".

The model is below (with all numbers changed to 999; those are just for blinding the actual results).

nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay <- function() {
  ini({
    tf_sc <- log(999)
    tf_infilt <- log(999)
    tka_sc <- log(999)
    tka_infilt <- log(999)
    tcl_low <- log(999)
    tcl_high <- log(999)
    tcl_c50 <- log(3000)
    e_wt_cl <- fixed(999)
    tv <- log(999)
    tq1 <- log(999)
    tvp1 <- log(10)
    tq2 <- log(999)
    tvp2 <- log(20)
    eta_cl~999
    eta_v~999
    prop_err <- 999

    tg_bl <- log(999)
    eta_tg_bl~999
    tg_kel <- log(999)
    tg_ec50 <- log(5000)
    tg_emax_kel <- log(2)
    ktr_tg <- log(999)
    prop_err_tg <- 999
  })
  no_deps(model({
    # PK setup
    f_sc <- exp(tf_sc)
    f_infilt <- exp(tf_infilt)
    ka_sc <- exp(tka_sc)
    ka_infilt <- exp(tka_infilt)
    cl_low <- exp(tcl_low + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_high <- exp(tcl_high + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_c50 <- exp(tcl_c50)
    v <- exp(tv + eta_v)
    q1 <- exp(tq1)
    vp1 <- exp(tvp1)
    q2 <- exp(tq2)
    vp2 <- exp(tvp2)

    # PK micro-parameters
    ke_low <- cl_low/v
    ke_high <- cl_high/v
    kc_p1 <- q1/v
    kp1_c <- q1/vp1
    kc_p2 <- q2/v
    kp2_c <- q2/vp2

    # differential equations
    cp <- CENTRAL/v*1e3 # 1e3 is for unit conversion
    ke <- ke_low + (ke_high - ke_low)*cp/(cp + cl_c50)
    d/dt(IVINFILT) =             - ka_infilt * IVINFILT
    d/dt(SC)       = -ka_sc * SC
    d/dt(CENTRAL)  =  ka_sc * SC + ka_infilt * IVINFILT - ke*CENTRAL - kc_p1*CENTRAL + kp1_c*P1 - kc_p2*CENTRAL + kp2_c*P2
    d/dt(P1) =                                                         kc_p1*CENTRAL - kp1_c*P1
    d/dt(P2) =                                                                                    kc_p2*CENTRAL - kp2_c*P2

    f(SC) <- f_sc
    f(IVINFILT) <- f_infilt

    # TG Emax model
    tgbl <- exp(tg_bl + eta_tg_bl)
    kin_tg <- tgbl*exp(tg_kel)
    TG(0) <- tgbl
    ktr_TG <- exp(ktr_tg)
    d/dt(TG_TR) = ktr_tg*cp - ktr_tg*TG_TR
    kout_tg <- exp(tg_kel) + exp(tg_emax_kel)*TG_TR/(TG_TR + exp(tg_ec50))
    d/dt(TG) = kin_tg - kout_tg*TG

    # Residual error models
    cp ~ prop(prop_err)
    TG ~ prop(prop_err_tg)
  }))
}

Edit: the first version of this request indicated that the model had already run, but it was a previous model that just finished running.

[billdenney]

After modifying the model part to be what I think is mathematically identical, I get a different error:

  no_deps(model({
    # PK setup
    f_sc <- exp(tf_sc)
    f_infilt <- exp(tf_infilt)
    ka_sc <- exp(tka_sc)
    ka_infilt <- exp(tka_infilt)
    cl_low <- exp(tcl_low + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_high <- exp(tcl_high + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_c50 <- exp(tcl_c50)
    v <- exp(tv + eta_v)
    q1 <- exp(tq1)
    vp1 <- exp(tvp1)
    q2 <- exp(tq2)
    vp2 <- exp(tvp2)

    # PK micro-parameters
    ke_low <- cl_low/v
    ke_high <- cl_high/v
    kc_p1 <- q1/v
    kp1_c <- q1/vp1
    kc_p2 <- q2/v
    kp2_c <- q2/vp2

    # TG setup
    tgbl <- exp(tg_bl + eta_tg_bl)
    kin_tg <- tgbl*exp(tg_kel)
    ktr_TG <- exp(ktr_tg)
    TG(0) <- tgbl

    # differential equations
    cp <- CENTRAL/v*1e3 # 1e3 is for unit conversion
    ke <- ke_low + (ke_high - ke_low)*cp/(cp + cl_c50)
    kout_tg <- exp(tg_kel) + exp(tg_emax_kel)*TG_TR/(TG_TR + exp(tg_ec50))
    d/dt(IVINFILT) =             - ka_infilt * IVINFILT
    d/dt(SC)       = -ka_sc * SC
    d/dt(CENTRAL)  =  ka_sc * SC + ka_infilt * IVINFILT - ke*CENTRAL - kc_p1*CENTRAL + kp1_c*P1 - kc_p2*CENTRAL + kp2_c*P2
    d/dt(P1) =                                                         kc_p1*CENTRAL - kp1_c*P1
    d/dt(P2) =                                                                                    kc_p2*CENTRAL - kp2_c*P2

    f(SC) <- f_sc
    f(IVINFILT) <- f_infilt

    # TG transit model
    d/dt(TG_TR) = ktr_tg*cp - ktr_tg*TG_TR
    d/dt(TG) = kin_tg - kout_tg*TG

    # Residual error models
    cp ~ prop(prop_err)
    TG ~ prop(prop_err_tg)
  }))

The error now is:

'names' attribute [21] must be the same length as the vector [20]

[billdenney]

The error "'names' attribute [21] must be the same length as the vector [20]" is occurring on line 795 here:

https://github.com/nlmixrdevelopment/nlmixr/blob/d77e429ac54cac303bf6df94c43824641f3013c1/R/nlmixr.R#L791-L795

The number of fixed effects being estimated in the model is 20 (including both "normal" fixed effects and residual error estimated, which I think count as fixed effects). When detected a few lines above, nphi is 18 (I think that is the number of "normal" fixed effects that are not residual error estimates). On the match() in line 793, it fins the value to be 21, and I think that causes the problem, but I'm not certain what exactly this code is trying to do (other than set which values are supposed to be fixed).

[mattfidler]

I can reproduce the first here:

library(nlmixr)

nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay <- function() {
  ini({
    tf_sc <- log(999)
    tf_infilt <- log(999)
    tka_sc <- log(999)
    tka_infilt <- log(999)
    tcl_low <- log(999)
    tcl_high <- log(999)
    tcl_c50 <- log(3000)
    e_wt_cl <- fixed(999)
    tv <- log(999)
    tq1 <- log(999)
    tvp1 <- log(10)
    tq2 <- log(999)
    tvp2 <- log(20)
    eta_cl~999
    eta_v~999
    prop_err <- 999

    tg_bl <- log(999)
    eta_tg_bl~999
    tg_kel <- log(999)
    tg_ec50 <- log(5000)
    tg_emax_kel <- log(2)
    ktr_tg <- log(999)
    prop_err_tg <- 999
  })
  model({
    # PK setup
    f_sc <- exp(tf_sc)
    f_infilt <- exp(tf_infilt)
    ka_sc <- exp(tka_sc)
    ka_infilt <- exp(tka_infilt)
    cl_low <- exp(tcl_low + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_high <- exp(tcl_high + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_c50 <- exp(tcl_c50)
    v <- exp(tv + eta_v)
    q1 <- exp(tq1)
    vp1 <- exp(tvp1)
    q2 <- exp(tq2)
    vp2 <- exp(tvp2)

    # PK micro-parameters
    ke_low <- cl_low/v
    ke_high <- cl_high/v
    kc_p1 <- q1/v
    kp1_c <- q1/vp1
    kc_p2 <- q2/v
    kp2_c <- q2/vp2

    # differential equations
    cp <- CENTRAL/v*1e3 # 1e3 is for unit conversion
    ke <- ke_low + (ke_high - ke_low)*cp/(cp + cl_c50)
    d/dt(IVINFILT) =             - ka_infilt * IVINFILT
    d/dt(SC)       = -ka_sc * SC
    d/dt(CENTRAL)  =  ka_sc * SC + ka_infilt * IVINFILT - ke*CENTRAL - kc_p1*CENTRAL + kp1_c*P1 - kc_p2*CENTRAL + kp2_c*P2
    d/dt(P1) =                                                         kc_p1*CENTRAL - kp1_c*P1
    d/dt(P2) =                                                                                    kc_p2*CENTRAL - kp2_c*P2

    f(SC) <- f_sc
    f(IVINFILT) <- f_infilt

    # TG Emax model
    tgbl <- exp(tg_bl + eta_tg_bl)
    kin_tg <- tgbl*exp(tg_kel)
    TG(0) <- tgbl
    ktr_TG <- exp(ktr_tg)
    d/dt(TG_TR) = ktr_tg*cp - ktr_tg*TG_TR
    kout_tg <- exp(tg_kel) + exp(tg_emax_kel)*TG_TR/(TG_TR + exp(tg_ec50))
    d/dt(TG) = kin_tg - kout_tg*TG

    # Residual error models
    cp ~ prop(prop_err)
    TG ~ prop(prop_err_tg)
  })
}

nlmixr(nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay)
#> Error: multiple compartment models with expressions need to be conditioned by `|`
#> ie log(cp) ~ add(err) | cmt
#> The following endpoints need to be corrected: cp

^{Created on 2021-04-15 by the reprex package (v2.0.0)}

[mattfidler]

The second requires some dummy data to reproduce as it parses ok. It may work with focei

library(nlmixr)

nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2 <- function() {
  ini({
    tf_sc <- log(999)
    tf_infilt <- log(999)
    tka_sc <- log(999)
    tka_infilt <- log(999)
    tcl_low <- log(999)
    tcl_high <- log(999)
    tcl_c50 <- log(3000)
    e_wt_cl <- fixed(999)
    tv <- log(999)
    tq1 <- log(999)
    tvp1 <- log(10)
    tq2 <- log(999)
    tvp2 <- log(20)
    eta_cl~999
    eta_v~999
    prop_err <- 999

    tg_bl <- log(999)
    eta_tg_bl~999
    tg_kel <- log(999)
    tg_ec50 <- log(5000)
    tg_emax_kel <- log(2)
    ktr_tg <- log(999)
    prop_err_tg <- 999
  })
  model({
    # PK setup
    f_sc <- exp(tf_sc)
    f_infilt <- exp(tf_infilt)
    ka_sc <- exp(tka_sc)
    ka_infilt <- exp(tka_infilt)
    cl_low <- exp(tcl_low + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_high <- exp(tcl_high + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
    cl_c50 <- exp(tcl_c50)
    v <- exp(tv + eta_v)
    q1 <- exp(tq1)
    vp1 <- exp(tvp1)
    q2 <- exp(tq2)
    vp2 <- exp(tvp2)

    # PK micro-parameters
    ke_low <- cl_low/v
    ke_high <- cl_high/v
    kc_p1 <- q1/v
    kp1_c <- q1/vp1
    kc_p2 <- q2/v
    kp2_c <- q2/vp2

    # TG setup
    tgbl <- exp(tg_bl + eta_tg_bl)
    kin_tg <- tgbl*exp(tg_kel)
    ktr_TG <- exp(ktr_tg)
    TG(0) <- tgbl

    # differential equations
    cp <- CENTRAL/v*1e3 # 1e3 is for unit conversion
    ke <- ke_low + (ke_high - ke_low)*cp/(cp + cl_c50)
    kout_tg <- exp(tg_kel) + exp(tg_emax_kel)*TG_TR/(TG_TR + exp(tg_ec50))
    d/dt(IVINFILT) =             - ka_infilt * IVINFILT
    d/dt(SC)       = -ka_sc * SC
    d/dt(CENTRAL)  =  ka_sc * SC + ka_infilt * IVINFILT - ke*CENTRAL - kc_p1*CENTRAL + kp1_c*P1 - kc_p2*CENTRAL + kp2_c*P2
    d/dt(P1) =                                                         kc_p1*CENTRAL - kp1_c*P1
    d/dt(P2) =                                                                                    kc_p2*CENTRAL - kp2_c*P2

    f(SC) <- f_sc
    f(IVINFILT) <- f_infilt

    # TG transit model
    d/dt(TG_TR) = ktr_tg*cp - ktr_tg*TG_TR
    d/dt(TG) = kin_tg - kout_tg*TG

    # Residual error models
    cp ~ prop(prop_err)
    TG ~ prop(prop_err_tg)
  })
}

f <- nlmixr(nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2)

f
#> ▂▂ RxODE-based ODE model ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ 
#> ── Initialization: ───────────────────────────────────────────────────────────── 
#> Fixed Effects ($theta): 
#>       tf_sc   tf_infilt      tka_sc  tka_infilt     tcl_low    tcl_high 
#>   6.9067548   6.9067548   6.9067548   6.9067548   6.9067548   6.9067548 
#>     tcl_c50     e_wt_cl          tv         tq1        tvp1         tq2 
#>   8.0063676 999.0000000   6.9067548   6.9067548   2.3025851   6.9067548 
#>        tvp2       tg_bl      tg_kel     tg_ec50 tg_emax_kel      ktr_tg 
#>   2.9957323   6.9067548   6.9067548   8.5171932   0.6931472   6.9067548 
#> 
#> Omega ($omega): 
#>           eta_cl eta_v eta_tg_bl
#> eta_cl       999     0         0
#> eta_v          0   999         0
#> eta_tg_bl      0     0       999
#> 
#>  Covariates or Uninitialized Parameters ($all.covs)
#> [1] "WEIGHT_BL"
#> ── Multiple Endpoint Model ($multipleEndpoint): ──────────────────────────────── 
#> ┌─────────────────────────┬─────────────────────────┬─────────────────────────┐
#> │ variable                │ cmt                     │ dvid*                   │
#> ├─────────────────────────┼─────────────────────────┼─────────────────────────┤
#> │ variable                │ cmt                     │ dvid*                   │
#> ├─────────────────────────┼─────────────────────────┼─────────────────────────┤
#> │ cp ~ …                  │ cmt='cp' or cmt=8       │ dvid='cp' or dvid=1     │
#> ├─────────────────────────┼─────────────────────────┼─────────────────────────┤
#> │ TG ~ …                  │ cmt='TG' or cmt=7       │ dvid='TG' or dvid=2     │
#> └─────────────────────────┴─────────────────────────┴─────────────────────────┘
#>   * If dvids are outside this range, all dvids are re-numered                  
#>   sequentially, ie 1,7, 10 becomes 1,2,3 etc                                   
#> 
#> ── μ-referencing ($muRefTable): ──────────────────────────────────────────────── 
#>                             ┌──────────┬───────────┐
#>                             │ theta    │ eta       │
#>                             ├──────────┼───────────┤
#>                             │ tcl_high │ eta_cl    │
#>                             ├──────────┼───────────┤
#>                             │ tv       │ eta_v     │
#>                             ├──────────┼───────────┤
#>                             │ tg_bl    │ eta_tg_bl │
#>                             └──────────┴───────────┘
#> 
#> ── Model: ────────────────────────────────────────────────────────────────────── 
#>     # PK setup
#>     f_sc <- exp(tf_sc)
#>     f_infilt <- exp(tf_infilt)
#>     ka_sc <- exp(tka_sc)
#>     ka_infilt <- exp(tka_infilt)
#>     cl_low <- exp(tcl_low + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
#>     cl_high <- exp(tcl_high + eta_cl)*(WEIGHT_BL/85)^e_wt_cl
#>     cl_c50 <- exp(tcl_c50)
#>     v <- exp(tv + eta_v)
#>     q1 <- exp(tq1)
#>     vp1 <- exp(tvp1)
#>     q2 <- exp(tq2)
#>     vp2 <- exp(tvp2)
#> 
#>     # PK micro-parameters
#>     ke_low <- cl_low/v
#>     ke_high <- cl_high/v
#>     kc_p1 <- q1/v
#>     kp1_c <- q1/vp1
#>     kc_p2 <- q2/v
#>     kp2_c <- q2/vp2
#> 
#>     # TG setup
#>     tgbl <- exp(tg_bl + eta_tg_bl)
#>     kin_tg <- tgbl*exp(tg_kel)
#>     ktr_TG <- exp(ktr_tg)
#>     TG(0) <- tgbl
#> 
#>     # differential equations
#>     cp <- CENTRAL/v*1e3 # 1e3 is for unit conversion
#>     ke <- ke_low + (ke_high - ke_low)*cp/(cp + cl_c50)
#>     kout_tg <- exp(tg_kel) + exp(tg_emax_kel)*TG_TR/(TG_TR + exp(tg_ec50))
#>     d/dt(IVINFILT) =             - ka_infilt * IVINFILT
#>     d/dt(SC)       = -ka_sc * SC
#>     d/dt(CENTRAL)  =  ka_sc * SC + ka_infilt * IVINFILT - ke*CENTRAL - kc_p1*CENTRAL + kp1_c*P1 - kc_p2*CENTRAL + kp2_c*P2
#>     d/dt(P1) =                                                         kc_p1*CENTRAL - kp1_c*P1
#>     d/dt(P2) =                                                                                    kc_p2*CENTRAL - kp2_c*P2
#> 
#>     f(SC) <- f_sc
#>     f(IVINFILT) <- f_infilt
#> 
#>     # TG transit model
#>     d/dt(TG_TR) = ktr_tg*cp - ktr_tg*TG_TR
#>     d/dt(TG) = kin_tg - kout_tg*TG
#> 
#>     # Residual error models
#>     cp ~ prop(prop_err)
#>     TG ~ prop(prop_err_tg) 
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂

f$saem.fixed
#> [1] 6

[billdenney]

Recognizing that the second error is separate, I made #502 for it.

[mattfidler]

Works for me. I merged into the master

[githubxr4869]

Hi @billdenney,

I'm a new in modelling with nlmirx2 and I think I'm facing with the same issue on "multiple compartment model detect DVID/CMT". But the error showed as "Error : 'dvid'->'cmt' or 'cmt' on observation record on a undefined compartment (use 'cmt()' 'dvid()') id: 24 row: 0". It will disappear if I only insert one prop.err. I'm wondering how you fix the issue. Many thanks.

[billdenney]

@githubxr4869 , the fixes here resolved the issue for me. Can you please open a new issue in the current nlmixr2 repository with a reproducible example?

https://github.com/nlmixr2/nlmixr2/issues