d_conc <- as.data.frame(Theoph) |> rename(time = Time, subject = Subject)
d_dose <- Theoph |> as.data.frame() |>
group_by(Subject) |>
summarise(dose = Dose[1] * Wt[1], weight = Wt[1], .groups = "drop") |>
rename(subject = Subject) |>
mutate(time = 0)
o_conc <- PKNCAconc(d_conc, conc ~ time | subject)
o_dose <- PKNCAdose(d_dose, dose ~ time | subject, route = "extravascular")
o_data <- PKNCAdata(o_conc, o_dose, intervals = data.frame(
start = 0, end = Inf,
auclast = TRUE, aucinf.obs = TRUE, cmax = TRUE,
half.life = TRUE, span.ratio = TRUE, adj.r.squared = TRUE,
aucpext.obs = TRUE
))
o_nca <- pk.nca(o_data)13 Post-Processing Results
13.1 Why post-process?
After pk.nca() you may need to:
- Extract results as a tidy data frame for further analysis
- Exclude specific results that failed QC (poor half-life fit, protocol deviation)
- Normalize parameters by dose (built-in
.dnvariants) - Summarize with custom statistics, respecting exclusions
- Report with the
business.*summary helper functions - Inspect which concentrations each λz regression used
Most of this operates on the PKNCAresults object without rerunning the analysis (requesting .dn parameters is the exception — they are computed by pk.nca()). The package’s own Post-Processing vignette walks the same workflow.
13.2 Example data and analysis
Every section below post-processes the same PKNCAresults object, built from the Theoph dataset. The interval requests the headline parameters plus the regression- and extrapolation-quality metrics that the exclusion demos rely on (span.ratio, adj.r.squared, aucpext.obs):
13.3 Extracting results
as.data.frame() returns the full tidy result including excluded rows:
df <- as.data.frame(o_nca)
# Columns:
# subject — group identifier
# start/end — interval boundaries
# PPTESTCD — parameter code
# PPORRES — numeric result
# exclude — NA if valid; character reason if excluded
df |>
filter(PPTESTCD == "auclast") |>
select(subject, PPORRES, exclude) |>
arrange(subject)# A tibble: 12 × 3
subject PPORRES exclude
<ord> <dbl> <chr>
1 6 71.7 <NA>
2 7 88.0 <NA>
3 8 86.8 <NA>
4 11 77.9 <NA>
5 3 95.9 <NA>
6 2 88.7 <NA>
7 4 103. <NA>
8 9 83.9 <NA>
9 12 115. <NA>
10 10 136. <NA>
11 1 147. <NA>
12 5 118. <NA> 13.4 Excluding results
13.4.1 Manual exclusion with exclude()
exclude() takes your results object plus a reason string and either a logical mask or a rule FUN. A manual mask excludes exactly the rows you flag — it does not propagate to other parameters. Rule-based exclusion (next section) does propagate: each rule excludes its affected parameters and the parameters that depend on them.
# Exclude subject 3's half-life manually (e.g. failed QC review)
o_nca_excl <- exclude(
o_nca,
reason = "half-life failed manual QC",
mask = as.data.frame(o_nca)$subject == "3" &
as.data.frame(o_nca)$PPTESTCD == "half.life"
)
# Only the masked half.life row is excluded
as.data.frame(o_nca_excl) |>
filter(subject == "3", !is.na(exclude)) |>
select(subject, PPTESTCD, PPORRES, exclude)# A tibble: 1 × 4
subject PPTESTCD PPORRES exclude
<ord> <chr> <dbl> <chr>
1 3 half.life 6.77 half-life failed manual QC13.4.2 Rule-based exclusion functions
PKNCA provides ready-made rule functions to pass to exclude(FUN=...). Unlike a manual mask, rule-based exclusion propagates: each rule excludes the parameters it affects (for the half-life rules, the whole λz family) and the parameters that depend on them, such as aucinf.obs and aucpext.obs. The built-in rules also supply their own reason text, which overrides the reason argument in the exclude column.
| Function | Excludes when… |
|---|---|
exclude_nca_span.ratio(min.span.ratio) |
span.ratio < min.span.ratio (default from PKNCA.options("min.span.ratio")) |
exclude_nca_min.hl.r.squared(min.hl.r.squared) |
half-life unadjusted R² < min.hl.r.squared (default from PKNCA.options("min.hl.r.squared")) |
exclude_nca_min.hl.adj.r.squared(min.hl.adj.r.squared = 0.9) |
half-life adjusted R² < min.hl.adj.r.squared (stricter: adjusted for the number of points) |
exclude_nca_max.aucinf.pext(max.aucinf.pext) |
%AUCextrap > max.aucinf.pext (default from PKNCA.options("max.aucinf.pext")) |
exclude_nca_tmax_0() |
Tmax at or before 0 (shortcut for exclude_nca_tmax_early(0)) |
exclude_nca_tmax_early(tmax_early = 0) |
Tmax at or before a threshold time (tmax_early, default 0) |
exclude_nca_count_conc_measured(min_count, exclude_param_pattern = ...) |
fewer than min_count measured concentrations; excludes the AUC-family parameters matching exclude_param_pattern |
exclude_nca_by_param(parameter, min_thr = NULL, max_thr = NULL, affected_parameters = parameter) (≥ 0.12.2) |
the value of parameter falls outside [min_thr, max_thr]; excludes affected_parameters (default: parameter itself) |
# Exclude half-lives where span ratio < 2
o_nca_span <- exclude(
o_nca,
reason = "span.ratio < 2",
FUN = exclude_nca_span.ratio(2)
)Loading required namespace: testthatas.data.frame(o_nca_span) |>
filter(PPTESTCD == "half.life") |>
select(subject, PPORRES, exclude) |>
arrange(subject)# A tibble: 12 × 3
subject PPORRES exclude
<ord> <dbl> <chr>
1 6 7.89 <NA>
2 7 7.85 <NA>
3 8 8.51 <NA>
4 11 7.26 <NA>
5 3 6.77 <NA>
6 2 6.66 <NA>
7 4 6.98 <NA>
8 9 8.41 span.ratio < 2
9 12 6.29 <NA>
10 10 9.25 span.ratio < 2
11 1 14.3 span.ratio < 2
12 5 8.00 <NA> # Propagation: subject 9 fails the span-ratio rule, so the entire lambda.z
# family and the parameters that depend on it are excluded together
as.data.frame(o_nca_span) |>
filter(subject == "9", !is.na(exclude)) |>
select(subject, PPTESTCD, exclude)# A tibble: 12 × 3
subject PPTESTCD exclude
<ord> <chr> <chr>
1 9 lambda.z span.ratio < 2
2 9 r.squared span.ratio < 2
3 9 adj.r.squared span.ratio < 2
4 9 lambda.z.corrxy span.ratio < 2
5 9 lambda.z.time.first span.ratio < 2
6 9 lambda.z.time.last span.ratio < 2
7 9 lambda.z.n.points span.ratio < 2
8 9 clast.pred span.ratio < 2
9 9 half.life span.ratio < 2
10 9 span.ratio span.ratio < 2
11 9 aucinf.obs span.ratio < 2
12 9 aucpext.obs span.ratio < 2# Chain multiple exclusion rules
o_nca_multi <- o_nca |>
exclude(reason = "span.ratio < 2", FUN = exclude_nca_span.ratio(2)) |>
exclude(reason = "%AUCextrap > 20%", FUN = exclude_nca_max.aucinf.pext(20))
as.data.frame(o_nca_multi) |>
filter(!is.na(exclude)) |>
select(subject, PPTESTCD, exclude) |>
arrange(subject, PPTESTCD)# A tibble: 36 × 3
subject PPTESTCD exclude
<ord> <chr> <chr>
1 9 adj.r.squared span.ratio < 2
2 9 aucinf.obs span.ratio < 2
3 9 aucpext.obs span.ratio < 2
4 9 clast.pred span.ratio < 2
5 9 half.life span.ratio < 2
6 9 lambda.z span.ratio < 2
7 9 lambda.z.corrxy span.ratio < 2
8 9 lambda.z.n.points span.ratio < 2
9 9 lambda.z.time.first span.ratio < 2
10 9 lambda.z.time.last span.ratio < 2
# ℹ 26 more rows13.4.3 Exclusions are respected in summary
Excluded values are dropped from the summary statistics: the parameters affected by the span-ratio rule are summarized over fewer subjects, and the caption defines the resulting n= annotation as the number of measurements included in the summary.
summary(o_nca_span) start end N auclast cmax half.life adj.r.squared
0 Inf 12 98.7 [22.5] 8.65 [17.0] 7.36 [0.742], n=9 0.997 [0.00331], n=9
span.ratio aucinf.obs aucpext.obs
2.34 [9.28], n=9 105 [16.4], n=9 11.3 [2.28], n=9
Caption: auclast, cmax, span.ratio, aucinf.obs: geometric mean and geometric coefficient of variation; half.life, adj.r.squared, aucpext.obs: arithmetic mean and standard deviation; N: number of subjects; n: number of measurements included in summary13.5 Dose-normalized parameters
PKNCA has built-in dose-normalized variants for several parameters. Request them in the interval — they are automatically computed as parameter / dose.
All available dose-normalized parameters, retrieved dynamically from the interval-column registry (get.interval.cols()):
cols <- get.interval.cols()
dn_params <- names(cols)[endsWith(names(cols), ".dn")]
data.frame(
parameter = dn_params,
description = sapply(cols[dn_params], function(x) x$desc),
normalized_from = sub("\\.dn$", "", dn_params)
) |> knitr::kable()| parameter | description | normalized_from | |
|---|---|---|---|
| auclast.dn | auclast.dn | Dose normalized auclast | auclast |
| aucall.dn | aucall.dn | Dose normalized aucall | aucall |
| aucinf.obs.dn | aucinf.obs.dn | Dose normalized aucinf.obs | aucinf.obs |
| aucinf.pred.dn | aucinf.pred.dn | Dose normalized aucinf.pred | aucinf.pred |
| aumclast.dn | aumclast.dn | Dose normalized aumclast | aumclast |
| aumcall.dn | aumcall.dn | Dose normalized aumcall | aumcall |
| aumcinf.obs.dn | aumcinf.obs.dn | Dose normalized aumcinf.obs | aumcinf.obs |
| aumcinf.pred.dn | aumcinf.pred.dn | Dose normalized aumcinf.pred | aumcinf.pred |
| cmax.dn | cmax.dn | Dose normalized cmax | cmax |
| cmin.dn | cmin.dn | Dose normalized cmin | cmin |
| clast.obs.dn | clast.obs.dn | Dose normalized clast.obs | clast.obs |
| clast.pred.dn | clast.pred.dn | Dose normalized clast.pred | clast.pred |
| cav.dn | cav.dn | Dose normalized cav | cav |
| ctrough.dn | ctrough.dn | Dose normalized ctrough | ctrough |
| clr.last.dn | clr.last.dn | Dose normalized clr.last | clr.last |
| clr.obs.dn | clr.obs.dn | Dose normalized clr.obs | clr.obs |
| clr.pred.dn | clr.pred.dn | Dose normalized clr.pred | clr.pred |
Each .dn parameter is computed as parameter / dose where dose is the dose from the PKNCAdose object for that interval. In PKNCA ≥ 0.12.2 this includes renal clearance variants (clr.last.dn, clr.obs.dn, clr.pred.dn).
Request the .dn variants in the interval like any other parameter:
dn_interval <- data.frame(
start = 0,
end = Inf,
auclast = TRUE,
auclast.dn = TRUE, # AUClast / dose
aucinf.obs = TRUE,
aucinf.obs.dn = TRUE, # AUCinf.obs / dose
aucinf.pred.dn = TRUE, # AUCinf.pred / dose
cmax = TRUE,
cmax.dn = TRUE, # Cmax / dose
cmin.dn = TRUE, # Cmin / dose
aumclast.dn = TRUE, # AUMClast / dose
clast.obs.dn = TRUE, # Clast.obs / dose
cav.dn = TRUE, # Cav / dose
half.life = TRUE # needed for aucinf.pred
)
o_nca_dn <- pk.nca(PKNCAdata(o_conc, o_dose, intervals = dn_interval))
as.data.frame(o_nca_dn) |>
filter(PPTESTCD %in% c("auclast", "auclast.dn", "cmax", "cmax.dn",
"aucinf.obs", "aucinf.obs.dn")) |>
select(subject, PPTESTCD, PPORRES) |>
tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
arrange(subject)# A tibble: 12 × 7
subject auclast cmax aucinf.obs auclast.dn aucinf.obs.dn cmax.dn
<ord> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 6 71.7 6.44 82.2 0.224 0.257 0.0201
2 7 88.0 7.09 101. 0.275 0.316 0.0222
3 8 86.8 7.56 102. 0.272 0.320 0.0237
4 11 77.9 8 86.9 0.244 0.272 0.0250
5 3 95.9 8.2 106. 0.300 0.332 0.0257
6 2 88.7 8.33 97.4 0.279 0.306 0.0261
7 4 103. 8.6 114. 0.321 0.357 0.0269
8 9 83.9 9.03 97.5 0.313 0.364 0.0337
9 12 115. 9.75 126. 0.359 0.392 0.0304
10 10 136. 10.2 168. 0.424 0.524 0.0319
11 1 147. 10.5 215. 0.460 0.672 0.0328
12 5 118. 11.4 136. 0.369 0.426 0.035613.5.1 Custom normalization with normalize() (≥ 0.12.2)
For normalization by columns other than dose (e.g. body weight, body surface area), PKNCA 0.12.2 adds normalize(). You supply a norm_table data frame containing the grouping columns, a normalization value, and a unit for each group:
# Normalize AUClast and Cmax by body weight
norm_tbl <- as.data.frame(Theoph) |>
group_by(Subject) |>
summarise(normalization = Wt[1], unit = "kg", .groups = "drop") |>
rename(subject = Subject)
df_results <- as.data.frame(o_nca)
df_bw_norm <- normalize(
df_results,
norm_table = norm_tbl,
parameters = c("auclast", "cmax"),
suffix = ".bw"
)
df_bw_norm |>
select(subject, PPTESTCD, PPORRES) |>
arrange(subject, PPTESTCD) |>
head(8)# A tibble: 8 × 3
subject PPTESTCD PPORRES
<ord> <chr> <dbl>
1 6 auclast.bw 0.896
2 6 cmax.bw 0.0805
3 7 auclast.bw 1.36
4 7 cmax.bw 0.110
5 8 auclast.bw 1.23
6 8 cmax.bw 0.107
7 11 auclast.bw 1.20
8 11 cmax.bw 0.123
normalize()returns only the normalized rows, under new parameter names —auclastbecomesauclast.bw,cmaxbecomescmax.bw. To keep both original and normalized values, bind the normalized rows onto the original results:bind_rows(df_results, df_bw_norm).
normalize_by_col() (≥ 0.12.2) is the low-level internal helper that pulls normalization values from a column already present in the PKNCAconc data frame, rather than from a separate table.
13.6 Custom summary statistics
PKNCA.set.summary() controls how each parameter is presented in summary().
# Arithmetic mean ± SD for AUClast
PKNCA.set.summary(
"auclast",
description = "mean ± SD",
point = mean,
spread = sd
)
# Median [min, max] for Tmax
# spread must return a 2-element numeric vector
PKNCA.set.summary(
"tmax",
description = "median [min, max]",
point = median,
spread = function(x) c(min(x), max(x))
)
summary(o_nca) start end N auclast cmax half.life adj.r.squared span.ratio
0 Inf 12 101 [23.5] 8.65 [17.0] 8.18 [2.12] 0.998 [0.00301] 2.07 [26.2]
aucinf.obs aucpext.obs
115 [28.4] 13.8 [6.34]
Caption: auclast: mean ± SD; cmax, span.ratio, aucinf.obs: geometric mean and geometric coefficient of variation; half.life, adj.r.squared, aucpext.obs: arithmetic mean and standard deviation; N: number of subjectsRestore the defaults — business.geomean/business.geocv for auclast and business.median/business.range for tmax (the business.* helpers are covered in the next section):
PKNCA.set.summary(
"auclast",
description = "geometric mean and geometric coefficient of variation",
point = business.geomean,
spread = business.geocv
)
PKNCA.set.summary(
"tmax",
description = "median and range",
point = business.median,
spread = business.range
)13.7 Business summary helpers
PKNCA exports the full set of business.* summary functions. These are the defaults summary() uses for most parameters, and they are also available directly for custom reporting pipelines.
| Function | Returns | Description |
|---|---|---|
business.geomean(x) |
scalar | Geometric mean: exp(mean(log(x))) |
business.geocv(x) |
scalar | Geometric CV%: sqrt(exp(var(log(x)))−1) × 100 |
business.mean(x) |
scalar | Arithmetic mean |
business.sd(x) |
scalar | Standard deviation |
business.cv(x) |
scalar | Coefficient of variation %: sd/mean × 100 |
business.median(x) |
scalar | Median |
business.min(x) |
scalar | Minimum |
business.max(x) |
scalar | Maximum |
business.range(x) |
length-2 vector | c(min, max) |
x <- c(10.2, 11.8, 9.6, 12.1, 10.5)
data.frame(
statistic = c("geomean", "geocv", "mean", "sd", "cv%", "median", "min", "max"),
value = c(
business.geomean(x), business.geocv(x),
business.mean(x), business.sd(x), business.cv(x),
business.median(x), business.min(x), business.max(x)
)
) statistic value
1 geomean 10.798057
2 geocv 9.859166
3 mean 10.840000
4 sd 1.069112
5 cv% 9.862655
6 median 10.500000
7 min 9.600000
8 max 12.100000Use these as point and spread arguments to PKNCA.set.summary() — rerunning summary() shows the new cmax presentation:
# Arithmetic mean ± CV% for Cmax
PKNCA.set.summary(
"cmax",
description = "mean [CV%]",
point = business.mean,
spread = business.cv
)
summary(o_nca) start end N auclast cmax half.life adj.r.squared span.ratio
0 Inf 12 98.7 [22.5] 8.76 [16.8] 8.18 [2.12] 0.998 [0.00301] 2.07 [26.2]
aucinf.obs aucpext.obs
115 [28.4] 13.8 [6.34]
Caption: auclast, span.ratio, aucinf.obs: geometric mean and geometric coefficient of variation; cmax: mean [CV%]; half.life, adj.r.squared, aucpext.obs: arithmetic mean and standard deviation; N: number of subjects13.8 Inspecting the terminal regression window
get_halflife_points() returns a logical vector aligned with the original concentration data — TRUE if that observation was used in the λz regression, FALSE if it was available but excluded, NA if half-life was not computed for that interval. In PKNCA ≥ 0.12.2 it also accepts a PKNCAdata object directly and handles start ≠ 0 intervals correctly. For how PKNCA selects the regression points, see the Half-Life Calculation vignette.
Because the example analysis requested half.life, o_nca already carries the λz fits:
# Logical vector (TRUE = used in lambda.z fit), aligned with the rows of d_conc
hl_used <- get_halflife_points(o_nca)
table(hl_used, useNA = "always")hl_used
FALSE TRUE <NA>
86 46 0 # >= 0.12.2: a PKNCAdata object is accepted directly (no need to run pk.nca first)
head(get_halflife_points(o_data))[1] FALSE FALSE FALSE FALSE FALSE FALSE# Attach to concentration data for plotting
d_conc_annotated <- d_conc |> mutate(in_lambda_z = hl_used)
subj1_conc <- d_conc_annotated |> filter(subject == "1")
ggplot(subj1_conc, aes(x = time, y = conc, colour = in_lambda_z)) +
geom_line(colour = "grey70") +
geom_point(size = 3) +
scale_colour_manual(values = c("FALSE" = "grey50", "TRUE" = "firebrick"),
labels = c("not used", "used in lambda.z"), name = NULL,
na.value = "grey80") +
scale_y_log10() +
labs(title = "Subject 1 - points used in lambda.z regression (red)",
x = "Time (h)", y = "Concentration (mg/L) [log scale]") +
theme_minimal()
13.9 Imputation methods
Imputation is a pre-calculation step, not a post-processing step. It is specified before pk.nca() runs via the impute argument to PKNCAdata() — either a method name applied to every interval, or the name of a column in the intervals data frame holding per-interval methods. See the Concentration Imputation chapter for the built-in strategies, chaining, custom methods, and a side-by-side comparison.