Generating Parameter Estimate Tables
parameter_tables.RmdOnce an NLME run has been imported into an xpose_data
object, most reporting tasks come down to pulling parameter estimates
and fit statistics back out in a tidy, table-ready form.
Certara.Xpose.NLME provides a family of get_*
extractors for a single model and
compare_prmNlme() for lining several
models up side by side.
| Function | Scope | Returns |
|---|---|---|
get_prmNlme() |
one model | raw parameter table (value, SE, %RSE, optional CI) |
get_overallNlme() |
one model | overall fit statistics (-2LL, AIC, BIC, condition
number, …) |
get_summaryNlme() |
one model | publication-style summary with transforms, %RSE, and shrinkage |
compare_prmNlme() |
many models | wide comparison table + run diagnostics |
To have something concrete to tabulate, we first fit two candidate
models – a one- and a two-compartment PK model – to the built-in
Certara.RsNLME::pkData dataset (16 subjects, single bolus
dose).
Fit two candidate models
We give each model its own working directory under the current
directory so that its NLME output files are kept separate. That working
directory is what we later hand to
xposeNlme(dir = ...).
pkData <- Certara.RsNLME::pkData
runs_dir <- file.path(getwd(), "runs")
# One-compartment model
oneCptModel <- pkmodel(
numCompartments = 1,
data = pkData,
ID = "Subject", Time = "Act_Time", A1 = "Amount", CObs = "Conc",
modelName = "OneCpt",
workingDir = file.path(runs_dir, "OneCpt")
) %>%
fixedEffect(effect = c("tvV", "tvCl"), value = c(15, 5)) %>%
randomEffect(effect = c("nV", "nCl"), value = c(0.1, 0.1)) %>%
residualError(predName = "C", SD = 0.2)
# Two-compartment model
twoCptModel <- pkmodel(
numCompartments = 2,
data = pkData,
ID = "Subject", Time = "Act_Time", A1 = "Amount", CObs = "Conc",
modelName = "TwoCpt",
workingDir = file.path(runs_dir, "TwoCpt")
) %>%
structuralParameter(paramName = "V2", hasRandomEffect = FALSE) %>%
fixedEffect(
effect = c("tvV", "tvCl", "tvV2", "tvCl2"),
value = c(15, 5, 40, 15)
) %>%
randomEffect(effect = c("nV", "nCl", "nCl2"), value = rep(0.1, 3)) %>%
residualError(predName = "C", SD = 0.2)
# Fit both models
oneCptFit <- fitmodel(oneCptModel)
twoCptFit <- fitmodel(twoCptModel)Import each run with xposeNlme(dir = ...)
xposeNlme()’s dir argument points at the
directory holding a run’s output files. Because we set each model’s
workingDir explicitly, we can pass those paths
directly.
xp_one <- xposeNlme(dir = oneCptModel@modelInfo@workingDir)
xp_two <- xposeNlme(dir = twoCptModel@modelInfo@workingDir)
xp_one
#> overview:
#> - Software: phx/nlme 26.9.1.1
#> - Attached files (memory usage 15.6 Mb):
#> + obs tabs: $prob no.1: nlme
#> + sim tabs: <none>
#> + output files: ConvergenceData, Overall, prmTable
#> + special: nlme eta_subject (#1)
#> - gg_theme:
#> - xp_theme: (modified)
#> - Options: dir = NULL, quiet = TRUE
get_prmNlme(): the raw parameter table
get_prmNlme() returns one row per estimated parameter,
carrying the type (the / ome /
sig), the PML label, the point estimate, its standard
error, and the relative standard error (rse, a
fraction).
get_prmNlme(xp_one)
#> # A tibble: 5 × 12
#> type name label value se rse fixed diagonal m n `2.5% CI`
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <int> <int> <dbl>
#> 1 the THETA(… tvV 32.8 2.16 0.0658 FALSE NA 1 NA 28.6
#> 2 the THETA(… tvCl 5.07 0.467 0.0920 FALSE NA 2 NA 4.15
#> 3 sig SIGMA(… CEps 0.737 0.0355 0.0482 FALSE TRUE 1 1 0.667
#> 4 ome OMEGA(… nV 0.0157 0.0203 1.29 FALSE TRUE 1 1 -0.0245
#> 5 ome OMEGA(… nCl 0.0759 0.0380 0.502 FALSE TRUE 2 2 0.000437
#> # ℹ 1 more variable: `97.5% CI` <dbl>Use digits to control the significant figures and
level to append a confidence interval (computed from
Student’s t on the residual degrees of freedom):
get_prmNlme(xp_one, digits = 4, level = 0.95)
#> # A tibble: 5 × 12
#> type name label value se rse fixed diagonal m n `2.5% CI`
#> <chr> <chr> <chr> <dbl> <dbl> <dbl> <lgl> <lgl> <int> <int> <dbl>
#> 1 the THETA(… tvV 32.8 2.16 0.0658 FALSE NA 1 NA 28.6
#> 2 the THETA(… tvCl 5.07 0.466 0.0920 FALSE NA 2 NA 4.15
#> 3 sig SIGMA(… CEps 0.737 0.0355 0.0482 FALSE TRUE 1 1 0.667
#> 4 ome OMEGA(… nV 0.0157 0.0203 1.29 FALSE TRUE 1 1 -0.0245
#> 5 ome OMEGA(… nCl 0.0759 0.0380 0.501 FALSE TRUE 2 2 0.000437
#> # ℹ 1 more variable: `97.5% CI` <dbl>Off-diagonal OMEGA elements that are frozen at zero are
dropped by default; pass show_all = TRUE to keep them.
get_overallNlme(): fit statistics
get_overallNlme() returns the Overall.csv
fit-statistics row for the run: the return code, log-likelihood and
-2LL, information criteria, the number of parameters /
observations / subjects, and the condition number.
get_overallNlme(xp_one)
#> # A tibble: 1 × 10
#> RetCode logLik `-2LL` AIC BIC nParm nObs nSub Condition ConditionBasis
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
#> 1 1 -780. 1559. 1569. 1583. 5 112 16 10.6 Covariance (fix…The ConditionBasis column records which basis the
condition number reflects; see ?get_overallNlme for how to
force a specific basis with the conditionNumber
argument.
get_summaryNlme(): a publication-style single-model
summary
get_summaryNlme() fuses fixed effects, random effects,
residual error, and shrinkage into one tibble, and – unlike the raw
table – reports each parameter on an interpretable scale. By default
random effects are shown as CV% (lognormal_cv)
and the %RSE is propagated onto that scale with the delta
method.
get_summaryNlme(xp_one, digits = 3)
#> # A tibble: 5 × 5
#> Section Parameter Estimate `%RSE` `Shrinkage (%)`
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 Fixed effects tvV 32.8 6.58 NA
#> 2 Fixed effects tvCl 5.07 9.2 NA
#> 3 Residual error CEps (CV%) 73.7 4.82 5.61
#> 4 Random effects nV (CV%) 12.6 65.1 59.9
#> 5 Random effects nCl (CV%) 28.1 26 19.6The default residual transform is multiplicative_cv
(100 * sigma), which is only exact for a genuinely
proportional error model. Override any parameter’s scale via
transform; for example, to report the residual error as the
standard deviation the engine actually estimated:
get_summaryNlme(
xp_one,
transform = list(CEps = "raw"),
digits = 3
)
#> # A tibble: 5 × 5
#> Section Parameter Estimate `%RSE` `Shrinkage (%)`
#> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 Fixed effects tvV 32.8 6.58 NA
#> 2 Fixed effects tvCl 5.07 9.2 NA
#> 3 Residual error CEps 0.737 4.82 5.61
#> 4 Random effects nV (CV%) 12.6 65.1 59.9
#> 5 Random effects nCl (CV%) 28.1 26 19.6See ?get_summaryNlme for the full preset catalog, custom
transform specs (list(fn = ..., dfn = ..., name = ...)),
and the shrinkage options.
get_bootSummaryNlme() is the bootstrap companion – see the
Bootstrap parameter summaries vignette.
compare_prmNlme(): comparing several models
compare_prmNlme() is the multi-model sibling of
get_summaryNlme(). It lines parameter estimates up across
runs and prepends a block of run-level diagnostics (-2LL,
OFV diff, method, RetCode,
condition, condition basis, nSub,
nObs, and total runtime), which makes it a convenient
run-record for model selection.
Pass a named list of xpose_data
objects; the list names become the column headers.
compare_prmNlme(list(OneCpt = xp_one, TwoCpt = xp_two))
#> # A tibble: 17 × 3
#> Description OneCpt TwoCpt
#> <chr> <chr> <chr>
#> 1 -2LL "1559.478" 1265.590
#> 2 OFV diff "" -293.888
#> 3 method "FOCE-ELS" FOCE-ELS
#> 4 RetCode "1" 1
#> 5 condition "10.6234" 4.83716
#> 6 condition basis "Covariance (fixed effects)" Covariance (fixed effects)
#> 7 nSub "16" 16
#> 8 nObs "112" 112
#> 9 total runtime (sec) "0.078" 0.078
#> 10 tvV "32.8402 (6.6%)" 15.3978 (7.4%)
#> 11 tvCl "5.07167 (9.2%)" 6.61267 (10.9%)
#> 12 tvV2 NA 41.2019 (2.6%)
#> 13 tvCl2 NA 14.0301 (7.4%)
#> 14 nV "0.0156936 (129.1%)" 0.0694048 (35.9%)
#> 15 nCl "0.0758571 (50.2%)" 0.182197 (29.7%)
#> 16 nCl2 NA 0.0427782 (82.6%)
#> 17 CEps "0.737328 (4.8%)" 0.161251 (12.2%)A few things to note in the output:
-
OFV diffis blank for the first (reference) model and shows the change in-2LLfor the others – here the two-compartment model drops the objective function substantially. - Parameters that exist in only one model (
tvV2,tvCl2,nCl2) areNAfor the model that lacks them, so common parameters (tvV,tvCl, …) still line up on a single row. - The estimate cells carry the estimate and its
%RSEinline. Watch for parameters with a very large%RSE– exactly the kind of over-parameterization this table is meant to surface.
The total runtime (sec) row is engine-reported
CPU time: the sum of the runtime and
covtime rows stored in each xpdb’s summary (parsed from
nlme7engine.log at import). That can differ substantially
from the wall-clock elapsed time shown by print.rsnlme_fit
/ a fit object’s runTime, especially on multi-core runs.
These two demo models fit in a fraction of a second, so both read near
0; on a real analysis this row reflects the actual engine
CPU total.
RSE in its own rows
Set rse_separate = TRUE to move each %RSE
onto its own row (suffixed " (RSE)"). This makes it easier
to scan across models for parameters whose precision changes materially
between fits.
compare_prmNlme(
list(OneCpt = xp_one, TwoCpt = xp_two),
rse_separate = TRUE
)
#> # A tibble: 25 × 3
#> Description OneCpt TwoCpt
#> <chr> <chr> <chr>
#> 1 -2LL "1559.478" 1265.590
#> 2 OFV diff "" -293.888
#> 3 method "FOCE-ELS" FOCE-ELS
#> 4 RetCode "1" 1
#> 5 condition "10.6234" 4.83716
#> 6 condition basis "Covariance (fixed effects)" Covariance (fixed effects)
#> 7 nSub "16" 16
#> 8 nObs "112" 112
#> 9 total runtime (sec) "0.078" 0.078
#> 10 tvV "32.8402" 15.3978
#> # ℹ 15 more rowsTransforming the random effects
By default diagonal OMEGA values are shown on the
estimated variance scale. transform rescales them (diagonal
SIGMA is always left on the reported scale):
-
"sqrt_om2"– standard deviation,sqrt(omega^2) -
"sqrt_exp_om2_minus_1"– approximate CV,sqrt(exp(omega^2) - 1)
The %RSE is propagated onto the chosen scale with the
delta method, so for "sqrt_om2" the OMEGA %RSE
is exactly halved.
compare_prmNlme(
list(OneCpt = xp_one, TwoCpt = xp_two),
transform = "sqrt_om2",
rse_separate = TRUE
)
#> # A tibble: 25 × 3
#> Description OneCpt TwoCpt
#> <chr> <chr> <chr>
#> 1 -2LL "1559.478" 1265.590
#> 2 OFV diff "" -293.888
#> 3 method "FOCE-ELS" FOCE-ELS
#> 4 RetCode "1" 1
#> 5 condition "10.6234" 4.83716
#> 6 condition basis "Covariance (fixed effects)" Covariance (fixed effects)
#> 7 nSub "16" 16
#> 8 nObs "112" 112
#> 9 total runtime (sec) "0.078" 0.078
#> 10 tvV "32.8402" 15.3978
#> # ℹ 15 more rowsOther useful arguments: param_order = "alphabetical" to
sort labels within each section, drop_dOFV = TRUE to omit
the OFV diff row, and max_runs to cap how many
models are included.
Auto-detecting runs in a directory
For users who run many models, compare_prmNlme() can
also scan a directory instead of taking an explicit list. With no
x supplied it looks in dir for model files
(*.mdl or *.mmdl) whose same-named run output
folder exists, imports each completed run in alphanumeric order, and
includes them. Any run that started but did not produce usable output (a
missing nlme7engine.log, or an import error) is skipped and
recorded in a log file rather than aborting the table.
This mode expects a project directory laid out with one model file
per run (name each model so its output folder matches the model file)
alongside a matching <run>/ output folder:
project/
run001.mmdl run001/ (nlme7engine.log, dmp.txt, ...)
run002.mmdl run002/ (nlme7engine.log, dmp.txt, ...)
# Scan the project directory, include every completed run, and write the
# full table to CSV.
tbl <- compare_prmNlme(
dir = "path/to/project",
output_file = "TableofParameters.csv"
)
# Runs that failed to load are listed here:
readLines(file.path("path/to/project", "Table_log.txt"))output_file writes the complete table (all models) as a
CSV in either format = "column" (models as columns, the
default) or format = "row" (models as rows,
transposed).
Rendering as a formatted table
The object returned by compare_prmNlme() is a tibble, so
it prints and subsets like any other. For a formatted report table it
also has an as_flextable() method that draws separators
between the diagnostic, parameter, and RSE blocks and reports the number
of models in the footer. This method requires the flextable
and officer packages (both Suggests).
tbl <- compare_prmNlme(
list(OneCpt = xp_one, TwoCpt = xp_two)
)
flextable::as_flextable(tbl)Description |
OneCpt |
TwoCpt |
|---|---|---|
-2LL |
1559.478 |
1265.590 |
OFV diff |
-293.888 |
|
method |
FOCE-ELS |
FOCE-ELS |
RetCode |
1 |
1 |
condition |
10.6234 |
4.83716 |
condition basis |
Covariance (fixed effects) |
Covariance (fixed effects) |
nSub |
16 |
16 |
nObs |
112 |
112 |
total runtime (sec) |
0.078 |
0.078 |
tvV |
32.8402 (6.6%) |
15.3978 (7.4%) |
tvCl |
5.07167 (9.2%) |
6.61267 (10.9%) |
tvV2 |
41.2019 (2.6%) |
|
tvCl2 |
14.0301 (7.4%) |
|
nV |
0.0156936 (129.1%) |
0.0694048 (35.9%) |
nCl |
0.0758571 (50.2%) |
0.182197 (29.7%) |
nCl2 |
0.0427782 (82.6%) |
|
CEps |
0.737328 (4.8%) |
0.161251 (12.2%) |
nModels: 2 | ||
Pass format = "row" to as_flextable() (or
to compare_prmNlme()’s output_file CSV) to
transpose the layout so models become rows instead of columns.
See also
-
get_summaryNlme()– single-model publication summary. -
get_bootSummaryNlme()– bootstrap parameter summaries. -
xposeNlme()/xposeNlmeModel()– importing runs intoxpose_data.