nlmixrdevelopment / nlmixr

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

CMT information in fitted table #471

Closed tnzo12 closed 3 years ago

tnzo12 commented 3 years ago

Sorry to come up with another issue,

I found that the nlmixr fitted table in multiple-endpoints is not generating the column 'CMT' which can inform that each 'DV' is belonged to...

I attached the code below

library(nlmixr)
library(RxODE)

df <- data.frame(
  ID = 123,
  MDV = c(0,1,0),
  CMT = c(2,1,1),
  DV = c(9,NA,15),
  AMT = c(NA,15,NA),
  RATE = c(NA,0,NA),
  ADDL = c(NA,1,NA),
  II = c(NA,NA,NA),
  TIME = c(0,0,1),
  PCA = c(NA,NA,32),
  WT = c(NA,NA,2),
  CRPZERO = 5
)
df

f <- function(){
  ini({
    # thetas
    theta1  <- c(0.387)     # Add err PK
    theta2  <- c(log(50.9))      # Vd
    theta3  <- c(log(3.42))      # Cl
    theta4  <- c(31.2)      # TM50
    theta5  <- c(3.68)      # Hill
    theta6  <- c(1.46)      # K growth
    theta7  <- c(0.187)     # K death
    theta8  <- c(1.52)      # Emax bact
    theta9  <- c(0.304)     # EC50 bact
    theta10 <- c(4.99)     # Gamma bact
    theta11 <- c(0.162)    # Add err PD
    theta12 <- c(0.274)    # Prop err PD
    theta13 <- c(log(0.276))    # base
    theta14 <- c(log(0.0431))   # Kout
    theta15 <- c(1.22)     # EC50/10^3
    theta16 <- c(0.134)    # Emax
    # ETAs
    eta1 ~ c(1.33)                     # Base
    eta2 + eta3 ~ c(0.194,          
                    -0.0159, 0.306) # Vd, Cl
    eta4 + eta5 ~ c(0.521,          
                    -0.435, 0.83)   # Emax, Kout
  })
  model({
    # maturation params
    TM50 <- theta4
    HILL <- theta5
    # central compartmental params
    V <- (WT/70)*exp(theta2 + eta2)
    Cl <- ((WT/70)**0.75) * (PCA**HILL) / ((PCA**HILL)+(TM50**HILL))*exp(theta3 + eta3)
    # bacterial params
    KGROWTH <- theta6
    KDEATH <- theta7
    EMAXBACT <- theta8
    EC50BACT <- theta9
    GAMMABACT <- theta10
    BACINIT <- 10**6
    # crp level params
    BASE <- exp(theta13 + eta1)
    KOUT <- exp(theta14 + eta5)
    KIN <- BASE*KOUT
    EC50 <- theta15*1000
    EMAX <- theta16*(1 + eta4)
    # initial conditions
    bact(0) <- BACINIT
    crp(0) <- CRPZERO
    # algebraic expressions
    K = Cl/V
    cp = centr/V
    # differential equations
    d/dt(centr) = -K*centr
    d/dt(crp)   = KIN + EMAX*bact/(EC50 + bact) - KOUT*crp
    d/dt(bact)  = KGROWTH*bact - KDEATH*bact - EMAXBACT*(centr/V)**GAMMABACT/(EC50BACT**GAMMABACT+(centr/V)**GAMMABACT)*bact
    # error model
    cp ~ add(theta1) | centr
    crp ~ prop(theta12) + add(theta11) | crp
  })
}

fit.s <- nlmixr(
  object = f,
  data = df,
  est='focei',
  control = foceiControl(
    covMethod="r,s",
    interaction = TRUE,
    maxOuterIterations = 0,
    iter.max=0)
  )
fit.s
mattfidler commented 3 years ago

Indeed; It also is missing the covariates:

"CRPZERO" "WT"      "PCA"
mattfidler commented 3 years ago

Should be fixed:

