mattfidler
commented
6 years ago It was very difficult for me to read this since you pasted your R session as output in the comment. I edited the comment so I could read it. I hope you don't mind. In the future please use ```R on the first line and then close it with ``` on the last line to enclose R code.
Also I'm unsure exactly what your questions were; I think these are your questions:
- Why do the liklihoods differ between nlme and saem?
- Why are the parameter estimates different between the different estimation methods?
- Why are the objects and outputs different between the two algoritnms?
Why are the liklihoods different between nlme and SAEM?
Each method makes it own likelhood approximation. The nlme approximation is not the same as the FOCEi approximation. The SAEM algorithm stotastically draws from the probabilitiy of the model and then uses for its maximizing the stoastic probability (See https://www.math.u-bordeaux.fr/MAS10/DOC/PDF-PRES/Lavielle.pdf). Therefore, the conventional wisdom is that liklihood between two different methods should not be compared because the approxmations are different and the optimizations between the algorithms are different. In fact for nlme when using REML you cannot really even compare two models at all (See https://stats.stackexchange.com/questions/116770/reml-or-ml-to-compare-two-mixed-effects-models-with-differing-fixed-effects-but)
Why are the parameter estimates different?
This is because of the different algorithms used to optimize the problem.
Why are the objects different?
This is manily due to a CRAN limitation on running vignettes in a timely manner and having software that is available on CRAN machines. I don't believe they have a recent python or sympy installation on those systems.
The changing of a traditional nlme object to a nlmixr-style option is controlled by the calc.resid option. If you decide to turn it on by calc.resid=TRUE (which is the default) or even drop it all together, you get much closer objects. Running the solved system gives:
> uif <- function() {
ini({
tka <- .5
tcl <- -3.2
tv <- -1
eta.ka ~ 1
eta.cl ~ 2
eta.v ~ 1
add.err <- 0.1
})
model({
ka <- exp(tka + eta.ka)
cl <- exp(tcl + eta.cl)
v <- exp(tv + eta.v)
linCmt() ~ add(add.err)
})
}
> fit <- nlmixr(uif, theo_sd, est="nlme")
**Iteration 1
LME step: Loglik: -182.2318, nlminb iterations: 1
reStruct parameters:
ID1 ID2 ID3
0.2783009 1.0022617 2.0758984
PNLS step: RSS = 63.14918
fixed effects: 0.4479976 -3.211924 -0.785932
iterations: 7
Convergence crit. (must all become <= tolerance = 1e-05):
fixed reStruct
0.2723747 3.2674463
... And even more output
> fit
nlmixr nlme fit by maximum likelihood (Solved)
FOCEi-based goodness of fit metrics:
OBJF AIC BIC Log-likelihood
296.5985 310.5985 330.7781 -148.2992
nlme-based goodness of fit metrics:
AIC BIC Log-likelihood
372.6488 392.8285 -179.3244
Time (sec; $time):
nlme setup FOCEi Evaulate covariance table
elapsed 2.17 21.44 0.5 0 0.96
Parameters ($par.fixed):
Estimate SE CV Untransformed (95%CI)
tka 0.445 0.198 44.4% 1.56 (1.06, 2.30)
tcl -3.21 0.0845 2.6% 0.0403 (0.0341, 0.0475)
tv -0.786 0.0463 5.9% 0.456 (0.416, 0.499)
add.err 0.692 0.692
Omega ($omega):
eta.ka eta.cl eta.v
eta.ka 0.412483 0.0000000 0.0000000
eta.cl 0.000000 0.0699386 0.0000000
eta.v 0.000000 0.0000000 0.0181094
Fit Data (object is a modified data.frame):
# A tibble: 132 x 15
ID TIME DV IPRED PRED IRES RES IWRES WRES CWRES CPRED
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0 0.740 0 0 0.740 0.740 1.07 0.0000207 -1.00 48128
2 1 0.250 2.84 3.86 2.82 -1.02 0.0218 -1.47 0.0316 -1.47 3.86
3 1 0.570 6.57 6.78 5.05 -0.207 1.52 -0.299 2.19 -0.299 6.78
# ... with 129 more rows, and 4 more variables: CRES <dbl>, eta.ka <dbl>,
# eta.cl <dbl>, eta.v <dbl>
> Note that:
- The 95% confidence intervals between each of the two model's parameters overlap, so they are not statistically different.
- When using the FOCEi likelihood approximation that NONMEM uses (which is the likelihood of the
SAEMand is output in the above object) gives anAICof310.5985which is close toSAEM's310.8482. - As noted previously the
nlmeliklihood is not the same as theFOCEilikelihood. In fact theAICvalues differ by >50 - To me there is no substiantial difference in either model.
Hopefully this answers your questions.
JSC67
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 (and by the way, is there a good way to make all of the following pasted text appear in a uniform typewriter font?): tutorial.txt tutorial output.txt
`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: