knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE, out.width = "100%" ) ## These options cache the models and the model simulations in R ## To run the actual models on your system, take the save options off. nlmixrVersion <- as.character(packageVersion("nlmixr")); ## options(nlmixr.save=TRUE, ## nlmixr.save.dir=file.path(system.file(package="nlmixr"), nlmixrVersion)); ## if (!dir.exists(getOption("nlmixr.save.dir"))) ## dir.create(getOption("nlmixr.save.dir"))
broom
and broom.mixed
are packages that attempt to put standard
model outputs into data frames. nlmixr supports the tidy
and
glance
methods but does not support augment
at this time.
Using a model with a covariance term, the Phenobarbital model, we can explore the different types of output that is used in the tidy functions.
To explore this, first we run the model:
library(nlmixr) library(broom.mixed) pheno <- function() { # Pheno with covariance ini({ tcl <- log(0.008) # typical value of clearance tv <- log(0.6) # typical value of volume ## var(eta.cl) eta.cl + eta.v ~ c(1, 0.01, 1) ## cov(eta.cl, eta.v), var(eta.v) # interindividual variability on clearance and volume add.err <- 0.1 # residual variability }) model({ cl <- exp(tcl + eta.cl) # individual value of clearance v <- exp(tv + eta.v) # individual value of volume ke <- cl / v # elimination rate constant d/dt(A1) = - ke * A1 # model differential equation cp = A1 / v # concentration in plasma cp ~ add(add.err) # define error model }) } ## We will run it two ways to allow comparisons fit.s <- nlmixr(pheno, pheno_sd, "saem", control=list(logLik=TRUE, print=0), table=list(cwres=TRUE, npde=TRUE)) fit.f <- nlmixr(pheno, pheno_sd, "focei", control=list(print=0), table=list(cwres=TRUE, npde=TRUE))
Often in fitting data, you would want to glance
at the fit to see
how well it fits. In broom
, glance
will give a summary of the fit
metrics of goodness of fit:
glance(fit.s)
Note in nlmixr it is possible to have more than one fit metric (based
on different quadratures, FOCEi approximation etc). However, the
glance
only returns the fit metrics that are current.
If you wish you can set the objective function to the focei objective function (which was already calculated with CWRES).
setOfv(fit.s,"gauss3_1.6")
Now the glance gives the gauss3_1.6
values.
glance(fit.s)
Of course you can always change the type of objective function that nlmixr uses:
setOfv(fit.s,"FOCEi") # Setting objective function to focei
By setting it back to the SAEM default objective function of
FOCEi
, the glance(fit.s)
has the same values again:
glance(fit.s)
For convenience, you can do this while you glance
at the objects:
glance(fit.s, type="FOCEi")
You can also tidy the model estimates into a data frame with broom for processing. This can be useful when integrating into 3rd parting modeling packages. With a consistent parameter format, tasks for multiple types of models can be automated and applied.
The default function for this is tidy
, which when applied to the
fit
object provides the overall parameter information in a tidy
dataset:
tidy(fit.s)
Note by default these are the parameters that are actually estimated in nlmixr, not the back-transformed values in the table from the printout. Of course, with mu-referenced models, you may want to exponentiate some of the terms. The broom package allows you to apply exponentiation on all the parameters, that is:
## Transformation applied on every parameter tidy(fit.s, exponentiate=TRUE)
Note:, in accordance with the rest of the broom package, when the
parameters with the exponentiated, the standard errors are transformed
to an approximate standard error by the formula: $\textrm{se}(\exp(x))
\approx \exp(\textrm{model estimate}_x)\times \textrm{se}_x$. This
can be confusing because the confidence intervals (described later)
are using the actual standard error and back-transforming to the
exponentiated scale. This is the reason why the default for nlmixr's
broom
interface is exponentiate=FALSE
, that is:
tidy(fit.s, exponentiate=FALSE) ## No transformation applied
If you want, you can also use the parsed back-transformation that is
used in nlmixr tables (ie fit$parFixedDf
). Please note that this
uses the approximate back-transformation for standard errors on the
log-scaled back-transformed values.
This is done by:
## Transformation applied to log-scaled population parameters tidy(fit.s, exponentiate=NA)
Also note, at the time of this writing the default separator between
variables is .
, which doesn't work well with this model giving
cor__eta.v.eta.cl
. You can easily change this by:
options(broom.mixed.sep2="..") tidy(fit.s)
This gives an easier way to parse value: cor__eta.v..eta.cl
The default R method confint
works with nlmixr fit objects:
confint(fit.s)
This transforms the variables as described above. You can still use
the exponentiate
parameter to control the display of the confidence
interval:
confint(fit.s, exponentiate=FALSE)
However, broom has also implemented it own way to make these data a tidy dataset. The easiest way to get these values in a nlmixr dataset is to use:
tidy(fit.s, conf.level=0.9)
The confidence interval is on the scale specified by exponentiate
, by default the estimated scale.
If you want to have the confidence on the adaptive back-transformed scale, you would simply use the following:
tidy(fit.s, conf.level=0.9, exponentiate=NA)
tidy
The type of information that is extracted can be controlled by the effects
argument.
The fixed effect parameters can be extracted by effects="fixed"
tidy(fit.s, effects="fixed")
The random standard deviations can be extracted by
effects="ran_pars"
:
tidy(fit.s, effects="ran_pars")
The random values, or in NONMEM the ETAs, can be extracted by
effects="ran_vals"
or effects="random"
head(tidy(fit.s, effects="ran_vals"))
This duplicate method of running effects
is because the broom
package supports effects="random"
while the broom.mixed
package
supports effects="ran_vals"
.
Random coefficients are the population fixed effect parameter + the random effect parameter, possibly transformed to the correct scale.
In this case we can extract this information from a nlmixr fit object by:
head(tidy(fit.s, effects="ran_coef"))
This can also be changed by the exponentiate
argument:
head(tidy(fit.s, effects="ran_coef", exponentiate=NA)) head(tidy(fit.s, effects="ran_coef", exponentiate=TRUE))
As explained above, this standard format makes it easier for tidyverse
packages to interact with model information. An example of this is
piping the tidy information to dplyr to filter the effects and then to
the dotwhisker
package to plot the model parameter confidence
intervals.
options(broom.mixed.sep2=": ", broom.mixed.sep2=", ") library(ggplot2) library(dotwhisker) library(dplyr) fit.s %>% tidy(exponentiate=NA) %>% filter(effect=="fixed") %>% dwplot()
This allows easy creation of report ready tables in many formats including word.
Huxtable relies on the broom
implementation
library(huxtable) tbl <- huxreg('Phenobarbitol'=fit.s) tbl
You can also use huxtable
to compare runs:
huxreg('SAEM'=fit.s, 'FOCEi'=fit.f)
A word-based table can also be easily created with the tool:
library(officer) library(flextable) ft <- huxtable::as_flextable(tbl); read_docx() %>% flextable::body_add_flextable(ft) %>% print(target="pheno.docx")
Which produces the following word document.
Happy tidying!
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.