library(nlmixr)
df <- data.frame(
  ID = c(123, 123, 123, 124, 124, 124),
  MDV = c(0,1,0, 0,1,0),
  CMT = c(2,1,1, 2,1,1),
  DV = c(9,NA,15, 10,NA,14),
  AMT = c(NA,15,NA, NA,15,NA),
  RATE = c(NA,0,NA, NA,0,NA),
  ADDL = c(NA,1,NA, NA,1,NA),
  II = c(NA,NA,NA, NA,NA,NA),
  TIME = c(0,0,1, 0,0,1),
  PCA = c(NA,NA,32, NA,NA,32),
  WT = c(NA,NA,2, NA,NA,2),
  CRPZERO = 5
)

f <- function(){
  ini({
    # thetas
    theta1  <- c(0.387)     # Add err PK
    theta2  <- c(log(50.9))      # Vd
    theta3  <- c(log(3.42))      # Cl
    theta4  <- c(31.2)      # TM50
    theta5  <- c(3.68)      # Hill
    theta6  <- c(1.46)      # K growth
    theta7  <- c(0.187)     # K death
    theta8  <- c(1.52)      # Emax bact
    theta9  <- c(0.304)     # EC50 bact
    theta10 <- c(4.99)     # Gamma bact
    theta11 <- c(0.162)    # Add err PD
    theta12 <- c(0.274)    # Prop err PD
    theta13 <- c(log(0.276))    # base
    theta14 <- c(log(0.0431))   # Kout
    theta15 <- c(1.22)     # EC50/10^3
    theta16 <- c(0.134)    # Emax
    # ETAs
    eta1 ~ c(1.33)                     # Base
    eta2 + eta3 ~ c(0.194,
                    -0.0159, 0.306) # Vd, Cl
                    eta4 + eta5 ~ c(0.521,
                                    -0.435, 0.83)   # Emax, Kout
  })
  model({
    # maturation params
    TM50 <- theta4
    HILL <- theta5
    # central compartmental params
    V <- (WT/70)*exp(theta2 + eta2)
    Cl <- ((WT/70)**0.75) * (PCA**HILL) / ((PCA**HILL)+(TM50**HILL))*exp(theta3 + eta3)
    # bacterial params
    KGROWTH <- theta6
    KDEATH <- theta7
    EMAXBACT <- theta8
    EC50BACT <- theta9
    GAMMABACT <- theta10
    BACINIT <- 10**6
    # crp level params
    BASE <- exp(theta13 + eta1)
    KOUT <- exp(theta14 + eta5)
    KIN <- BASE*KOUT
    EC50 <- theta15*1000
    EMAX <- theta16*(1 + eta4)
    # initial conditions
    bact(0) <- BACINIT
    crp(0) <- CRPZERO
    # algebraic expressions
    K = Cl/V
    cp = centr/V
    # differential equations
    d/dt(centr) = -K*centr
    d/dt(crp)   = KIN + EMAX*bact/(EC50 + bact) - KOUT*crp
    d/dt(bact)  = KGROWTH*bact - KDEATH*bact - EMAXBACT*(centr/V)**GAMMABACT/(EC50BACT**GAMMABACT+(centr/V)**GAMMABACT)*bact
    # error model
    cp ~ add(theta1) | centr
    crp ~ prop(theta12) + add(theta11) | crp
  })
}

