search

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

simultaneous fit of data in CMT=1 and CMT=2 #519

Closed sliao999 closed 3 years ago

sliao999 commented 3 years ago

Are there any example of simultaneous fit of data? and use of log-transformed data I am really interested to learn how to use nlmixr in pk modeling.

billdenney commented 3 years ago

For simultaneous fitting, I would recommend starting from https://nlmixrdevelopment.github.io/nlmixr/articles/multiple-endpoints.html and also specifically looking at how to structure your dataset (specifically looking at the use of CMT and DVID) on the same page.

For log-transform, I'm not sure what the best practices are. But there is some info here about what you may be able to try: https://github.com/nlmixrdevelopment/nlmixr/blob/master/vignettes/residualErrors.csv

mattfidler commented 3 years ago

Thanks Bill. I would agree with the suggestion by Bill, the multiple endpoints vignettes is a good place to start. I think that the best way to deal with log-normally distributed data is to use the lnorm distribution on untransformed data, which is a bit different than the typical NONMEM approach. This is because:

Here is an example using the warfarin dataset:

library(nlmixr)

pk.turnover.emax3 <- function() {
  ini({
    tktr <- log(1)
    tka <- log(1)
    tcl <- log(0.1)
    tv <- log(10)
    ##
    eta.ktr ~ 1
    eta.ka ~ 1
    eta.cl ~ 2
    eta.v ~ 1
    lnorm.sd <- 0.1
    ##
    temax <- logit(0.8)
    tec50 <- log(0.5)
    tkout <- log(0.05)
    te0 <- log(100)
    ##
    eta.emax ~ .5
    eta.ec50  ~ .5
    eta.kout ~ .5
    eta.e0 ~ .5
    ##
    pdadd.err <- 10
  })
  model({
    ktr <- exp(tktr + eta.ktr)
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    emax = expit(temax+eta.emax)
    ec50 =  exp(tec50 + eta.ec50)
    kout = exp(tkout + eta.kout)
    e0 = exp(te0 + eta.e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    effect(0) = e0
    kin = e0*kout
    ##
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center
    d/dt(effect) = kin*PD -kout*effect
    ##
    cp = center / v
    cp ~ lnorm(lnorm.sd)
    effect ~ add(pdadd.err) | pca
  })
}

fit.TOS <- nlmixr(pk.turnover.emax3, warfarin, "saem", control=list(print=0),
                  table=list(cwres=TRUE, npde=TRUE))
#> ℹ parameter labels from comments will be replaced by 'label()'
#> → generate SAEM model
#> ✔ done
#> → creating RxODE include directory
#> → getting R compile options
#> → precompiling headers
#> ✔ done
#> RxODE 1.0.8.1 using 4 threads (see ?getRxThreads)
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → creating full model...
#> → pruning branches (`if`/`else`)...
#> ✔ done
#> → loading into symengine environment...
#> ✔ done
#> → calculate jacobian
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate sensitivities
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(f)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → calculate ∂(R²)/∂(η)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling inner model...
#> ✔ done
#> → compiling EBE model...
#> ✔ done
#> Needed Covariates:
#> [1] "CMT"
#> Calculating residuals/tables
#> Compiling NPDE model...done
#> done

print(fit.TOS)
#> ── nlmixr SAEM(ODE); OBJF by FOCEi approximation fit ───────────────────────────
#> 
#>           OBJF      AIC      BIC Log-likelihood Condition Number
#> FOCEi 3546.027 4469.722 4544.962      -2216.861         708191.9
#> 
#> ── Time (sec $time): ───────────────────────────────────────────────────────────
#> 
#>           saem  setup table  cwres covariance    other
#> elapsed 30.235 11.073 1.289 12.418      0.027 1.870001
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ───────────────────────────
#> 
#>             Est.     SE  %RSE    Back-transformed(95%CI) BSV(CV% or SD)
#> tktr       -1.62   6.63   410 0.198 (4.48e-07, 8.76e+04)           30.1
#> tka        -1.67   6.54   391 0.188 (5.05e-07, 6.96e+04)           34.2
#> tcl         -1.9  0.148  7.77         0.149 (0.112, 0.2)           32.5
#> tv           1.9  0.161  8.45          6.69 (4.88, 9.17)           15.3
#> lnorm.sd    2.08                                    2.08               
#> temax       2.61  0.306  11.7       0.932 (0.882, 0.961)          0.110
#> tec50     -0.292  0.212  72.8        0.747 (0.493, 1.13)           38.3
#> tkout      -2.82 0.0494  1.75     0.0594 (0.054, 0.0655)           7.01
#> te0         4.57 0.0113 0.246            96.1 (94, 98.2)           5.25
#> pdadd.err   3.34                                    3.34               
#>           Shrink(SD)%
#> tktr           74.8% 
#> tka            73.4% 
#> tcl            13.4% 
#> tv             60.7% 
#> lnorm.sd             
#> temax          80.8% 
#> tec50          23.4% 
#> tkout          46.1% 
#> te0            14.1% 
#> pdadd.err            
#>  
#>   Covariance Type ($covMethod): linFim
#>   Some strong fixed parameter correlations exist ($cor) :
#>        cor:tka,tktr    cor:tcl,tktr     cor:tv,tktr  cor:temax,tktr  cor:tec50,tktr 
#>         -0.998          0.0845         -0.0161         -0.0773        -0.00517  
#>  cor:tkout,tktr    cor:te0,tktr     cor:tcl,tka      cor:tv,tka   cor:temax,tka 
#>          0.190         -0.0132         -0.0735          0.0343           0.101  
#>   cor:tec50,tka   cor:tkout,tka     cor:te0,tka      cor:tv,tcl   cor:temax,tcl 
#>        0.00331          -0.231          0.0156           0.786          0.0594  
#>   cor:tec50,tcl   cor:tkout,tcl     cor:te0,tcl    cor:temax,tv    cor:tec50,tv 
#>         -0.589          -0.143         0.00746           0.397          -0.275  
#>    cor:tkout,tv      cor:te0,tv cor:tec50,temax cor:tkout,temax   cor:te0,temax 
#>         -0.180          0.0365           0.590          -0.537          0.0215  
#> cor:tkout,tec50   cor:te0,tec50   cor:te0,tkout 
#>        -0.0345         -0.0230          0.0385  
#>  
#> 
#>   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 
#> 
#> ── Fit Data (object is a modified tibble): ─────────────────────────────────────
#> # A tibble: 483 x 42
#>   ID     TIME CMT      DV    EPRED      ERES   NPDE    NPD   PRED     RES   WRES
#>   <fct> <dbl> <fct> <dbl>    <dbl>     <dbl>  <dbl>  <dbl>  <dbl>   <dbl>  <dbl>
#> 1 1       0.5 cp      0    1.03e-7  -8.85e-8 -1.16  -1.16  0.0648 -0.0648 -7.19 
#> 2 1       1   cp      1.9  1.47e-6   1.90e+0  2.13   2.13  0.242   1.66    0.969
#> 3 1       2   cp      3.3  7.54e-5   3.30e+0  0.854  0.854 0.849   2.45    0.641
#> # … with 480 more rows, and 31 more variables: IPRED <dbl>, IRES <dbl>,
#> #   IWRES <dbl>, CPRED <dbl>, CRES <dbl>, CWRES <dbl>, eta.ktr <dbl>,
#> #   eta.ka <dbl>, eta.cl <dbl>, eta.v <dbl>, eta.emax <dbl>, eta.ec50 <dbl>,
#> #   eta.kout <dbl>, eta.e0 <dbl>, depot <dbl>, gut <dbl>, center <dbl>,
#> #   effect <dbl>, ktr <dbl>, ka <dbl>, cl <dbl>, v <dbl>, emax <dbl>,
#> #   ec50 <dbl>, kout <dbl>, e0 <dbl>, DCP <dbl>, PD <dbl>, kin <dbl>,
#> #   tad <dbl>, dosenum <dbl>

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

While you can log transform the data, it gives the same answers as shown in the following example:

https://nlmixrdevelopment.github.io/nlmixr/articles/mavoglurant.html#lognoral-1

sliao999 commented 3 years ago

Thanks, Bill, and Mat, This is very helpful.

Sam