nlmixrdevelopment / nlmixr

nlmixr: an R package for population PKPD modeling
https://nlmixrdevelopment.github.io/nlmixr/
GNU General Public License v2.0
114 stars 45 forks source link

Names error setting fixed with saem #502

Closed billdenney closed 3 years ago

billdenney commented 3 years ago

What I reported in #501 was actually two separate issues and not a single issue. I'm separating out the second issue here, and here is a reprex for it.

library(nlmixr)
d_blind <-
  structure(list(
    ID = c(
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
      2, 2, 2, 2
    ), CMT = c(
      "TG", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "TG", "CENTRAL", "TG", "TG", "TG"
    ), DVID = c(
      "TG",
      "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp",
      "TG", "cp", "TG", "TG", "TG"
    ), TIME = c(
      -672, 0, 1, 4, 8, 24,
      48, 168, 336, 672, 1008, 1680, -24, 2, 48, 336, 1008
    ), DV = c(
      1.23857621002093,
      1.61552201665003, 1.0918010859138, 1.55870992424631, 0.695699182924043,
      1.130515846091, 0.42153806881519, 1.63170301959007, 0.69480995799257,
      0.274105555786868, 0.474369348047014, 2.51318371567032, 2.11711507953987,
      0.0813858368249748, 0.0477902280089111, 1.00026583587068, 0.674341247526479
    ), MDV = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
    WEIGHT_BL = c(
      73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4,
      73.4, 73.4, 73.4, 73.4, 95.5, 95.5, 95.5, 95.5, 95.5
    )
  ), row.names = c(
    NA,
    -17L
  ), class = c("tbl_df", "tbl", "data.frame"))

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)
  })
}

fit <-
  nlmixr(
    object=nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2,
    data=d_blind,
    est="saem",
    control=list(print=1), 
    table=list(cwres=TRUE, npde=TRUE)
  )
#> > generate SAEM model
#> v done
#> Error in names(inits$theta)[fixid] <- "FIXED": 'names' attribute [21] must be the same length as the vector [20]

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

mattfidler commented 3 years ago

Of note, the option singleOde=FALSE does not work either:

library(nlmixr)

d_mask <-
  structure(list(
    ID = c(
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
      2, 2, 2, 2
    ), CMT = c(
      "TG", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "TG", "CENTRAL", "TG", "TG", "TG"
    ), DVID = c(
      "TG",
      "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp",
      "TG", "cp", "TG", "TG", "TG"
    ), TIME = c(
      -672, 0, 1, 4, 8, 24,
      48, 168, 336, 672, 1008, 1680, -24, 2, 48, 336, 1008
    ), DV = c(
      1.23857621002093,
      1.61552201665003, 1.0918010859138, 1.55870992424631, 0.695699182924043,
      1.130515846091, 0.42153806881519, 1.63170301959007, 0.69480995799257,
      0.274105555786868, 0.474369348047014, 2.51318371567032, 2.11711507953987,
      0.0813858368249748, 0.0477902280089111, 1.00026583587068, 0.674341247526479
    ), MDV = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
    WEIGHT_BL = c(
      73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4,
      73.4, 73.4, 73.4, 73.4, 95.5, 95.5, 95.5, 95.5, 95.5
    )
  ), row.names = c(
    NA,
    -17L
  ), class = c("data.frame"))

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)
  })
}

fit <- nlmixr(
  object=nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2,
  data=d_mask,
  est="saem",
  control=list(print=1, nEm = 3, nBurn = 3, singleOde=FALSE)
)
#> → generate SAEM model
#> ✔ done
#> Error in names(inits$theta)[fixid] <- "FIXED": 'names' attribute [21] must be the same length as the vector [20]

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

mattfidler commented 3 years ago

Simply fixing "by hand" may have another insight but it doesn't work either:

library(nlmixr)

d_mask <-
  structure(list(
    ID = c(
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
      2, 2, 2, 2
    ), CMT = c(
      "TG", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "TG", "CENTRAL", "TG", "TG", "TG"
    ), DVID = c(
      "TG",
      "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp",
      "TG", "cp", "TG", "TG", "TG"
    ), TIME = c(
      -672, 0, 1, 4, 8, 24,
      48, 168, 336, 672, 1008, 1680, -24, 2, 48, 336, 1008
    ), DV = c(
      1.23857621002093,
      1.61552201665003, 1.0918010859138, 1.55870992424631, 0.695699182924043,
      1.130515846091, 0.42153806881519, 1.63170301959007, 0.69480995799257,
      0.274105555786868, 0.474369348047014, 2.51318371567032, 2.11711507953987,
      0.0813858368249748, 0.0477902280089111, 1.00026583587068, 0.674341247526479
    ), MDV = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
    WEIGHT_BL = c(
      73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4,
      73.4, 73.4, 73.4, 73.4, 95.5, 95.5, 95.5, 95.5, 95.5
    )
  ), row.names = c(
    NA,
    -17L
  ), class = c("data.frame"))

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)
    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)
    e_wt_cl <- 1
    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)
  })
}

fit <- nlmixr(
  object=nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2,
  data=d_mask,
  est="saem",
  control=list(print=1, nEm = 3, nBurn = 3, singleOde=FALSE)
)
#> Error in configsaem(model = model, data = dat, inits = inits, mcmc = .mcmc, : nphi and covstruct dim mismatch

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

mattfidler commented 3 years ago

Nor is using saem's prefered mu-referencing for covariates (math for calculation is different).

library(nlmixr)

d_mask <-
  structure(list(
    ID = c(
      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2,
      2, 2, 2, 2
    ), CMT = c(
      "TG", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL", "CENTRAL",
      "CENTRAL", "TG", "CENTRAL", "TG", "TG", "TG"
    ), DVID = c(
      "TG",
      "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp", "cp",
      "TG", "cp", "TG", "TG", "TG"
    ), TIME = c(
      -672, 0, 1, 4, 8, 24,
      48, 168, 336, 672, 1008, 1680, -24, 2, 48, 336, 1008
    ), DV = c(
      1.23857621002093,
      1.61552201665003, 1.0918010859138, 1.55870992424631, 0.695699182924043,
      1.130515846091, 0.42153806881519, 1.63170301959007, 0.69480995799257,
      0.274105555786868, 0.474369348047014, 2.51318371567032, 2.11711507953987,
      0.0813858368249748, 0.0477902280089111, 1.00026583587068, 0.674341247526479
    ), MDV = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
    WEIGHT_BL = c(
      73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4, 73.4,
      73.4, 73.4, 73.4, 73.4, 95.5, 95.5, 95.5, 95.5, 95.5
    )
  ), row.names = c(
    NA,
    -17L
  ), class = c("data.frame"))

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

d_mask <- d_mask %>% mutate(CWT=(WEIGHT_BL/85))

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)
    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
    e_wt_cl <- 1

  })
  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 + e_wt_cl * CWT)
    cl_high <- exp(tcl_high + eta_cl + e_wt_cl * CWT)
    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)
  })
}

fit <- nlmixr(
  object=nlmixr_threecmt_mm_no_add_wtcl_pdtg_kout_delay2,
  data=d_mask,
  est="saem",
  control=list(print=1, nEm = 3, nBurn = 3)
)
#> → generate SAEM model
#> ✔ done
#> Error in configsaem(model = model, data = dat, inits = inits, mcmc = .mcmc, : nphi and covstruct dim mismatch

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

mattfidler commented 3 years ago

This fix breaks #307

mattfidler commented 3 years ago

Fixed on my end. Merging into master.

billdenney commented 3 years ago

Thank you!

mattfidler commented 3 years ago

I'm glad it works for you