fit.s <- nlmixr(
  object = f,
  data = df,
  est='focei',
  control = foceiControl(
    covMethod="r,s",
    interaction = TRUE,
    maxOuterIterations = 0,
    iter.max=0)
)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Needed Covariates:
#> [1] "CRPZERO" "WT"      "PCA"     "CMT"
#> RxODE 1.0.2 using 4 threads (see ?getRxThreads)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> Calculating residuals/tables
#> done
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : initial ETAs were nudged; (can control by foceiControl(etaNudge=.,
#> etaNudge2=))
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod = FALSE, : parameter estimate near boundary; covariance not calculated
#>  use 'getVarCov' to calculate anyway

print(fit.s)
#> ── nlmixr FOCEi (outer: nlminb) posthoc estimation ─────────────────────────────
#> 
#>           OBJF      AIC      BIC Log-likelihood
#> FOCEi 32.05525 85.40675 71.29152      -19.70338
#> 
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#> 
#>            setup optimize covariance table    other
#> elapsed 1.018792 0.000131   0.000133 0.009 0.287944
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#> 
#>           Parameter  Est. Back-transformed BSV(CV%) Shrink(SD)%
#> theta1   Add err PK 0.387            0.387                     
#> theta2           Vd  3.93             50.9     46.3      87.8% 
#> theta3           Cl  1.23             3.42     59.8      97.6% 
#> theta4         TM50  31.2             31.2                     
#> theta5         Hill  3.68             3.68                     
#> theta6     K growth  1.46             1.46                     
#> theta7      K death 0.187            0.187                     
#> theta8    Emax bact  1.52             1.52                     
#> theta9    EC50 bact 0.304            0.304                     
#> theta10  Gamma bact  4.99             4.99                     
#> theta11  Add err PD 0.162            0.162                     
#> theta12 Prop err PD 0.274            0.274                     
#> theta13        base -1.29            0.276     167.      100.% 
#> theta14        Kout -3.14           0.0431     114.      100.% 
#> theta15   EC50/10^3  1.22             1.22                     
#> theta16        Emax 0.134            0.134                     
#>  
#> ── BSV Covariance ($omega): ────────────────────────────────────────────────────
#>      eta1    eta2    eta3   eta4   eta5
#> eta1 1.33  0.0000  0.0000  0.000  0.000
#> eta2 0.00  0.1940 -0.0159  0.000  0.000
#> eta3 0.00 -0.0159  0.3060  0.000  0.000
#> eta4 0.00  0.0000  0.0000  0.521 -0.435
#> eta5 0.00  0.0000  0.0000 -0.435  0.830
#> 
#>   Not all variables are μ-referenced.
#>   Can also see BSV Correlation ($omegaR; diagonals=SDs)
#>   Covariance Type ($covMethod): Boundary issue; Get SEs with getVarCov
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Minimization message ($message):  
#>     Posthoc prediction with provided THETAs 
#> 
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 4 x 42
#>   ID     TIME CMT      DV  PRED   RES  WRES IPRED   IRES  IWRES CPRED  CRES
#>   <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl>
#> 1 123       0 crp       9  5     4     2.90   5   4      2.90    5     4   
#> 2 123       1 centr    15  9.47  5.53  1.44  15.0 0.0302 0.0781  8.27  6.73
#> 3 124       0 crp      10  5     5     3.62   5   5      3.62    5     5   
#> # … with 1 more row, and 30 more variables: CWRES <dbl>, eta1 <dbl>,
#> #   eta2 <dbl>, eta3 <dbl>, eta4 <dbl>, eta5 <dbl>, centr <dbl>, crp <dbl>,
#> #   bact <dbl>, TM50 <dbl>, HILL <dbl>, V <dbl>, Cl <dbl>, KGROWTH <dbl>,
#> #   KDEATH <dbl>, EMAXBACT <dbl>, EC50BACT <dbl>, GAMMABACT <dbl>, BASE <dbl>,
#> #   KOUT <dbl>, KIN <dbl>, EC50 <dbl>, EMAX <dbl>, K <dbl>, cp <dbl>,
#> #   tad <dbl>, dosenum <dbl>, WT <dbl>, PCA <dbl>, CRPZERO <dbl>

Created on 2021-01-20 by the reprex package (v0.3.0)

mattfidler commented 3 years ago

Note this is different than your reprex; You only have 1 ID, which is tecnically not allowed;

Let me know if this works for you too.

tnzo12 commented 3 years ago

I'm afraid your code doesn't reproduce the same result as yours even I tried updating the RxODE and nlmixr

I guess the packages has not been pushed may be...?

I will follow up with your progress #472, thank you very much for your efforts!

mattfidler commented 3 years ago

The packages were pushed. Do you only have one id? When I worked with that it didn't work.

tnzo12 commented 3 years ago

Sorry about confusion, After reboot I got CMT column (with 2 IDs). Thank you!

mattfidler commented 3 years ago

This sometimes happens to me too (especially on windows).

mattfidler commented 3 years ago

The compartment information was incorrect for the multiple endpoint vignette.

mattfidler commented 3 years ago

The CMT values are now correct for the multiple endpoint vignette