--- title: "Integration with nlmixr2" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Integration with nlmixr2} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = FALSE ) ``` ## Overview `nlmixr2autoinit` is designed to slot naturally into the `nlmixr2` ecosystem. Once initial estimates are generated with `getPPKinits()`, they can be extracted and used to write a model function with the initial values pre-filled automatically, ready to pass to `nlmixr2` for fitting. --- ## Generating Initial Estimates The starting point is `getPPKinits()`, which analyses the dataset and returns recommended initial estimates: ```{r generate-inits} library(nlmixr2autoinit) inits.out <- getPPKinits(pheno_sd) inits.out$Recommended_initial_estimates # Parameters Methods Values # 1 Ka IV NA # 2 CL Adaptive single-point method 0.0087 # 3 Vd Adaptive single-point method 1.2500 # 4 Vmax Parameter sweeping 1.0572 # 5 Km Parameter sweeping 120.0145 # 6 Vc(2CMPT) Parameter sweeping 1.2500 # 7 Vp(2CMPT) Parameter sweeping 0.1250 # 8 Q(2CMPT) Parameter sweeping 0.0174 # 9 Vc(3CMPT) Parameter sweeping 1.2500 # 10 Vp(3CMPT) Parameter sweeping 0.1250 # 11 Vp2(3CMPT) Parameter sweeping 0.1250 # 12 Q(3CMPT) Parameter sweeping 0.0174 # 13 Q2(3CMPT) Parameter sweeping 0.0174 # 14 Sigma additive Fallback (fixed fraction) 5.0434 # 15 Sigma proportional Fallback (fixed fraction) 0.2000 ``` For `pheno_sd`, the primary structural parameters are `CL = 0.0087`, `Vd = 1.25` and `Sigma additive = 5.0434`. --- ## Writing the Model The values from `Recommended_initial_estimates` can be read directly into the `ini` block. For `pheno_sd`, the one-compartment IV bolus model is: ```{r ini-manual, eval=F} ini({ tcl <- log(0.0087) # typical value of clearance (log scale) tv <- log(1.25) # typical value of volume (log scale) eta.cl + eta.v ~ c(1, 0.01, 1) add.err <- 5.0434 # additive residual variability }) ``` You can also extract the estimates into named variables first, then reference them directly in the `ini` block: ```{r one-cmpt, eval=F} est <- inits.out$Recommended_initial_estimates CL = as.numeric(est$Values[est$Parameters == "CL"]) Vd = as.numeric(est$Values[est$Parameters == "Vd"]) add_err = as.numeric(est$Values[est$Parameters == "Sigma additive"]) pheno <- function() { ini({ tcl <- log(CL) tv <- log(Vd) eta.cl + eta.v ~ c(0.1, 0.01, 0.1) add.err <- add_err }) model({ cl <- exp(tcl + eta.cl) v <- exp(tv + eta.v) ke <- cl / v d/dt(A1) = -ke * A1 cp = A1 / v cp ~ add(add.err) }) } fit <- nlmixr2(pheno, pheno_sd, est = "saem", control = saemControl(print = 0)) ``` Print a summary of the model fit results: ```{r one-cmpt-result, eval=F} print(fit) # ── nlmixr² SAEM OBJF by FOCEi approximation ── # # Gaussian/Laplacian Likelihoods: AIC(fit) or fit$objf etc. # FOCEi CWRES & Likelihoods: addCwres(fit) # # ── Time (sec fit$time): ── # # setup covariance saem table compress other # elapsed 0.002899 0.025013 8.602 0.117 0.046 1.968088 # # ── Population Parameters (fit$parFixed or fit$parFixedDf): ── # # Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)% # tcl -5.02 0.0763 1.52 0.00662 (0.0057, 0.00768) 52.7 4.18% # tv 0.349 0.0548 15.7 1.42 (1.27, 1.58) 41.9 1.55% # add.err 2.74 2.74 # # Covariance Type (fit$covMethod): linFim # Correlations in between subject variability (BSV) matrix: # cor:eta.v,eta.cl # 0.927 # # # Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs) # Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink # Censoring (fit$censInformation): No censoring # # ── Fit Data (object fit is a modified tibble): ── # # A tibble: 155 × 17 # ID TIME DV PRED RES IPRED IRES IWRES eta.cl eta.v cp A1 cl v ke tad dosenum # # 1 1 2 17.3 17.5 -0.173 18.5 -1.16 -0.421 -0.0909 -0.0544 18.5 24.8 0.00604 1.34 0.00450 2 1 # 2 1 112. 31 28.1 2.94 30.0 1.02 0.371 -0.0909 -0.0544 30.0 40.3 0.00604 1.34 0.00450 4 10 # 3 2 2 9.7 10.5 -0.784 12.4 -2.69 -0.982 -0.238 -0.167 12.4 14.9 0.00521 1.20 0.00435 2 1 # # ℹ 152 more rows # # ℹ Use `print(n = ...)` to see more rows ``` If you plan to run multiple model structures, you can extract all parameters upfront at once and reuse them across models without repeating the lookup: ```{r two-cmt-iv-env, eval=F} CL <- as.numeric(est$Values[est$Parameters == "CL"]) Vd <- as.numeric(est$Values[est$Parameters == "Vd"]) Ka <- as.numeric(est$Values[est$Parameters == "Ka"]) Vmax <- as.numeric(est$Values[est$Parameters == "Vmax"]) Km <- as.numeric(est$Values[est$Parameters == "Km"]) Vc2cmt <- as.numeric(est$Values[est$Parameters == "Vc(2CMPT)"]) Vp2cmt <- as.numeric(est$Values[est$Parameters == "Vp(2CMPT)"]) Q2cmt <- as.numeric(est$Values[est$Parameters == "Q(2CMPT)"]) Vc3cmt <- as.numeric(est$Values[est$Parameters == "Vc(3CMPT)"]) Vp3cmt <- as.numeric(est$Values[est$Parameters == "Vp(3CMPT)"]) Vp23cmt <- as.numeric(est$Values[est$Parameters == "Vp2(3CMPT)"]) Q3cmt <- as.numeric(est$Values[est$Parameters == "Q(3CMPT)"]) Q23cmt <- as.numeric(est$Values[est$Parameters == "Q2(3CMPT)"]) add_err <- as.numeric(est$Values[est$Parameters == "Sigma additive"]) prop_err <- as.numeric(est$Values[est$Parameters == "Sigma proportional"]) two_cmt_iv <- function() { ini({ tcl <- CL tvc <- Vc2cmt tvp <- Vp2cmt tq <- Q2cmt eta.cl + eta.vc ~ c(1, 0.01, 1) add.err <- add_err }) model({ cl <- exp(tcl + eta.cl) vc <- exp(tvc + eta.vc) vp <- exp(tvp) q <- exp(tq) d/dt(central) = -(cl/vc) * central - (q/vc) * central + (q/vp) * peripheral d/dt(peripheral) = (q/vc) * central - (q/vp) * peripheral cp = central / vc cp ~ add(add.err) }) } fit_2cmt <- nlmixr2( two_cmt_iv, pheno_sd, est="saem",control=saemControl(print = 0,logLik = T)) ``` Print a summary of the model fit results: ```{r two-cmpt-result, eval=F} # print(fit_2cmt) # setup covariance saem table compress other # elapsed 0.003354 0.042016 85.371 0.141 0.042 2.58763 # # ── Population Parameters (fit_2cmt$parFixed or fit_2cmt$parFixedDf): ── # # Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)% # tcl -5.1 0.0879 1.72 0.00609 (0.00513, 0.00724) 57.8 -6.89% # tvc -0.554 0.156 28.1 0.575 (0.423, 0.78) 68.1 -9.03% # tvp -0.215 0.0859 39.9 0.806 (0.681, 0.954) # tq 0.103 0.633 617 1.11 (0.32, 3.83) # add.err 3.03 3.03 # # Covariance Type (fit_2cmt$covMethod): linFim # Correlations in between subject variability (BSV) matrix: # cor:eta.vc,eta.cl # 0.942 # # # Full BSV covariance (fit_2cmt$omega) or correlation (fit_2cmt$omegaR; diagonals=SDs) # Distribution stats (mean/skewness/kurtosis/p-value) available in fit_2cmt$shrink # Censoring (fit_2cmt$censInformation): No censoring # # ── Fit Data (object fit_2cmt is a modified tibble): ── # # A tibble: 155 × 19 # ID TIME DV PRED RES IPRED IRES IWRES eta.cl eta.vc cp central peripheral cl vc vp q tad dosenum # # 1 1 2 17.3 17.9 -0.611 18.5 -1.18 -0.390 -0.0706 -0.0769 18.5 9.83 14.9 0.00568 0.532 0.806 1.11 2 1 # 2 1 112. 31 29.2 1.76 30.5 0.457 0.151 -0.0706 -0.0769 30.5 16.3 24.7 0.00568 0.532 0.806 1.11 4 10 # 3 2 2 9.7 10.7 -1.05 12.1 -2.37 -0.782 -0.285 -0.304 12.1 5.12 9.75 0.00458 0.424 0.806 1.11 2 1 # # ℹ 152 more rows # # ℹ Use `print(n = ...)` to see more rows ``` ## See Also - [Integration with nlmixr2auto](workflow_nlmixr2auto.html) — automated model selection using `nlmixr2auto` - [Parameter Sweeping and Visualisation](parameter_sweeping.html) — inspect and visualise the parameter sweep results