The tutorial (https://github.com/nlmixrdevelopment/nlmixr/blob/master/vignettes/running_nlmixr.Rmd) has code for compiling a model for closed form solution using linCmt() and for ODE solution, and also shows how one might use NMLE and SAEM algorithms. With the single dose theophylline data at hand, this gives 4 possibilities. I followed along the tutorial but made a few changes for clarity and for comparison and found that the NMLE result for closed form and for ODE agree very well (but not exactly), whereas the SAEM objective function values, AIC, BIC, and log-likelihood are very different for the explicit versus the ODE models and do not agree with the NMLE results either. The SAEM results seem to be way off. Below is the code and results: R version 3.4.3 (2017-11-30) -- "Kite-Eating Tree" Copyright (C) 2017 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > setwd("C:/Users/James/Documents/pharmacometrics/practice/nlmixr") > #https://github.com/nlmixrdevelopment/nlmixr/blob/master/vignettes/running_nlmixr.Rmd > > # Preliminaries > library(ggplot2) > library(nlmixr) Loading required package: nlme Loading required package: RxODE > str(theo_sd) 'data.frame': 144 obs. of 7 variables: $ ID : int 1 1 1 1 1 1 1 1 1 1 ... $ TIME: num 0 0 0.25 0.57 1.12 2.02 3.82 5.1 7.03 9.05 ... $ DV : num 0 0.74 2.84 6.57 10.5 9.66 8.58 8.36 7.47 6.89 ... $ AMT : num 4.02 0 0 0 0 0 0 0 0 0 ... $ EVID: int 101 0 0 0 0 0 0 0 0 0 ... $ CMT : int 1 2 2 2 2 2 2 2 2 2 ... $ WT : num 79.6 79.6 79.6 79.6 79.6 79.6 79.6 79.6 79.6 79.6 ... > ggplot(theo_sd, aes(TIME, DV)) + + geom_line(aes(group=ID), col="red") + + scale_x_continuous("Time (h)") + + scale_y_continuous("Concentration") + + labs(title="Theophylline single-dose", + subtitle="Concentration vs. time by individual") > # Try fitting a simple one-compartment PK model to this small dataset. > # 'uif' for Unified Interface Function, lin for linCmt() > uif.lin <- function() { + ini({ + logtcl <- -3.2 # log typical value Cl (L/hr) + logtka <- 0.5 # log typical value Ka (1/hr) + logtv <- -1 # log typical value V (L) + + # error model + add.err <- 0.1 + + # Initial estimates for IIV variances + # Labels work for single parameters. + eta.cl ~ 2 + eta.ka ~ 1 ## BSV Ka + eta.v ~ 1 + }) + model({ + cl <- exp(logtcl + eta.cl) + ka <- exp(logtka + eta.ka) + v <- exp(logtv + eta.v) + linCmt() ~ add(add.err) + }) + } > # We can alternatively express the same model by ODEs: > uif.ode <- function() { + ini({ + logtcl <- -3.2 # log typical value Cl (L/hr) + logtka <- 0.5 # log typical value Ka (1/hr) + logtv <- -1 # log typical value V (L) + + # error model + add.err <- 0.1 + + # Initial estimates for IIV variances + # Labels work for single parameters. + eta.cl ~ 2 + eta.ka ~ 1 ## BSV Ka + eta.v ~ 1 + }) + model({ + cl <- exp(logtcl + eta.cl) + ka <- exp(logtka + eta.ka) + v <- exp(logtv + eta.v) + d/dt(depot) = -ka * depot + d/dt(center) = ka * depot - cl / v * center + cp = center / v + cp ~ add(add.err) + }) + } > # NLME fittings > > fit.nlme.lin <- nlmixr(uif.lin, theo_sd, est="nlme", calc.resid=FALSE) **Iteration 1 LME step: Loglik: -182.2318, nlminb iterations: 1 reStruct parameters: ID1 ID2 ID3 1.0022617 0.2783009 2.0758984 PNLS step: RSS = 63.14922 fixed effects: -3.211937 0.4479979 -0.7859318 iterations: 7 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.272375 3.267309 **Iteration 2 LME step: Loglik: -179.291, nlminb iterations: 9 reStruct parameters: ID1 ID2 ID3 0.96385527 0.08333435 1.63552295 PNLS step: RSS = 63.28221 fixed effects: -3.211535 0.4441505 -0.7863391 iterations: 7 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.008662381 0.172135330 **Iteration 3 LME step: Loglik: -179.3381, nlminb iterations: 8 reStruct parameters: ID1 ID2 ID3 0.96129928 0.07085834 1.63885629 PNLS step: RSS = 63.2196 fixed effects: -3.211732 0.4458486 -0.7861743 iterations: 7 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.003808629 0.082919683 **Iteration 4 LME step: Loglik: -179.3201, nlminb iterations: 7 reStruct parameters: ID1 ID2 ID3 0.96243507 0.07618717 1.63736945 PNLS step: RSS = 63.22991 fixed effects: -3.211823 0.4451494 -0.7862432 iterations: 7 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.001570583 0.013884126 **Iteration 5 LME step: Loglik: -179.3277, nlminb iterations: 4 reStruct parameters: ID1 ID2 ID3 0.96197111 0.07396744 1.63796900 PNLS step: RSS = 63.24112 fixed effects: -3.211722 0.4454457 -0.7862145 iterations: 6 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.0006651255 0.0147544693 **Iteration 6 LME step: Loglik: -179.3246, nlminb iterations: 1 reStruct parameters: ID1 ID2 ID3 0.96215433 0.07487426 1.63772430 PNLS step: RSS = 63.24112 fixed effects: -3.211722 0.4454457 -0.7862145 iterations: 1 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.000000e+00 3.910435e-10 Warning message: In nlmixrUI.nlme.var(obj) : Initial condition for additive error ignored with nlme > fit.nlme.ode <- nlmixr(uif.ode, theo_sd, est="nlme", calc.resid=FALSE) **Iteration 1 LME step: Loglik: -182.2318, nlminb iterations: 1 reStruct parameters: ID1 ID2 ID3 1.0022618 0.2783008 2.0759012 PNLS step: RSS = 63.24481 fixed effects: -3.211675 0.4451503 -0.7862385 iterations: 7 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.2718787 2.7010441 **Iteration 2 LME step: Loglik: -179.3267, nlminb iterations: 6 reStruct parameters: ID1 ID2 ID3 0.96197025 0.07396818 1.63806766 PNLS step: RSS = 63.23773 fixed effects: -3.211876 0.4453302 -0.7862234 iterations: 2 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.0004040574 0.0088556330 **Iteration 3 LME step: Loglik: -179.3253, nlminb iterations: 1 reStruct parameters: ID1 ID2 ID3 0.96208503 0.07453515 1.63782577 PNLS step: RSS = 63.23773 fixed effects: -3.211876 0.4453302 -0.7862234 iterations: 1 Convergence crit. (must all become <= tolerance = 1e-05): fixed reStruct 0.000000e+00 3.183879e-10 Warning message: In nlmixrUI.nlme.var(obj) : Initial condition for additive error ignored with nlme > # SAEM fittings > > # For SAEM, it's nice to avoid printing each C++ iteration to the monitor. > sink("saem output") > fit.saem.lin <- nlmixr(uif.lin, theo_sd, est="saem") # works with spaces in path Compiling SAEM user function...C:/RTools/3.4/mingw_64/bin/g++ -I"C:/R/R-34~1.3/include" -DNDEBUG -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/nlmixr/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/STANHE~1/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/Rcpp/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/RCPPAR~1/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/RCPPEI~1/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/BH/include -O2 -Wall -mtune=generic -c saem137829e957ffx64.cpp -o saem137829e957ffx64.o C:/RTools/3.4/mingw_64/bin/g++ -shared -s -static-libgcc -o saem137829e957ffx64.dll tmp.def saem137829e957ffx64.o -LC:/R/R-34~1.3/bin/x64 -lRblas -LC:/R/R-34~1.3/bin/x64 -lRlapack -lgfortran -lm -lquadmath -LC:/R/R-34~1.3/bin/x64 -lR done. Using sympy via SnakeCharmR C:/RTools/3.4/mingw_64/bin/gcc -I"C:/R/R-34~1.3/include" -DNDEBUG -O2 -Wall -std=gnu99 -mtune=generic -c rx_920cdb7f8c8c7b5e42b1a8db66863398_x64.c -o rx_920cdb7f8c8c7b5e42b1a8db66863398_x64.o C:/RTools/3.4/mingw_64/bin/gcc -shared -s -static-libgcc -o rx_920cdb7f8c8c7b5e42b1a8db66863398_x64.dll tmp.def rx_920cdb7f8c8c7b5e42b1a8db66863398_x64.o -LC:/R/R-34~1.3/bin/x64 -lRblas -LC:/R/R-34~1.3/bin/x64 -lRlapack -lgfortran -lm -lquadmath -LC:/R/R-34~1.3/bin/x64 -lR The model has been setup to calculate residuals. It will be cached for future runs. Calculating Table Variables... done nlmixr UI combined dataset and properties $ par.hist : Parameter history (if available) $ par.hist.stacked : Parameter history in stacked form for easy plotting (if available) $ par.fixed : Fixed Effect Parameter Table $ eta : Individual Parameter Estimates $ seed : Seed (if applicable) $ model.name : Model name (from R function) $ data.name : Name of R data input > # Using RStudio, in Environment tab, what is object curi with value 0 has appeared. > # What is curi? > # NB: The following won't work with spaces in file full file path, don't know why: > fit.saem.ode <- nlmixr(uif.ode, theo_sd, est="saem") # fails with spaces in path Compiling RxODE differential equations...done. C:/RTools/3.4/mingw_64/bin/g++ -I"C:/R/R-34~1.3/include" -DNDEBUG -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/nlmixr/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/STANHE~1/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/Rcpp/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/RCPPAR~1/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/RCPPEI~1/include -IC:/Users/James/DOCUME~1/R/WIN-LI~1/3.4/BH/include -O2 -Wall -mtune=generic -c saem13783d7445b6x64.cpp -o saem13783d7445b6x64.o C:/RTools/3.4/mingw_64/bin/g++ -shared -s -static-libgcc -o saem13783d7445b6x64.dll tmp.def saem13783d7445b6x64.o C:/Users/James/Documents/pharmacometrics/practice/nlmixr/rx_7cbddd1a680b22377f3ee375ef2ddda8_x64.dll -LC:/R/R-34~1.3/bin/x64 -lRblas -LC:/R/R-34~1.3/bin/x64 -lRlapack -lgfortran -lm -lquadmath -LC:/R/R-34~1.3/bin/x64 -lR done. ## Calculate ETA-based prediction and error derivatives: Calculate Jacobian...................done. Calculate sensitivities....... done. ## Calculate d(f)/d(eta) ## ... ## done ## ... ## done C:/RTools/3.4/mingw_64/bin/gcc -I"C:/R/R-34~1.3/include" -DNDEBUG -O2 -Wall -std=gnu99 -mtune=generic -c rx_32a3413ee26a7ebdb0539062d9fc430e_x64.c -o rx_32a3413ee26a7ebdb0539062d9fc430e_x64.o C:/RTools/3.4/mingw_64/bin/gcc -shared -s -static-libgcc -o rx_32a3413ee26a7ebdb0539062d9fc430e_x64.dll tmp.def rx_32a3413ee26a7ebdb0539062d9fc430e_x64.o -LC:/R/R-34~1.3/bin/x64 -lRblas -LC:/R/R-34~1.3/bin/x64 -lRlapack -lgfortran -lm -lquadmath -LC:/R/R-34~1.3/bin/x64 -lR C:/RTools/3.4/mingw_64/bin/gcc -I"C:/R/R-34~1.3/include" -DNDEBUG -O2 -Wall -std=gnu99 -mtune=generic -c rx_a2baf9a7f76ed6696b7c9c5507ef49d7_x64.c -o rx_a2baf9a7f76ed6696b7c9c5507ef49d7_x64.o C:/RTools/3.4/mingw_64/bin/gcc -shared -s -static-libgcc -o rx_a2baf9a7f76ed6696b7c9c5507ef49d7_x64.dll tmp.def rx_a2baf9a7f76ed6696b7c9c5507ef49d7_x64.o -LC:/R/R-34~1.3/bin/x64 -lRblas -LC:/R/R-34~1.3/bin/x64 -lRlapack -lgfortran -lm -lquadmath -LC:/R/R-34~1.3/bin/x64 -lR The model-based sensitivities have been calculated. It will be cached for future runs. Calculating Table Variables... done nlmixr UI combined dataset and properties $ par.hist : Parameter history (if available) $ par.hist.stacked : Parameter history in stacked form for easy plotting (if available) $ par.fixed : Fixed Effect Parameter Table $ eta : Individual Parameter Estimates $ seed : Seed (if applicable) $ model.name : Model name (from R function) $ data.name : Name of R data input > sink() > # Now compare fits for explicit (i.e., linCmt() ) model vs. ODE, and by algorithms: > fit.nlme.lin # Log-likelihood: -179.3246 Nonlinear mixed-effects model fit by maximum likelihood Model: DV ~ (nlmeModList("user_fn"))(logtcl, eta.cl, logtka, eta.ka, logtv, eta.v, TIME, ID) Log-likelihood: -179.3246 Fixed: logtcl + logtka + logtv ~ 1 logtcl logtka logtv -3.2117217 0.4454457 -0.7862145 Random effects: Formula: list(eta.cl ~ 1, eta.ka ~ 1, eta.v ~ 1) Level: ID Structure: Diagonal eta.cl eta.ka eta.v Residual StdDev: 0.2644566 0.6422369 0.134573 0.6921699 Number of Observations: 132 Number of Groups: 12 > fit.nlme.ode # Log-likelihood: -179.3253 Nonlinear mixed-effects model fit by maximum likelihood Model: DV ~ (nlmeModList("user_fn"))(logtcl, eta.cl, logtka, eta.ka, logtv, eta.v, TIME, ID) Log-likelihood: -179.3253 Fixed: logtcl + logtka + logtv ~ 1 logtcl logtka logtv -3.2118758 0.4453302 -0.7862234 Random effects: Formula: list(eta.cl ~ 1, eta.ka ~ 1, eta.v ~ 1) Level: ID Structure: Diagonal eta.cl eta.ka eta.v Residual StdDev: 0.2644679 0.6424376 0.1345558 0.6921515 Number of Observations: 132 Number of Groups: 12 > fit.saem.lin # Log-likelihood: -148.4444 nlmixr SAEM fit (Solved) OBJF AIC BIC Log-likelihood 296.8482 310.8482 331.0278 -148.4241 Time (sec; $time): saem setup Likelihood Calculation covariance table elapsed 53.9 65.11 2.09 0 1.5 Parameters ($par.fixed): Parameter Estimate SE CV Untransformed logtcl log typical value Cl (L/hr) -3.22 0.0816 2.5% 0.0401 logtka log typical value Ka (1/hr) 0.450 0.194 43.0% 1.57 logtv log typical value V (L) -0.784 0.0435 5.6% 0.457 add.err add.err 0.692 0.692 (95%CI) logtcl (0.0342, 0.0471) logtka (1.07, 2.29) logtv (0.419, 0.497) add.err Omega ($omega): eta.cl eta.ka eta.v eta.cl 0.0696184 0.0000000 0.00000000 eta.ka 0.0000000 0.4187747 0.00000000 eta.v 0.0000000 0.0000000 0.01819676 Fit Data (object is a modified data.frame): # A tibble: 132 x 15 ID TIME DV IPRED PRED IRES RES IWRES WRES CWRES CPRED CRES 1 1 0 0.740 0 0 0.740 0.740 1.07 2.06e-5 -0.999 4.81e4 -4.81e+4 2 1 0.250 2.84 3.85 2.82 -1.01 0.0178 -1.46 2.57e-2 -1.46 3.85e0 -1.01e+0 3 1 0.570 6.57 6.76 5.06 -0.191 1.51 -0.276 2.19e+0 -0.276 6.76e0 -1.91e-1 # ... with 129 more rows, and 3 more variables: eta.cl , eta.ka , # eta.v > fit.saem.ode # Log-likelihood: -58.07448 nlmixr SAEM fit (ODE) OBJF AIC BIC Log-likelihood 116.149 130.149 150.3286 -58.07448 Time (sec; $time): saem setup Likelihood Calculation covariance table elapsed 49.99 17.71 0.28 0 1.18 Parameters ($par.fixed): Parameter Estimate SE CV Untransformed logtcl log typical value Cl (L/hr) -3.22 0.0818 2.5% 0.0400 logtka log typical value Ka (1/hr) 0.458 0.192 41.9% 1.58 logtv log typical value V (L) -0.782 0.0437 5.6% 0.458 add.err add.err 0.694 0.694 (95%CI) logtcl (0.0341, 0.0470) logtka (1.09, 2.30) logtv (0.420, 0.499) add.err Omega ($omega): eta.cl eta.ka eta.v eta.cl 0.07190448 0.0000000 0.00000000 eta.ka 0.00000000 0.4174639 0.00000000 eta.v 0.00000000 0.0000000 0.01818897 Fit Data (object is a modified data.frame): # A tibble: 132 x 15 ID TIME DV IPRED PRED IRES RES IWRES WRES CWRES CPRED CRES 1 1 0 0.740 0 0 0.740 0.740 1.07 1.07 1.07 0 0.740 2 1 0.250 2.84 3.90 2.83 -1.06 0.00644 -1.53 0.00382 0.0813 2.66 0.177 3 1 0.570 6.57 6.83 5.07 -0.255 1.50 -0.368 0.675 0.632 4.83 1.74 # ... with 129 more rows, and 3 more variables: eta.cl , eta.ka , # eta.v > # Why are they so very different? > # Why does object x appear (1 obs. of 4 variables) appear after fit.saem.ode only? > > # Suggestion: Use same printing format for models with either algorithm.