diff --git a/DESCRIPTION b/DESCRIPTION
index 41f7cb12b..ff1200d60 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -45,6 +45,6 @@ LinkingTo: dparser(>= 0.1.8), Rcpp (>= 0.12.3),
 URL: https://github.com/nlmixrdevelopment/nlmixr
 LazyData: true
 VignetteBuilder: knitr
-RoxygenNote: 6.0.1
+RoxygenNote: 6.1.0
 Packaged: 2018-04-30 01:43:13 UTC; annie
 Maintainer: Wenping Wang <wwang8198@gmail.com>
diff --git a/R/foceiFit.R b/R/foceiFit.R
index edd2f5b3e..37886482a 100644
--- a/R/foceiFit.R
+++ b/R/foceiFit.R
@@ -168,6 +168,9 @@ is.latex <- function() {
 ##' @param noAbort Boolean to indicate if you should abort the FOCEi
 ##'     evaluation if it runs into troubles.  (default TRUE)
 ##'
+##' @param interaction Boolean indicate FOCEi should be used (TRUE)
+##'     instead of FOCE (FALSE)
+##'
 ##' @inheritParams RxODE::rxSolve
 ##'
 ##' @details
@@ -228,6 +231,7 @@ foceiControl <- function(sigdig=4,
                          boundTol=NULL,
                          calcTables=TRUE,
                          noAbort=TRUE,
+                         interaction=TRUE,
                          ..., stiff){
     if (is.null(boundTol)){
         boundTol <- 5 * 10 ^ (-sigdig + 1)
@@ -336,7 +340,8 @@ foceiControl <- function(sigdig=4,
                  boundTol=as.double(boundTol),
                  calcTables=calcTables,
                  printNcol=as.integer(printNcol),
-                 noAbort=as.integer(noAbort));
+                 noAbort=as.integer(noAbort),
+                 interaction=as.integer(interaction));
     class(.ret) <- "foceiControl"
     return(.ret);
 }
@@ -448,15 +453,20 @@ constructLinCmt <- function(fun){
 ##' w7$EVID <- 0
 ##' w7$AMT <- 0
 ##'
-##' ## Wang2007 prop error 39.458 for NONMEM, nlmixr matches.
-##' fitP <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
+##' ## Wang2007 prop error 39.458 for NONMEM FOCEi, nlmixr matches.
+##' fitPi <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
 ##'      control=foceiControl(maxOuterIterations=0,covMethod=""))
 ##'
+##' ## Wang2007 prop error 39.207 for NONMEM FOCE; nlmixr matches.
+##' fitP <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
+##'      control=foceiControl(maxOuterIterations=0,covMethod="",
+##'      interaction=FALSE))
+##'
 ##' ## Note if you have the etas you can evaluate the likelihood
 ##' ## of an arbitrary model.  It doesn't have to be solved by
 ##' ## FOCEi
 ##'
-##' etaMat <- matrix(fitP$eta[,-1])
+##' etaMat <- matrix(fitPi$eta[,-1])
 ##'
 ##' fitP2 <- foceiFit(w7, inits, mypar2,mod,pred,errProp, etaMat=etaMat,
 ##'       control=foceiControl(maxOuterIterations=0,maxInnerIterations=0,
@@ -786,7 +796,9 @@ foceiFit.data.frame0 <- function(data,
     if (!exists("noLik", envir=.ret)){
         .ret$model <- RxODE::rxSymPySetupPred(model, pred, PKpars, err, grad=(control$derivMethod == 2L),
                                               pred.minus.dv=TRUE, sum.prod=control$sumProd,
-                                              theta.derivs=FALSE, optExpression=control$optExpression, run.internal=TRUE);
+                                              theta.derivs=FALSE, optExpression=control$optExpression,
+                                              interaction=(control$interaction == 1L),
+                                              run.internal=TRUE);
 
         .covNames <- .parNames <- RxODE::rxParams(.ret$model$pred.only);
         .covNames <- .covNames[regexpr(rex::rex(start, or("THETA", "ETA"), "[", numbers, "]", end), .covNames) == -1];
@@ -1264,8 +1276,8 @@ print.nlmixrFitCore <- function(x, ...){
     }
     .posthoc <- (x$control$maxOuterIterations == 0L & x$control$maxInnerIterations > 0L)
     .posthoc <- ifelse(.posthoc, paste0(crayon::bold(" posthoc"), " estimation"), " fit");
-    message(cli::rule(paste0(crayon::bold$blue("nlmix"), crayon::bold$red("r"), " ", crayon::bold$yellow(x$method),.posthoc,
-                             x$extra)))
+    message(cli::rule(paste0(crayon::bold$blue("nlmix"), crayon::bold$red("r"), " ", crayon::bold$yellow(x$method),
+                             x$extra, .posthoc)))
     print(x$objDf)
     message(paste0("\n", cli::rule(paste0(crayon::bold("Time"), " (sec; ", crayon::yellow(.bound), crayon::bold$blue("$time"), "):"))));
     print(x$time)
diff --git a/R/simulate.R b/R/simulate.R
index da11fe635..e6ee09165 100644
--- a/R/simulate.R
+++ b/R/simulate.R
@@ -90,7 +90,8 @@
 ##' simulated from the uncertainty in the Omega/Sigma matrices based
 ##' on the number of subjects and observations the model was based on.
 ##'
-##' @inheritParams RxODE::rxSolve
+##' @param object nlmixr object
+##' @param ... Other arguments sent to \code{rxSolve}
 ##'
 ##' @export
 nlmixrSim <- function(object, ...){
@@ -208,26 +209,13 @@ plot.nlmixrSim <- function(x, y, ...){
 ##' @param ipred Flag to calculate individual predictions. When
 ##'     \code{ipred} is \code{TRUE}, calculate individual predictions.
 ##'     When \code{ipred} is \code{FALSE}, set calculate typical population predations.
-##'     When \code{ipred} is \code{NA}, calculateboth individual and
+##'     When \code{ipred} is \code{NA}, calculate both individual and
 ##'     population predictions.
 ##'
 ##' @inheritParams RxODE::rxSolve
 ##'
 ##' @export
-nlmixrPred <- function(object, params=NULL, events=NULL, inits = NULL, scale = NULL,
-                       covs = NULL, method = c("liblsoda", "lsoda", "dop853"),
-                       transitAbs = NULL, atol = 1.0e-6, rtol = 1.0e-4,
-                       maxsteps = 5000L, hmin = 0L, hmax = NULL, hini = 0L, maxordn = 12L, maxords = 5L, ...,
-                       cores, covsInterpolation = c("linear", "locf", "nocb", "midpoint"),
-                       addCov = FALSE, matrix = FALSE, sigma = NULL, sigmaDf = NULL,
-                       nCoresRV = 1L, sigmaIsChol = FALSE, nDisplayProgress=10000L,
-                       amountUnits = NA_character_, timeUnits = "hours", stiff,
-                       theta = NULL, eta = NULL, addDosing=FALSE, updateObject=FALSE,doSolve=TRUE,
-                       omega = NULL, omegaDf = NULL, omegaIsChol = FALSE,
-                       nSub = 1L, thetaMat = NULL, thetaDf = NULL, thetaIsChol = FALSE,
-                       nStud = 1L, dfSub=0.0, dfObs=0.0, returnType=c("data.frame", "rxSolve", "matrix"),
-                       seed=NULL, nsim=NULL,
-                       ipred=FALSE){
+nlmixrPred <- function(object, ..., ipred=FALSE){
     lst <- as.list(match.call()[-1]);
     if (RxODE::rxIs(lst$params, "rx.event")){
         if (!is.null(lst$events)){
diff --git a/man/configsaem.Rd b/man/configsaem.Rd
index b74b6c6f0..0396929ab 100644
--- a/man/configsaem.Rd
+++ b/man/configsaem.Rd
@@ -4,10 +4,10 @@
 \alias{configsaem}
 \title{Configure an SAEM model}
 \usage{
-configsaem(model, data, inits, mcmc = list(niter = c(200, 300), nmc = 3, nu =
-  c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1, transitAbs
-  = 0), distribution = c("normal", "poisson", "binomial"), seed = 99,
-  fixed = NULL, DEBUG = 0)
+configsaem(model, data, inits, mcmc = list(niter = c(200, 300), nmc = 3,
+  nu = c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1,
+  transitAbs = 0), distribution = c("normal", "poisson", "binomial"),
+  seed = 99, fixed = NULL, DEBUG = 0)
 }
 \arguments{
 \item{model}{a compiled saem model by gen_saem_user_fn()}
diff --git a/man/dot-collectWarnings.Rd b/man/dot-collectWarnings.Rd
new file mode 100644
index 000000000..c48cc564f
--- /dev/null
+++ b/man/dot-collectWarnings.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.R
+\name{.collectWarnings}
+\alias{.collectWarnings}
+\title{Collect warnings and just warn once.}
+\usage{
+.collectWarnings(expr)
+}
+\arguments{
+\item{expr}{R expression}
+}
+\value{
+The value of the expression
+}
+\description{
+Collect warnings and just warn once.
+}
+\author{
+Matthew L. Fidler
+}
diff --git a/man/dynmodel.Rd b/man/dynmodel.Rd
index 3e9089c80..76587371b 100644
--- a/man/dynmodel.Rd
+++ b/man/dynmodel.Rd
@@ -5,8 +5,8 @@
 \title{Fit a non-population dynamic model}
 \usage{
 dynmodel(system, model, evTable, inits, data, fixPars = NULL,
-  method = c("Nelder-Mead", "L-BFGS-B", "PORT"), control = list(ftol_rel =
-  1e-06, maxeval = 999), squared = T)
+  method = c("Nelder-Mead", "L-BFGS-B", "PORT"),
+  control = list(ftol_rel = 1e-06, maxeval = 999), squared = T)
 }
 \arguments{
 \item{system}{an RxODE object}
diff --git a/man/foceiControl.Rd b/man/foceiControl.Rd
index c404c7c3c..0aef15e85 100644
--- a/man/foceiControl.Rd
+++ b/man/foceiControl.Rd
@@ -10,14 +10,14 @@ foceiControl(sigdig = 4, epsilon = NULL, maxInnerIterations = 1000,
   maxstepsOde = 5000L, hmin = 0L, hmax = NULL, hini = 0,
   maxordn = 12L, maxords = 5L, cores, covsInterpolation = c("linear",
   "locf", "nocb", "midpoint"), printInner = 0L, print = 1L,
-  printNcol = 8L, scaleTo = 1, scaleObjective = 1, derivEps = c(1e-05,
-  1e-05), derivMethod = c("forward", "central"),
-  covDerivMethod = c("forward", "central"), covMethod = c("r,s", "r", "s",
-  ""), lbfgsLmm = 50L, lbfgsPgtol = 0, lbfgsFactr = NULL, eigen = TRUE,
-  addPosthoc = TRUE, diagXform = c("sqrt", "log", "identity"),
-  sumProd = FALSE, optExpression = TRUE, ci = 0.95,
+  printNcol = 8L, scaleTo = 1, scaleObjective = 1,
+  derivEps = c(1e-05, 1e-05), derivMethod = c("forward", "central"),
+  covDerivMethod = c("forward", "central"), covMethod = c("r,s", "r",
+  "s", ""), lbfgsLmm = 50L, lbfgsPgtol = 0, lbfgsFactr = NULL,
+  eigen = TRUE, addPosthoc = TRUE, diagXform = c("sqrt", "log",
+  "identity"), sumProd = FALSE, optExpression = TRUE, ci = 0.95,
   useColor = crayon::has_color(), boundTol = NULL, calcTables = TRUE,
-  noAbort = TRUE, ..., stiff)
+  noAbort = TRUE, interaction = TRUE, ..., stiff)
 }
 \arguments{
 \item{sigdig}{Optimization significant digits. This controls:
@@ -229,6 +229,9 @@ will calculate tables. By default this is \code{TRUE}}
 \item{noAbort}{Boolean to indicate if you should abort the FOCEi
 evaluation if it runs into troubles.  (default TRUE)}
 
+\item{interaction}{Boolean indicate FOCEi should be used (TRUE)
+instead of FOCE (FALSE)}
+
 \item{...}{Ignored parameters}
 
 \item{stiff}{a logical (\code{TRUE} by default) indicating whether
diff --git a/man/foceiFit.Rd b/man/foceiFit.Rd
index bc57c7ea3..8642fcd30 100644
--- a/man/foceiFit.Rd
+++ b/man/foceiFit.Rd
@@ -93,15 +93,20 @@ w7$DV <- w7$Y
 w7$EVID <- 0
 w7$AMT <- 0
 
-## Wang2007 prop error 39.458 for NONMEM, nlmixr matches.
-fitP <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
+## Wang2007 prop error 39.458 for NONMEM FOCEi, nlmixr matches.
+fitPi <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
      control=foceiControl(maxOuterIterations=0,covMethod=""))
 
+## Wang2007 prop error 39.207 for NONMEM FOCE; nlmixr matches.
+fitP <- foceiFit(w7, inits, mypar2,mod,pred,errProp,
+     control=foceiControl(maxOuterIterations=0,covMethod="",
+     interaction=FALSE))
+
 ## Note if you have the etas you can evaluate the likelihood
 ## of an arbitrary model.  It doesn't have to be solved by
 ## FOCEi
 
-etaMat <- matrix(fitP$eta[,-1])
+etaMat <- matrix(fitPi$eta[,-1])
 
 fitP2 <- foceiFit(w7, inits, mypar2,mod,pred,errProp, etaMat=etaMat,
       control=foceiControl(maxOuterIterations=0,maxInnerIterations=0,
diff --git a/man/gen_saem_user_fn.Rd b/man/gen_saem_user_fn.Rd
index df62ed743..7d56edb8d 100644
--- a/man/gen_saem_user_fn.Rd
+++ b/man/gen_saem_user_fn.Rd
@@ -4,8 +4,8 @@
 \alias{gen_saem_user_fn}
 \title{Generate an SAEM model}
 \usage{
-gen_saem_user_fn(model, PKpars = attr(model, "default.pars"), pred = NULL,
-  inPars = NULL)
+gen_saem_user_fn(model, PKpars = attr(model, "default.pars"),
+  pred = NULL, inPars = NULL)
 }
 \arguments{
 \item{model}{a compiled SAEM model by gen_saem_user_fn()}
diff --git a/man/lincmt.Rd b/man/lincmt.Rd
index 53988c76e..bc59c138f 100644
--- a/man/lincmt.Rd
+++ b/man/lincmt.Rd
@@ -4,7 +4,8 @@
 \alias{lincmt}
 \title{Parameters for a linear compartment model for SAEM}
 \usage{
-lincmt(ncmt, oral = T, tlag = F, infusion = F, parameterization = 1)
+lincmt(ncmt, oral = T, tlag = F, infusion = F,
+  parameterization = 1)
 }
 \arguments{
 \item{ncmt}{number of compartments}
diff --git a/man/nlme_lin_cmpt.Rd b/man/nlme_lin_cmpt.Rd
index 81c9c2def..822b3985f 100644
--- a/man/nlme_lin_cmpt.Rd
+++ b/man/nlme_lin_cmpt.Rd
@@ -7,16 +7,16 @@
 \title{Fit nlme-based linear compartment mixed-effect model using closed form solution}
 \usage{
 nlme_lin_cmpt(dat, parModel, ncmt, oral = TRUE, infusion = FALSE,
-  tlag = FALSE, parameterization = 1, parTrans = .getParfn(oral, ncmt,
-  parameterization, tlag), mcCores = 1, ...)
+  tlag = FALSE, parameterization = 1, parTrans = .getParfn(oral,
+  ncmt, parameterization, tlag), mcCores = 1, ...)
 
 nlmeLinCmpt(dat, parModel, ncmt, oral = TRUE, infusion = FALSE,
-  tlag = FALSE, parameterization = 1, parTrans = .getParfn(oral, ncmt,
-  parameterization, tlag), mcCores = 1, ...)
+  tlag = FALSE, parameterization = 1, parTrans = .getParfn(oral,
+  ncmt, parameterization, tlag), mcCores = 1, ...)
 
 nlmeLinCmt(dat, parModel, ncmt, oral = TRUE, infusion = FALSE,
-  tlag = FALSE, parameterization = 1, parTrans = .getParfn(oral, ncmt,
-  parameterization, tlag), mcCores = 1, ...)
+  tlag = FALSE, parameterization = 1, parTrans = .getParfn(oral,
+  ncmt, parameterization, tlag), mcCores = 1, ...)
 }
 \arguments{
 \item{dat}{data to be fitted}
diff --git a/man/nlme_ode.Rd b/man/nlme_ode.Rd
index 725ae0875..e71e93edb 100644
--- a/man/nlme_ode.Rd
+++ b/man/nlme_ode.Rd
@@ -5,15 +5,15 @@
 \alias{nlmeOde}
 \title{Fit nlme-based mixed-effect model using ODE implementation}
 \usage{
-nlme_ode(dat.o, model, parModel, parTrans, response, responseScaler = NULL,
-  transitAbs = FALSE, atol = 1e-06, rtol = 1e-04, maxsteps = 5000,
-  hmin = 0, hmax = NULL, hini = 0, maxordn = 12, maxords = 5,
-  debugODE = FALSE, mcCores = 1, ...)
-
-nlmeOde(dat.o, model, parModel, parTrans, response, responseScaler = NULL,
-  transitAbs = FALSE, atol = 1e-06, rtol = 1e-04, maxsteps = 5000,
-  hmin = 0, hmax = NULL, hini = 0, maxordn = 12, maxords = 5,
-  debugODE = FALSE, mcCores = 1, ...)
+nlme_ode(dat.o, model, parModel, parTrans, response,
+  responseScaler = NULL, transitAbs = FALSE, atol = 1e-06,
+  rtol = 1e-04, maxsteps = 5000, hmin = 0, hmax = NULL, hini = 0,
+  maxordn = 12, maxords = 5, debugODE = FALSE, mcCores = 1, ...)
+
+nlmeOde(dat.o, model, parModel, parTrans, response,
+  responseScaler = NULL, transitAbs = FALSE, atol = 1e-06,
+  rtol = 1e-04, maxsteps = 5000, hmin = 0, hmax = NULL, hini = 0,
+  maxordn = 12, maxords = 5, debugODE = FALSE, mcCores = 1, ...)
 }
 \arguments{
 \item{dat.o}{data to be fitted}
diff --git a/man/nlmixr.Rd b/man/nlmixr.Rd
index cd36dddcf..51b29a0af 100644
--- a/man/nlmixr.Rd
+++ b/man/nlmixr.Rd
@@ -6,13 +6,14 @@
 \alias{nlmixr.nlmixrUI}
 \title{nlmixr fits population PK and PKPD non-linear mixed effects models.}
 \usage{
-nlmixr(object, data, est = NULL, control = list(), table = tableControl(),
-  ...)
+nlmixr(object, data, est = NULL, control = list(),
+  table = tableControl(), ...)
 
 \method{nlmixr}{function}(object, data, est = NULL, control = list(),
   table = tableControl(), ...)
 
-\method{nlmixr}{nlmixrUI}(object, data, est = NULL, control = list(), ...)
+\method{nlmixr}{nlmixrUI}(object, data, est = NULL, control = list(),
+  ...)
 }
 \arguments{
 \item{object}{Fitted object or function specifying the model.}
diff --git a/man/nlmixrAugPred.Rd b/man/nlmixrAugPred.Rd
index 8e7179295..a04e9d211 100644
--- a/man/nlmixrAugPred.Rd
+++ b/man/nlmixrAugPred.Rd
@@ -5,12 +5,13 @@
 \alias{augPred.nlmixrFitData}
 \title{Augmented Prediction for nlmixr fit}
 \usage{
-nlmixrAugPred(object, ..., covsInterpolation = c("linear", "locf", "nocb",
-  "midpoint"), primary = NULL, minimum = NULL, maximum = NULL,
-  length.out = 51L)
+nlmixrAugPred(object, ..., covsInterpolation = c("linear", "locf",
+  "nocb", "midpoint"), primary = NULL, minimum = NULL,
+  maximum = NULL, length.out = 51L)
 
 \method{augPred}{nlmixrFitData}(object, primary = NULL,
-  minimum = min(primary), maximum = max(primary), length.out = 51, ...)
+  minimum = min(primary), maximum = max(primary), length.out = 51,
+  ...)
 }
 \arguments{
 \item{object}{Nlmixr fit object}
diff --git a/man/nlmixrBounds.focei.upper.lower.Rd b/man/nlmixrBounds.focei.upper.lower.Rd
index 51a5a859d..28bd2fc0e 100644
--- a/man/nlmixrBounds.focei.upper.lower.Rd
+++ b/man/nlmixrBounds.focei.upper.lower.Rd
@@ -4,7 +4,8 @@
 \alias{nlmixrBounds.focei.upper.lower}
 \title{Get upper/lower/names for THETAs}
 \usage{
-nlmixrBounds.focei.upper.lower(obj, type = c("upper", "lower", "name", "err"))
+nlmixrBounds.focei.upper.lower(obj, type = c("upper", "lower", "name",
+  "err"))
 }
 \arguments{
 \item{obj}{Bounds object}
diff --git a/man/nlmixrPred.Rd b/man/nlmixrPred.Rd
index 60f5f4546..fc4a63021 100644
--- a/man/nlmixrPred.Rd
+++ b/man/nlmixrPred.Rd
@@ -5,19 +5,7 @@
 \alias{predict.nlmixrFitData}
 \title{Predict a nlmixr solved system}
 \usage{
-nlmixrPred(object, params = NULL, events = NULL, inits = NULL,
-  scale = NULL, covs = NULL, method = c("liblsoda", "lsoda", "dop853"),
-  transitAbs = NULL, atol = 1e-06, rtol = 1e-04, maxsteps = 5000L,
-  hmin = 0L, hmax = NULL, hini = 0L, maxordn = 12L, maxords = 5L, ...,
-  cores, covsInterpolation = c("linear", "locf", "nocb", "midpoint"),
-  addCov = FALSE, matrix = FALSE, sigma = NULL, sigmaDf = NULL,
-  nCoresRV = 1L, sigmaIsChol = FALSE, nDisplayProgress = 10000L,
-  amountUnits = NA_character_, timeUnits = "hours", stiff, theta = NULL,
-  eta = NULL, addDosing = FALSE, updateObject = FALSE, doSolve = TRUE,
-  omega = NULL, omegaDf = NULL, omegaIsChol = FALSE, nSub = 1L,
-  thetaMat = NULL, thetaDf = NULL, thetaIsChol = FALSE, nStud = 1L,
-  dfSub = 0, dfObs = 0, returnType = c("data.frame", "rxSolve", "matrix"),
-  seed = NULL, nsim = NULL, ipred = FALSE)
+nlmixrPred(object, ..., ipred = FALSE)
 
 \method{predict}{nlmixrFitData}(object, ...)
 }
@@ -26,239 +14,15 @@ nlmixrPred(object, params = NULL, events = NULL, inits = NULL,
 with a RxODE model specification, or a string with a RxODE
 model specification.}
 
-\item{params}{a numeric named vector with values for every
-parameter in the ODE system; the names must correspond to the
-parameter identifiers used in the ODE specification;}
-
-\item{events}{an \code{eventTable} object describing the input
-(e.g., doses) to the dynamic system and observation sampling
-time points (see \code{\link{eventTable}});}
-
-\item{inits}{a vector of initial values of the state variables
-(e.g., amounts in each compartment), and the order in this
-vector must be the same as the state variables (e.g., PK/PD
-compartments);}
-
-\item{scale}{a numeric named vector with scaling for ode
-parameters of the system.  The names must correstond to the
-parameter identifiers in the ODE specification. Each of the
-ODE variables will be divided by the scaling factor.  For
-example \code{scale=(center=2)} will divide the center ODE
-variable by 2.}
-
-\item{covs}{a matrix or dataframe the same number of rows as the
-sampling points defined in the events \code{eventTable}.  This
-is for time-varying covariates.}
-
-\item{method}{The method for solving ODEs.  Currently this supports:
-
-\itemize{
-\item \code{"liblsoda"} thread safe lsoda.  This supports parallel
-           thread-based solving, and ignores user Jacobian specification.
-\item \code{"lsoda"} -- LSODA solver.  Does not support parallel thread-based
-      solving, but allows user Jacobian specification.
-\item \code{"dop853"} -- DOP853 solver.  Does not support parallel thread-based
-        solving nor user Jacobain specification
-}}
-
-\item{transitAbs}{boolean indicating if this is a transit
-compartment absorption}
-
-\item{atol}{a numeric absolute tolerance (1e-8 by default) used
-by the ODE solver to determine if a good solution has been
-achieved;  This is also used in the solved linear model to check
-if prior doses do not add anything to the solution.}
-
-\item{rtol}{a numeric relative tolerance (1e-6 by default) used
-by the ODE solver to determine if a good solution has been
-achieved. This is also used in the solved linear model to check
- if prior doses do not add anything to the solution.}
-
-\item{maxsteps}{maximum number of (internally defined) steps allowed
-during one call to the solver. (5000 by default)}
-
-\item{hmin}{The minimum absolute step size allowed. The default
-value is 0.}
-
-\item{hmax}{The maximum absolute step size allowed.  The default
-checks for the maximum difference in times in your sampling and
-events, and uses this value.  The value 0 is equivalent to
-infinite maximum absolute step size.}
-
-\item{hini}{The step size to be attempted on the first step. The
-default value is determined by the solver (when hini = 0)}
-
-\item{maxordn}{The maximum order to be allowed for the nonstiff
-(Adams) method.  The default is 12.  It can be between 1 and
-12.}
-
-\item{maxords}{The maximum order to be allowed for the stiff (BDF)
-method.  The default value is 5.  This can be between 1 and 5.}
-
 \item{...}{Other arguments including scaling factors for each
 compartment.  This includes S# = numeric will scale a compartment
 # by a dividing the compartment amount by the scale factor,
 like NONMEM.}
 
-\item{cores}{Number of cores used in parallel ODE solving.  This
-defaults to the number or system cores determined by
-\code{\link{rxCores}} for methods that support parallel
-solving (ie thread-safe methods like "liblsoda").}
-
-\item{covsInterpolation}{specifies the interpolation method for
-    time-varying covariates. When solving ODEs it often samples
-    times outside the sampling time specified in \code{events}.
-    When this happens, the time varying covariates are
-    interpolated.  Currently this can be:
-
-\itemize{
-\item \code{"linear"} interpolation (the default), which interpolates the covariate
-    by solving the line between the observed covariates and extrapolating the new
-    covariate value.
-\item \code{"constant"} -- Last observation carried forward.
-\item \code{"NOCB"} -- Next Observation Carried Backward.  This is the same method
-      that NONMEM uses.
-\item \code{"midpoint"} Last observation carried forward to midpoint; Next observation
-  carried backward to midpoint.
-}}
-
-\item{addCov}{A boolean indicating if covariates should be added
-to the output matrix or data frame. By default this is
-disabled.}
-
-\item{matrix}{A boolean inticating if a matrix should be returned
-instead of the RxODE's solved object.}
-
-\item{sigma}{Named sigma covariance or Cholesky decomposition of a
-covariance matrix.  The names of the columns indicate
-parameters that are simulated.  These are simulated for every
-observation in the solved system.}
-
-\item{sigmaDf}{Degrees of freedom of the sigma t-distribution.  By
-default it is equivalent to \code{Inf}, or a normal distribution.}
-
-\item{nCoresRV}{Number of cores used for the simulation of the
-sigma variables.  By default this is 1. This uses the package
-\code{\link[mvnfast]{rmvn}} and \code{\link[mvnfast]{rmvt}}.
-To reproduce the results you need to run on the same platform
-with the same number of cores. This is the reason this is set
-to be one, regardless of what the number of cores are used in
-threaded ODE solving.}
-
-\item{sigmaIsChol}{Boolean indicating if the sigma is in the
-Cholesky decomposition instead of a symmetric covariance}
-
-\item{nDisplayProgress}{An integer indicating the minimum number
-of c-based solves before a progress bar is shown.  By default
-this is 10,000.}
-
-\item{amountUnits}{This supplies the dose units of a data frame
-supplied instead of an event table.  This is for importing the
-data as an RxODE event table.}
-
-\item{timeUnits}{This supplies the time units of a data frame
-supplied instead of an event table.  This is for importing the
-data as an RxODE event table.}
-
-\item{stiff}{a logical (\code{TRUE} by default) indicating whether
-    the ODE system is stiff or not.
-
-    For stiff ODE sytems (\code{stiff = TRUE}), \code{RxODE} uses the
-    LSODA (Livermore Solver for Ordinary Differential Equations)
-    Fortran package, which implements an automatic method switching
-    for stiff and non-stiff problems along the integration
-    interval, authored by Hindmarsh and Petzold (2003).
-
-    For non-stiff systems (\code{stiff = FALSE}), \code{RxODE} uses
-    DOP853, an explicit Runge-Kutta method of order 8(5, 3) of
-    Dormand and Prince as implemented in C by Hairer and Wanner
-    (1993).}
-
-\item{theta}{A vector of parameters that will be named THETA[#] and
-added to parameters}
-
-\item{eta}{A vector of parameters that will be named ETA[#] and
-added to parameters}
-
-\item{addDosing}{Boolean indicating if the solve should add RxODE
-evid and amt columns.  This will also include dosing
-information and estimates at the doses.  Be default, RxODE
-only includes estimates at the observations. (default
-\code{FALSE}).}
-
-\item{updateObject}{This is an internally used flag to update the
-RxODE solved object (when supplying an RxODE solved object) as
-well as returning a new object.  You probably should not
-modify it's \code{FALSE} default unless you are willing to
-have unexpected results.}
-
-\item{doSolve}{Internal flag.  By default this is \code{TRUE},
-when \code{FALSE} a list of solving options is returned.}
-
-\item{omega}{Named omega matrix.}
-
-\item{omegaDf}{The degrees of freedom of a t-distribution for
-simulation.  By default this is \code{NULL} which is
-equivalent to \code{Inf} degrees, or to simulate from a normal
-distribution instead of a t-distribution.}
-
-\item{omegaIsChol}{Indicates if the \code{omega} supplied is a
-Cholesky decomposed matrix instead of the traditional
-symmetric matrix.}
-
-\item{nSub}{Number between subject variabilities (ETAs) simulated for every 
-realization of the parameters.}
-
-\item{thetaMat}{Named theta matrix.}
-
-\item{thetaDf}{The degrees of freedom of a t-distribution for
-simulation.  By default this is \code{NULL} which is
-equivalent to \code{Inf} degrees, or to simulate from a normal
-distribution instead of a t-distribution.}
-
-\item{thetaIsChol}{Indicates if the \code{theta} supplied is a
-Cholesky decomposed matrix instead of the traditional
-symmetric matrix.}
-
-\item{nStud}{Number virtual studies to characterize uncertainty in estimated 
-parameters.}
-
-\item{dfSub}{Degrees of freedom to sample the between subject variaiblity matrix from the 
-inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.}
-
-\item{dfObs}{Degrees of freedom to sample the unexplained variaiblity matrix from the 
-inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.}
-
-\item{returnType}{This tells what type of object is returned.  The currently supported types are:
-\itemize{
-\item \code{"rxSolve"} (default) will return a reactive data frame
-     that can change easily change different pieces of the solve and
-     update the data frame.  This is the currently standard solving
-     method in RxODE,  is used for \code{rxSolve(object, ...)}, \code{solev(object,...)},
-\item \code{"data.frame"} -- returns a plain, non-reactive data
-     frame; Currently very slightly Faster than \code{returnType=\"matrix\"}
-\item \code{"matrix"} -- returns a plain matrix with column names attached
-    to the solved object.  This is what is used \code{object$run} as well as \code{object$solve}
-}}
-
-\item{seed}{an object specifying if and how the random number
-    generator should be initialized (\sQuote{seeded}).\cr
-    For the "lm" method, either \code{NULL} or an integer that will be
-    used in a call to \code{set.seed} before simulating the response
-    vectors.  If set, the value is saved as the \code{"seed"} attribute
-    of the returned value.  The default, \code{NULL} will not change the
-    random generator state, and return \code{\link{.Random.seed}} as the
-    \code{"seed"} attribute, see \sQuote{Value}.
-  }
-
-\item{nsim}{represents the number of simulations.  For RxODE, if
-you supply single subject event tables (created with
-eventTable)}
-
 \item{ipred}{Flag to calculate individual predictions. When
 \code{ipred} is \code{TRUE}, calculate individual predictions.
 When \code{ipred} is \code{FALSE}, set calculate typical population predations.
-When \code{ipred} is \code{NA}, calculateboth individual and
+When \code{ipred} is \code{NA}, calculate both individual and
 population predictions.}
 }
 \description{
diff --git a/man/nlmixrSim.Rd b/man/nlmixrSim.Rd
index e8b7e1a40..82005d091 100644
--- a/man/nlmixrSim.Rd
+++ b/man/nlmixrSim.Rd
@@ -12,13 +12,14 @@ nlmixrSim(object, ...)
 \method{rxSolve}{nlmixrFitData}(object, params = NULL, events = NULL,
   inits = NULL, scale = NULL, covs = NULL, method = c("liblsoda",
   "lsoda", "dop853"), transitAbs = NULL, atol = 1e-06, rtol = 1e-04,
-  maxsteps = 5000L, hmin = 0L, hmax = NULL, hini = 0L, maxordn = 12L,
-  maxords = 5L, ..., cores, covsInterpolation = c("linear", "locf", "nocb",
-  "midpoint"), addCov = FALSE, matrix = FALSE, sigma = NULL,
-  sigmaDf = NULL, nCoresRV = 1L, sigmaIsChol = FALSE,
-  nDisplayProgress = 10000L, amountUnits = NA_character_,
-  timeUnits = "hours", stiff, theta = NULL, eta = NULL,
-  addDosing = FALSE, updateObject = FALSE, doSolve = TRUE, omega = NULL,
+  maxsteps = 5000L, hmin = 0L, hmax = NULL, hini = 0L,
+  maxordn = 12L, maxords = 5L, ..., cores,
+  covsInterpolation = c("linear", "locf", "nocb", "midpoint"),
+  addCov = FALSE, matrix = FALSE, sigma = NULL, sigmaDf = NULL,
+  nCoresRV = 1L, sigmaIsChol = FALSE, nDisplayProgress = 10000L,
+  amountUnits = NA_character_, timeUnits = "hours", stiff,
+  theta = NULL, eta = NULL, addDosing = FALSE,
+  updateObject = FALSE, doSolve = TRUE, omega = NULL,
   omegaDf = NULL, omegaIsChol = FALSE, nSub = 1L, thetaMat = NULL,
   thetaDf = NULL, thetaIsChol = FALSE, nStud = 1L, dfSub = 0,
   dfObs = 0, returnType = c("rxSolve", "matrix", "data.frame"),
@@ -29,246 +30,9 @@ nlmixrSim(object, ...)
 \method{solve}{nlmixrFitData}(a, b, ...)
 }
 \arguments{
-\item{object}{is a either a RxODE family of objects, or a file-name
-with a RxODE model specification, or a string with a RxODE
-model specification.}
+\item{object}{nlmixr object}
 
-\item{...}{Other arguments including scaling factors for each
-compartment.  This includes S# = numeric will scale a compartment
-# by a dividing the compartment amount by the scale factor,
-like NONMEM.}
-
-\item{params}{a numeric named vector with values for every
-parameter in the ODE system; the names must correspond to the
-parameter identifiers used in the ODE specification;}
-
-\item{events}{an \code{eventTable} object describing the input
-(e.g., doses) to the dynamic system and observation sampling
-time points (see \code{\link{eventTable}});}
-
-\item{inits}{a vector of initial values of the state variables
-(e.g., amounts in each compartment), and the order in this
-vector must be the same as the state variables (e.g., PK/PD
-compartments);}
-
-\item{scale}{a numeric named vector with scaling for ode
-parameters of the system.  The names must correstond to the
-parameter identifiers in the ODE specification. Each of the
-ODE variables will be divided by the scaling factor.  For
-example \code{scale=(center=2)} will divide the center ODE
-variable by 2.}
-
-\item{covs}{a matrix or dataframe the same number of rows as the
-sampling points defined in the events \code{eventTable}.  This
-is for time-varying covariates.}
-
-\item{method}{The method for solving ODEs.  Currently this supports:
-
-\itemize{
-\item \code{"liblsoda"} thread safe lsoda.  This supports parallel
-           thread-based solving, and ignores user Jacobian specification.
-\item \code{"lsoda"} -- LSODA solver.  Does not support parallel thread-based
-      solving, but allows user Jacobian specification.
-\item \code{"dop853"} -- DOP853 solver.  Does not support parallel thread-based
-        solving nor user Jacobain specification
-}}
-
-\item{transitAbs}{boolean indicating if this is a transit
-compartment absorption}
-
-\item{atol}{a numeric absolute tolerance (1e-8 by default) used
-by the ODE solver to determine if a good solution has been
-achieved;  This is also used in the solved linear model to check
-if prior doses do not add anything to the solution.}
-
-\item{rtol}{a numeric relative tolerance (1e-6 by default) used
-by the ODE solver to determine if a good solution has been
-achieved. This is also used in the solved linear model to check
- if prior doses do not add anything to the solution.}
-
-\item{maxsteps}{maximum number of (internally defined) steps allowed
-during one call to the solver. (5000 by default)}
-
-\item{hmin}{The minimum absolute step size allowed. The default
-value is 0.}
-
-\item{hmax}{The maximum absolute step size allowed.  The default
-checks for the maximum difference in times in your sampling and
-events, and uses this value.  The value 0 is equivalent to
-infinite maximum absolute step size.}
-
-\item{hini}{The step size to be attempted on the first step. The
-default value is determined by the solver (when hini = 0)}
-
-\item{maxordn}{The maximum order to be allowed for the nonstiff
-(Adams) method.  The default is 12.  It can be between 1 and
-12.}
-
-\item{maxords}{The maximum order to be allowed for the stiff (BDF)
-method.  The default value is 5.  This can be between 1 and 5.}
-
-\item{cores}{Number of cores used in parallel ODE solving.  This
-defaults to the number or system cores determined by
-\code{\link{rxCores}} for methods that support parallel
-solving (ie thread-safe methods like "liblsoda").}
-
-\item{covsInterpolation}{specifies the interpolation method for
-    time-varying covariates. When solving ODEs it often samples
-    times outside the sampling time specified in \code{events}.
-    When this happens, the time varying covariates are
-    interpolated.  Currently this can be:
-
-\itemize{
-\item \code{"linear"} interpolation (the default), which interpolates the covariate
-    by solving the line between the observed covariates and extrapolating the new
-    covariate value.
-\item \code{"constant"} -- Last observation carried forward.
-\item \code{"NOCB"} -- Next Observation Carried Backward.  This is the same method
-      that NONMEM uses.
-\item \code{"midpoint"} Last observation carried forward to midpoint; Next observation
-  carried backward to midpoint.
-}}
-
-\item{addCov}{A boolean indicating if covariates should be added
-to the output matrix or data frame. By default this is
-disabled.}
-
-\item{matrix}{A boolean inticating if a matrix should be returned
-instead of the RxODE's solved object.}
-
-\item{sigma}{Named sigma covariance or Cholesky decomposition of a
-covariance matrix.  The names of the columns indicate
-parameters that are simulated.  These are simulated for every
-observation in the solved system.}
-
-\item{sigmaDf}{Degrees of freedom of the sigma t-distribution.  By
-default it is equivalent to \code{Inf}, or a normal distribution.}
-
-\item{nCoresRV}{Number of cores used for the simulation of the
-sigma variables.  By default this is 1. This uses the package
-\code{\link[mvnfast]{rmvn}} and \code{\link[mvnfast]{rmvt}}.
-To reproduce the results you need to run on the same platform
-with the same number of cores. This is the reason this is set
-to be one, regardless of what the number of cores are used in
-threaded ODE solving.}
-
-\item{sigmaIsChol}{Boolean indicating if the sigma is in the
-Cholesky decomposition instead of a symmetric covariance}
-
-\item{nDisplayProgress}{An integer indicating the minimum number
-of c-based solves before a progress bar is shown.  By default
-this is 10,000.}
-
-\item{amountUnits}{This supplies the dose units of a data frame
-supplied instead of an event table.  This is for importing the
-data as an RxODE event table.}
-
-\item{timeUnits}{This supplies the time units of a data frame
-supplied instead of an event table.  This is for importing the
-data as an RxODE event table.}
-
-\item{stiff}{a logical (\code{TRUE} by default) indicating whether
-    the ODE system is stiff or not.
-
-    For stiff ODE sytems (\code{stiff = TRUE}), \code{RxODE} uses the
-    LSODA (Livermore Solver for Ordinary Differential Equations)
-    Fortran package, which implements an automatic method switching
-    for stiff and non-stiff problems along the integration
-    interval, authored by Hindmarsh and Petzold (2003).
-
-    For non-stiff systems (\code{stiff = FALSE}), \code{RxODE} uses
-    DOP853, an explicit Runge-Kutta method of order 8(5, 3) of
-    Dormand and Prince as implemented in C by Hairer and Wanner
-    (1993).}
-
-\item{theta}{A vector of parameters that will be named THETA[#] and
-added to parameters}
-
-\item{eta}{A vector of parameters that will be named ETA[#] and
-added to parameters}
-
-\item{addDosing}{Boolean indicating if the solve should add RxODE
-evid and amt columns.  This will also include dosing
-information and estimates at the doses.  Be default, RxODE
-only includes estimates at the observations. (default
-\code{FALSE}).}
-
-\item{updateObject}{This is an internally used flag to update the
-RxODE solved object (when supplying an RxODE solved object) as
-well as returning a new object.  You probably should not
-modify it's \code{FALSE} default unless you are willing to
-have unexpected results.}
-
-\item{doSolve}{Internal flag.  By default this is \code{TRUE},
-when \code{FALSE} a list of solving options is returned.}
-
-\item{omega}{Named omega matrix.}
-
-\item{omegaDf}{The degrees of freedom of a t-distribution for
-simulation.  By default this is \code{NULL} which is
-equivalent to \code{Inf} degrees, or to simulate from a normal
-distribution instead of a t-distribution.}
-
-\item{omegaIsChol}{Indicates if the \code{omega} supplied is a
-Cholesky decomposed matrix instead of the traditional
-symmetric matrix.}
-
-\item{nSub}{Number between subject variabilities (ETAs) simulated for every 
-realization of the parameters.}
-
-\item{thetaMat}{Named theta matrix.}
-
-\item{thetaDf}{The degrees of freedom of a t-distribution for
-simulation.  By default this is \code{NULL} which is
-equivalent to \code{Inf} degrees, or to simulate from a normal
-distribution instead of a t-distribution.}
-
-\item{thetaIsChol}{Indicates if the \code{theta} supplied is a
-Cholesky decomposed matrix instead of the traditional
-symmetric matrix.}
-
-\item{nStud}{Number virtual studies to characterize uncertainty in estimated 
-parameters.}
-
-\item{dfSub}{Degrees of freedom to sample the between subject variaiblity matrix from the 
-inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.}
-
-\item{dfObs}{Degrees of freedom to sample the unexplained variaiblity matrix from the 
-inverse Wishart distribution (scaled) or scaled inverse chi squared distribution.}
-
-\item{returnType}{This tells what type of object is returned.  The currently supported types are:
-\itemize{
-\item \code{"rxSolve"} (default) will return a reactive data frame
-     that can change easily change different pieces of the solve and
-     update the data frame.  This is the currently standard solving
-     method in RxODE,  is used for \code{rxSolve(object, ...)}, \code{solev(object,...)},
-\item \code{"data.frame"} -- returns a plain, non-reactive data
-     frame; Currently very slightly Faster than \code{returnType=\"matrix\"}
-\item \code{"matrix"} -- returns a plain matrix with column names attached
-    to the solved object.  This is what is used \code{object$run} as well as \code{object$solve}
-}}
-
-\item{seed}{an object specifying if and how the random number
-    generator should be initialized (\sQuote{seeded}).\cr
-    For the "lm" method, either \code{NULL} or an integer that will be
-    used in a call to \code{set.seed} before simulating the response
-    vectors.  If set, the value is saved as the \code{"seed"} attribute
-    of the returned value.  The default, \code{NULL} will not change the
-    random generator state, and return \code{\link{.Random.seed}} as the
-    \code{"seed"} attribute, see \sQuote{Value}.
-  }
-
-\item{nsim}{represents the number of simulations.  For RxODE, if
-you supply single subject event tables (created with
-eventTable)}
-
-\item{a}{when using \code{solve}, this is equivalent to the
-\code{object} argument.  If you specify \code{object} later in
-the argument list it overwrites this parameter.}
-
-\item{b}{when using \code{solve}, this is equivalent to the
-\code{params} argument.  If you specify \code{params} as a
-named argument, this overwrites the output}
+\item{...}{Other arguments sent to \code{rxSolve}}
 }
 \description{
 This takes the uncertainty in the model parameter estimates and to
diff --git a/man/nlmixr_fit.Rd b/man/nlmixr_fit.Rd
index d3654197e..c8acd21a4 100644
--- a/man/nlmixr_fit.Rd
+++ b/man/nlmixr_fit.Rd
@@ -4,8 +4,8 @@
 \alias{nlmixr_fit}
 \title{Fit a nlmixr model}
 \usage{
-nlmixr_fit(uif, data, est = NULL, control = list(), ..., sum.prod = FALSE,
-  table = tableControl())
+nlmixr_fit(uif, data, est = NULL, control = list(), ...,
+  sum.prod = FALSE, table = tableControl())
 }
 \arguments{
 \item{uif}{Parsed nlmixr model (by \code{nlmixr(mod.fn)}).}
diff --git a/man/reexports.Rd b/man/reexports.Rd
index 766b4d2c1..936c417d4 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -4,88 +4,46 @@
 \name{reexports}
 \alias{reexports}
 \alias{\%>\%}
-\alias{reexports}
 \alias{RxODE}
-\alias{reexports}
 \alias{rxSolve}
-\alias{reexports}
 \alias{rxClean}
-\alias{reexports}
 \alias{rxCat}
-\alias{reexports}
 \alias{nlme}
-\alias{reexports}
 \alias{ACF}
-\alias{reexports}
 \alias{VarCorr}
-\alias{reexports}
 \alias{getVarCov}
-\alias{reexports}
 \alias{augPred}
-\alias{reexports}
 \alias{fixef}
-\alias{reexports}
 \alias{fixed.effects}
-\alias{reexports}
 \alias{ranef}
-\alias{reexports}
 \alias{random.effects}
-\alias{reexports}
 \alias{nlmeControl}
-\alias{reexports}
 \alias{eventTable}
-\alias{reexports}
 \alias{add.dosing}
-\alias{reexports}
 \alias{add.sampling}
-\alias{reexports}
 \alias{rxPermissive}
-\alias{reexports}
 \alias{pdDiag}
-\alias{reexports}
 \alias{pdSymm}
-\alias{reexports}
 \alias{pdLogChol}
-\alias{reexports}
 \alias{pdIdent}
-\alias{reexports}
 \alias{pdCompSymm}
-\alias{reexports}
 \alias{pdBlocked}
-\alias{reexports}
 \alias{pdNatural}
-\alias{reexports}
 \alias{pdConstruct}
-\alias{reexports}
 \alias{pdFactor}
-\alias{reexports}
 \alias{pdMat}
-\alias{reexports}
 \alias{pdMatrix}
-\alias{reexports}
 \alias{reStruct}
-\alias{reexports}
 \alias{varWeights}
-\alias{reexports}
 \alias{varPower}
-\alias{reexports}
 \alias{varFixed}
-\alias{reexports}
 \alias{varFunc}
-\alias{reexports}
 \alias{varExp}
-\alias{reexports}
 \alias{varConstPower}
-\alias{reexports}
 \alias{varIdent}
-\alias{reexports}
 \alias{varComb}
-\alias{reexports}
 \alias{groupedData}
-\alias{reexports}
 \alias{getData}
-\alias{reexports}
-\alias{rxSolve}
 \title{Objects exported from other packages}
 \keyword{internal}
 \description{
diff --git a/man/saem.fit.Rd b/man/saem.fit.Rd
index d0b00e93f..b3144212a 100644
--- a/man/saem.fit.Rd
+++ b/man/saem.fit.Rd
@@ -10,10 +10,10 @@
 \alias{saem.fit.default}
 \title{Fit an SAEM model}
 \usage{
-saem.fit(model, data, inits, PKpars = NULL, pred = NULL, covars = NULL,
-  mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
-  ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1, transitAbs = 0),
-  distribution = c("normal", "poisson", "binomial"), seed = 99)
+saem.fit(model, data, inits, PKpars = NULL, pred = NULL,
+  covars = NULL, mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2,
+  2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1, transitAbs =
+  0), distribution = c("normal", "poisson", "binomial"), seed = 99)
 
 saem(model, data, inits, PKpars = NULL, pred = NULL, covars = NULL,
   mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
@@ -21,30 +21,34 @@ saem(model, data, inits, PKpars = NULL, pred = NULL, covars = NULL,
   distribution = c("normal", "poisson", "binomial"), seed = 99)
 
 \method{saem}{fit.nlmixr.ui.nlme}(model, data, inits, PKpars = NULL,
-  pred = NULL, covars = NULL, mcmc = list(niter = c(200, 300), nmc = 3, nu
-  = c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1,
-  transitAbs = 0), distribution = c("normal", "poisson", "binomial"),
+  pred = NULL, covars = NULL, mcmc = list(niter = c(200, 300), nmc =
+  3, nu = c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff =
+  1, transitAbs = 0), distribution = c("normal", "poisson", "binomial"),
   seed = 99)
 
-\method{saem}{fit.function}(model, data, inits, PKpars = NULL, pred = NULL,
-  covars = NULL, mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
-  ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1, transitAbs = 0),
-  distribution = c("normal", "poisson", "binomial"), seed = 99)
+\method{saem}{fit.function}(model, data, inits, PKpars = NULL,
+  pred = NULL, covars = NULL, mcmc = list(niter = c(200, 300), nmc =
+  3, nu = c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff =
+  1, transitAbs = 0), distribution = c("normal", "poisson", "binomial"),
+  seed = 99)
 
-\method{saem}{fit.nlmixrUI}(model, data, inits, PKpars = NULL, pred = NULL,
-  covars = NULL, mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
-  ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1, transitAbs = 0),
-  distribution = c("normal", "poisson", "binomial"), seed = 99)
+\method{saem}{fit.nlmixrUI}(model, data, inits, PKpars = NULL,
+  pred = NULL, covars = NULL, mcmc = list(niter = c(200, 300), nmc =
+  3, nu = c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff =
+  1, transitAbs = 0), distribution = c("normal", "poisson", "binomial"),
+  seed = 99)
 
-\method{saem}{fit.RxODE}(model, data, inits, PKpars = NULL, pred = NULL,
-  covars = NULL, mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
-  ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1, transitAbs = 0),
-  distribution = c("normal", "poisson", "binomial"), seed = 99)
+\method{saem}{fit.RxODE}(model, data, inits, PKpars = NULL,
+  pred = NULL, covars = NULL, mcmc = list(niter = c(200, 300), nmc =
+  3, nu = c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff =
+  1, transitAbs = 0), distribution = c("normal", "poisson", "binomial"),
+  seed = 99)
 
-\method{saem}{fit.default}(model, data, inits, PKpars = NULL, pred = NULL,
-  covars = NULL, mcmc = list(niter = c(200, 300), nmc = 3, nu = c(2, 2, 2)),
-  ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff = 1, transitAbs = 0),
-  distribution = c("normal", "poisson", "binomial"), seed = 99)
+\method{saem}{fit.default}(model, data, inits, PKpars = NULL,
+  pred = NULL, covars = NULL, mcmc = list(niter = c(200, 300), nmc =
+  3, nu = c(2, 2, 2)), ODEopt = list(atol = 1e-06, rtol = 1e-04, stiff =
+  1, transitAbs = 0), distribution = c("normal", "poisson", "binomial"),
+  seed = 99)
 }
 \arguments{
 \item{model}{an RxODE model or lincmt()}
diff --git a/man/saemControl.Rd b/man/saemControl.Rd
index b39956c3b..c2f98a345 100644
--- a/man/saemControl.Rd
+++ b/man/saemControl.Rd
@@ -4,9 +4,10 @@
 \alias{saemControl}
 \title{Control Options for SAEM}
 \usage{
-saemControl(seed = 99, nBurn = 200, nEm = 300, nmc = 3, nu = c(2, 2,
-  2), atol = 1e-06, rtol = 1e-04, stiff = TRUE, transitAbs = FALSE,
-  print = 1, covMethod = c("fim", "r,s", "r", "s"), ...)
+saemControl(seed = 99, nBurn = 200, nEm = 300, nmc = 3, nu = c(2,
+  2, 2), atol = 1e-06, rtol = 1e-04, stiff = TRUE,
+  transitAbs = FALSE, print = 1, trace = 0, covMethod = c("fim",
+  "r,s", "r", "s"), logLik = FALSE, ...)
 }
 \arguments{
 \item{seed}{Random Seed for SAEM step.  (Needs to be set for
@@ -84,6 +85,9 @@ anything is printed to the console.  By default, this is 1.}
 
  "" Does not calculate the covariance step.}
 
+\item{logLik}{boolean indicating that log-likelihood should be
+calculate by Gaussian quadrature.}
+
 \item{...}{Other arguments to control SAEM.}
 }
 \value{
diff --git a/man/vpc_ui.Rd b/man/vpc_ui.Rd
index c98313141..0d09006a1 100644
--- a/man/vpc_ui.Rd
+++ b/man/vpc_ui.Rd
@@ -8,12 +8,12 @@
 \title{VPC based on ui model}
 \usage{
 vpc_ui(fit, data = NULL, n = 100, bins = "jenks", n_bins = "auto",
-  bin_mid = "mean", show = NULL, stratify = NULL, pred_corr = FALSE,
-  pred_corr_lower_bnd = 0, pi = c(0.05, 0.95), ci = c(0.05, 0.95),
-  uloq = NULL, lloq = NULL, log_y = FALSE, log_y_min = 0.001,
-  xlab = NULL, ylab = NULL, title = NULL, smooth = TRUE,
-  vpc_theme = NULL, facet = "wrap", labeller = NULL, vpcdb = FALSE,
-  verbose = FALSE, ...)
+  bin_mid = "mean", show = NULL, stratify = NULL,
+  pred_corr = FALSE, pred_corr_lower_bnd = 0, pi = c(0.05, 0.95),
+  ci = c(0.05, 0.95), uloq = NULL, lloq = NULL, log_y = FALSE,
+  log_y_min = 0.001, xlab = NULL, ylab = NULL, title = NULL,
+  smooth = TRUE, vpc_theme = NULL, facet = "wrap", labeller = NULL,
+  vpcdb = FALSE, verbose = FALSE, ...)
 
 \method{vpc}{nlmixrFitData}(sim, ...)