Summary and recommendations

For sparse sampling, where each animal gives one or a few samples:

• Use the delta-method standard error for the AUC to infinity (the default). When little of the AUC is extrapolated, Yuan’s standard error, which treats \(\lambda_z\) as known, and the delta method both cover close to the nominal 95% (Yuan 94 to 96%). When about a fifth of the AUC is extrapolated, Yuan’s intervals are too narrow (81 to 89% coverage), and the delta method is closer to nominal (90 to 94%). Set sparse_lambda_z_se = "none" for Yuan’s method.
• Give the AUC a concentration at time 0. As for any PKNCA AUC, the sparse AUC is NA without one. Sacrificial designs rarely sample at the dose, so after an extravascular dose impute zero with impute = "start_conc0". After an IV bolus, request the IV parameters (aucivlast, aucivinf.obs): they back-extrapolate \(C_0\) from the mean profile and include its uncertainty.
• Treat the AUMC, and the MRT and Vss that come from it, with caution. The AUMC to infinity was biased low in most scenarios, and its confidence intervals under-covered in every scenario, even with the delta method (79 to 92%), most when elimination was slow.
• The half-life of the mean profile is not the typical animal’s half-life. It was too long when the terminal phase was well sampled (by up to 11%) and too short when the profile was cut short by BLQ samples or a two-compartment distribution (down to -10%). Sample long enough to see at least three terminal half-lives.
• Cmax of the mean profile reads high with few animals per time. Taking the maximum of noisy means near a broad peak overestimated Cmax by up to 10%.
• Geometric-mean profiles were evaluated and not adopted. They were no more biased and were as precise or modestly more precise, but they estimate a different quantity (the typical animal rather than the population mean curve), need a rule for BLQ samples that drives their BLQ results, and fall outside the variance methods for sparse AUCs, which are for arithmetic means; adopting them would change all sparse calculations for a modest gain.
• The bootstrap (sparse_bootstrap()) gives intervals for any parameter, but they are too narrow with few animals at each time. With three or four animals per time, its 95% intervals for the AUC to infinity covered 82 to 90% of the time. Prefer the analytical standard errors for the AUCs, and use the bootstrap for parameters and comparisons that have none, such as Cmax and ratios between groups. With 20 or more animals per time, its standard error was accurate and its intervals matched the delta method’s.
• With many animals, the bias of the sparse estimate dominates. BLQ samples and the extrapolation of the mean profile biased the AUC to infinity low (by about 7% with BLQ samples here), whatever the number of animals and whichever half-life was used, so with many animals every confidence interval missed the truth more often than its nominal level.
• The Tobit half-life was evaluated with BLQ samples and not adopted for the sparse AUC. The Tobit fit of the mean profile (hl_method = "tobit", with the LLOQ carried through the mean profile) gave nearly the same AUC to infinity as the log-linear fit, and it has no delta-method standard error. A (mixed-effects) Tobit fit of the individual samples was slightly less biased and more precise, but its standard error added to Yuan’s covered less well than the delta method with slow elimination (87 to 92% against 91 to 94%), and neither Tobit fit removed the bias that BLQ samples cause.

Introduction

With sparse sampling, each animal contributes only one or a few samples to a concentration-time profile, so there is no individual AUC to summarize. Instead, the AUC is estimated from the mean concentration at each time, and its standard error comes from the variability between animals at each time. Bailer (1988) gave the estimator and its variance for the AUC from time 0 to the last sampling time with serial sacrifice (one sample per animal). Nedelman and Jia (1998) extended it to designs where animals give more than one sample, with a Satterthwaite approximation for the degrees of freedom, and Holder (2001) gave an unbiased estimator of the covariance between times. Yuan (1993) extended Bailer’s method to the AUC from 0 to infinity by adding the extrapolation \(\bar{y}_n/\lambda_z\), treating the terminal elimination rate \(\lambda_z\) as known.

In practice, \(\lambda_z\) is estimated from the same data, so its uncertainty adds to that of the AUC to infinity. This vignette describes how PKNCA calculates the sparse AUC and AUMC to infinity with their standard errors and degrees of freedom, gives a delta-method extension of Yuan’s method that includes the uncertainty in \(\lambda_z\), extends both to IV bolus dosing, and compares the methods in a simulation-estimation study. A second study compares arithmetic and geometric mean profiles.

The sparse AUC to the last measurement

At each sampling time \(t_i\) (\(i = 0, \ldots, n\)), \(r_i\) animals are sampled and \(\bar{y}_i\) is the arithmetic mean concentration. \(r_{ij}\) animals are sampled at both \(t_i\) and \(t_j\) (\(r_{ii} = r_i\)). The linear trapezoidal AUC of the mean profile is a weighted sum of the means,

\[\widehat{AUC}_{0-t_{last}} = \sum_{i=0}^{n} w_i \bar{y}_i, \quad w_0 = \frac{t_1 - t_0}{2}, \quad w_i = \frac{t_{i+1} - t_{i-1}}{2}, \quad w_n = \frac{t_n - t_{n-1}}{2},\]

where \(t_n = t_{last}\) is the last time with a positive mean. Since the estimate is linear in the means, its variance is

\[\text{Var}\left(\widehat{AUC}_{0-t_{last}}\right) = \sum_{i}\sum_{j} w_i w_j \frac{r_{ij}\,\sigma_{ij}}{r_i r_j},\]

with \(\sigma_{ij}\) the covariance of the concentrations at \(t_i\) and \(t_j\) within an animal. PKNCA estimates \(\sigma_{ii}\) with the sample variance and \(\sigma_{ij}\) (\(i \neq j\)) with the unbiased estimator of Holder (2001), setting it to 0 when fewer than two animals are sampled at both times. With serial sacrifice, \(r_{ij} = 0\) for \(i \neq j\) and this is Bailer’s (1988) variance, \(\sum_i w_i^2 s_i^2/r_i\).

Only times up to \(t_{last}\) have weight. When more than half of the samples at a time are below the limit of quantification (BLQ), PKNCA sets that time’s mean to zero; times after the last positive mean are not part of \(AUC_{0-t_{last}}\), and they contribute nothing to its variance either.

As for any AUC in PKNCA, the area starts at the start of the interval, so a concentration is needed there, measured or imputed; without one, the sparse AUC is NA. Animals are rarely sampled at the dose in sacrificial designs, so after an extravascular dose the usual choice is to impute a concentration of zero with impute = "start_conc0". An imputed zero is known rather than measured, so it has weight in the estimate (as \(\bar{y}_0 = 0\)) but no variance.

Extrapolation to infinity (Yuan 1993)

With first-order terminal elimination, the area after \(t_{last}\) is \(\bar{y}_{n}/\lambda_z\), so

\[\widehat{AUC}_{0-\infty} = \sum_{i=0}^{n} w_i \bar{y}_i + \frac{\bar{y}_n}{\lambda_z}.\]

This is still linear in the means if \(\lambda_z\) is known, with the last weight increased by \(1/\lambda_z\). Yuan’s (1993) variance (his equation 4) is Bailer’s with that weight, \(g_n = w_n + 1/\lambda_z\):

\[\text{Var}\left(\widehat{AUC}_{0-\infty}\right) = \sum_{i}\sum_{j} g_i g_j \frac{r_{ij}\,\sigma_{ij}}{r_i r_j}, \qquad g_i = w_i + \frac{\delta_{i,n}}{\lambda_z}.\]

Written this way, the method applies to any sampling design, not only serial sacrifice.

The AUMC is the AUC of the moment curve \(t\,C(t)\). Its estimate uses the moment means \(\bar{m}_i = t_i \bar{y}_i\), and the extrapolation \(t_{last} C_{last}/\lambda_z + C_{last}/\lambda_z^2\) is \(\bar{m}_n\left(1/\lambda_z + 1/(t_{last}\lambda_z^2)\right)\), so the AUMC weights are \(g_i = w_i + \delta_{i,n}\left(1/\lambda_z + 1/(t_{last}\lambda_z^2)\right)\) applied to the moment means, with the covariances of the moments \(t_i t_j \sigma_{ij}\).

IV bolus: the back-extrapolated C0

After an IV bolus, the concentration at time 0 is not zero, and it cannot be measured. As for dense data, PKNCA estimates it from the mean profile with the methods of pk.calc.c0(): log-linear back-extrapolation from the first two means when they decline,

\[\hat{C}_0 = \bar{y}_1 \left(\frac{\bar{y}_1}{\bar{y}_2}\right)^{t_1/(t_2 - t_1)},\]

and otherwise \(\bar{y}_1\). The AUC adds the trapezoid from time 0 to \(t_1\), \(t_1(\hat{C}_0 + \bar{y}_1)/2\). \(\hat{C}_0\) is a smooth function of the first two means, so its uncertainty enters the variance through its gradient, as in the next section: the weights of \(\bar{y}_1\) and \(\bar{y}_2\) gain \((t_1/2)\,\partial \hat{C}_0/\partial \bar{y}_j\), with

\[\frac{\partial \hat{C}_0}{\partial \bar{y}_1} = \frac{\hat{C}_0\, t_2}{(t_2 - t_1)\,\bar{y}_1}, \qquad \frac{\partial \hat{C}_0}{\partial \bar{y}_2} = -\frac{\hat{C}_0\, t_1}{(t_2 - t_1)\,\bar{y}_2}.\]

With the linear trapezoidal rule, the AUMC from time 0 to \(t_1\) is \(t_1^2\bar{y}_1/2\) whatever \(C_0\) is, so \(C_0\) does not affect the AUMC. These are the sparse aucivlast, aucivall, and aucivinf.obs (and their AUMC equivalents), each with its standard error and degrees of freedom.

The uncertainty of lambda.z (delta method)

\(\lambda_z\) is estimated from the mean profile: PKNCA fits a line to the log mean concentrations at the half-life points, chosen automatically by the adjusted \(r^2\) (see vignette("v06-half-life-calculation")). For a given set of half-life points \(H\), the estimate is minus the least-squares slope,

\[\hat{\lambda}_z = -\sum_{j \in H} c_j \log \bar{y}_j, \qquad c_j = \frac{t_j - \bar{t}_H}{\sum_{k \in H}\left(t_k - \bar{t}_H\right)^2},\]

a smooth function of the means. So \(\widehat{AUC}_{0-\infty}\) is a smooth function of the means, and the delta method gives its variance from the same covariance of the means with the gradient in place of the weights:

\[g_j = \frac{\partial \widehat{AUC}_{0-\infty}}{\partial \bar{y}_j} = w_j + \frac{\delta_{j,n}}{\lambda_z} + \frac{\bar{y}_n}{\lambda_z^2}\,\frac{c_j}{\bar{y}_j}\,\mathbb{1}_{j \in H}.\]

The last term is the uncertainty in \(\lambda_z\). Because the gradient is taken with respect to all of the means together, the variance includes the covariance of \(\hat{\lambda}_z\) with \(\bar{y}_n\) and with \(\widehat{AUC}_{0-t_{last}}\). Those covariances matter: a high mean at an early half-life point raises the AUC to \(t_{last}\) but steepens the slope, which lowers the extrapolated area, so the two partly cancel. Adding a separately estimated variance of \(\hat{\lambda}_z\) to Yuan’s variance would miss that.

For the AUMC, with the moment means, \(\hat{\lambda}_z = -\sum_{j \in H} c_j \left(\log \bar{m}_j - \log t_j\right)\) and

\[g_j = w_j + \delta_{j,n}\left(\frac{1}{\lambda_z} + \frac{1}{t_{last}\lambda_z^2}\right) + \bar{m}_n\left(\frac{1}{\lambda_z^2} + \frac{2}{t_{last}\lambda_z^3}\right)\frac{c_j}{\bar{m}_j}\,\mathbb{1}_{j \in H}.\]

The method conditions on the selected half-life points: the variability from choosing them is not included. It is a first-order approximation, which is most accurate when the coefficient of variation of the means at the half-life points is small.

Degrees of freedom

Confidence intervals use a \(t\) distribution. The variance estimate is a quadratic form in the individual measurements \(y\), \(\hat{V} = y^T M y\), where \(M\) centers the measurements at each time and weights the products of each animal’s measurements by \(g_i g_j r_{ij}/(r_i r_j h_{ij})\), with \(h_{ij}\) the divisor of the Holder covariance. With \(\Omega\) the covariance of the measurements, Nedelman and Jia’s (1998) equation 6 gives the Satterthwaite degrees of freedom,

\[\nu = \frac{2\,E^2}{V}, \qquad E = \text{tr}(M\Omega), \qquad V = 2\,\text{tr}(M\Omega M\Omega).\]

With serial sacrifice, this is equation 6a of Nedelman, Gibiansky, and Lau (1995). PKNCA uses the same expression with the weights \(g\) for the AUC and AUMC to infinity, with or without the delta-method terms.

An alternative: lambda.z from the individual samples

Instead of the log means, \(\lambda_z\) could be estimated from every individual sample in the terminal phase, with its standard error from that regression. Samples below the limit of quantification are censored, so the regression is a Tobit regression on the log concentrations; when animals give more than one sample, a random intercept for each animal accounts for their correlation. The variance of \(\hat{\lambda}_z\) is then added to Yuan’s variance.

Two properties of this approach differ from the delta method. A line through the individual log concentrations follows the geometric-mean profile, while the AUC and \(C_{last}\) come from the arithmetic means; the two have the same terminal slope only when the between-animal coefficient of variation is the same at every terminal time. And the regression does not give the covariance of \(\hat{\lambda}_z\) with the means, so it is added as though it were independent. The simulation below compares both.

Simulation-estimation study

Arithmetic or geometric mean profiles

The sparse AUC uses the arithmetic mean at each time, whose profile is the population mean concentration curve. A geometric mean profile instead approximates the typical animal, and it may be less affected by animals with high concentrations. The second study compared the two for the AUClast, AUCinf, Cmax (\(C_0\) for IV bolus), AUMClast, AUMCinf, and half-life, using the same design with 2,000 replicates per scenario (data-raw/sparse_mean_type_simulation.R).

For each replicate, both mean profiles were formed, and pk.nca() calculated the parameters from each with the linear trapezoidal rule (imputing zero at time 0 for oral dosing and back-extrapolating \(C_0\) for IV bolus). For the geometric mean, a time with more than half of its samples BLQ was BLQ (as for the arithmetic sparse mean), and otherwise BLQ samples counted as half of the LLOQ. Each profile was compared with its own target (the arithmetic or geometric mean of the individual parameters of 50,000 virtual animals) and with the typical animal. Precision is the robust coefficient of variation (the interquartile range divided by 1.349, relative to the median).

Findings. Against their own targets, the two profiles were similarly accurate for the AUCs and Cmax. The geometric profile was as precise or more precise, most clearly for the half-life with BLQ samples and for the AUMC to infinity, but the gains were modest elsewhere. Both underestimated the AUMC to infinity (arithmetic -30 to -6%, geometric -19 to 1%). The half-life of the arithmetic mean profile was not consistently too short: it was too long when the terminal phase was well sampled, because the late mean is dominated by the slowly eliminating animals, and too short with BLQ samples (counted as zero in the late means) and in the two-compartment model with slow elimination, where the fitted points reach back toward the distribution phase. The geometric profile was generally a few percent short, because the geometric mean of exponentials declines at the average elimination rate.

The geometric profile was not adopted. Its gains were modest; it estimates a different quantity (the typical animal rather than the population mean curve); its BLQ rule (here, half the LLOQ) is a choice that drives its results with BLQ samples; and the sparse variance methods (Bailer, Nedelman and Jia, Holder) are for arithmetic means, so adopting it would change every sparse calculation.

Calculating it with PKNCA

Request aucinf.obs (or aumcinf.obs) with sparse data. The standard error and degrees of freedom are reported with it as aucinf.obs_se and aucinf.obs_df, and the summary shows the estimate with its standard error. These animals were not sampled at the dose, so the concentration at time 0 is imputed as zero.

d_sparse <-
  data.frame(
    time = rep(c(0.5, 1, 2, 4, 6, 8, 12, 24), each = 4),
    conc =
      c(3.1, 4.4, 3.6, 5.2,  5.8, 6.9, 5.1, 7.4,  7.2, 6.1, 8.3, 6.6,
        5.0, 6.2, 4.4, 5.6,  4.1, 3.4, 4.8, 3.7,  2.9, 3.5, 2.4, 3.1,
        1.6, 1.9, 1.3, 1.8,  0.33, 0.27, 0.41, 0.30)
  )
d_sparse$animal <- seq_len(nrow(d_sparse))
o_data <-
  PKNCAdata(
    PKNCAconc(d_sparse, conc ~ time | animal, sparse = TRUE),
    intervals = data.frame(start = 0, end = Inf, aucinf.obs = TRUE, aumcinf.obs = TRUE),
    impute = "start_conc0"
  )
o_nca <- pk.nca(o_data)
#> The sparse estimators use the linear trapezoidal rule, so the auc.method option
#> ("lin up/log down") does not apply to: aucinf.obs, aumcinf.obs
d_result <- as.data.frame(o_nca)
d_result[grepl("^au.*inf|^lambda.z$", d_result$PPTESTCD), c("PPTESTCD", "PPORRES")]
#> # A tibble: 7 × 2
#>   PPTESTCD       PPORRES
#>   <chr>            <dbl>
#> 1 lambda.z         0.139
#> 2 aucinf.obs      62.4  
#> 3 aucinf.obs_se    1.85 
#> 4 aucinf.obs_df   15.6  
#> 5 aumcinf.obs    468.   
#> 6 aumcinf.obs_se  20.2  
#> 7 aumcinf.obs_df   7.14
summary(o_nca)
#>  start end  aucinf.obs aumcinf.obs
#>      0 Inf 62.4 [1.85]  468 [20.2]
#> 
#> Caption: aucinf.obs, aumcinf.obs: estimate and standard error

The sparse_lambda_z_se option (set globally with PKNCA.options() or for one analysis with the options argument of PKNCAdata()) is "delta" (the default) for the delta method or "none" for Yuan’s method. The method used is recorded in the PPANMETH column.

After an IV bolus, request the IV parameters; \(C_0\) is back-extrapolated from the mean profile, and no imputation is needed.

d_sparse_iv <-
  data.frame(
    time = rep(c(0.25, 0.5, 1, 2, 4, 8, 24), each = 4),
    conc =
      c(9.1, 10.4, 8.7, 9.9,  8.8, 9.5, 8.1, 9.2,  8.2, 7.6, 8.9, 7.9,
        6.4, 7.1, 6.0, 6.9,  4.6, 4.1, 5.0, 4.3,  2.1, 1.8, 2.3, 1.9,
        0.10, 0.08, 0.12, 0.09)
  )
d_sparse_iv$animal <- seq_len(nrow(d_sparse_iv))
o_data_iv <-
  PKNCAdata(
    PKNCAconc(d_sparse_iv, conc ~ time | animal, sparse = TRUE),
    PKNCAdose(data.frame(time = 0, dose = 100), dose ~ time, route = "intravascular"),
    intervals = data.frame(start = 0, end = Inf, c0 = TRUE, aucivinf.obs = TRUE)
  )
d_result_iv <- as.data.frame(suppressWarnings(pk.nca(o_data_iv)))
#> The sparse estimators use the linear trapezoidal rule, so the auc.method option
#> ("lin up/log down") does not apply to: aucinf.obs, aucivinf.obs
d_result_iv[grepl("^c0$|^aucivinf", d_result_iv$PPTESTCD), c("PPTESTCD", "PPORRES")]
#> # A tibble: 4 × 2
#>   PPTESTCD        PPORRES
#>   <chr>             <dbl>
#> 1 c0                10.2 
#> 2 aucivinf.obs      58.0 
#> 3 aucivinf.obs_se    1.35
#> 4 aucivinf.obs_df    5.97

The warnings suppressed here say that the AUCs from time 0 without \(C_0\) (aucinf.obs, which aucivinf.obs depends on) are not calculated.

Discussion

When \(\lambda_z\) is estimated from the same samples, its uncertainty adds to that of the AUC to infinity, and treating it as known, as Yuan’s (1993) standard error does, under-covers once the extrapolated area is a substantial part of the AUC. The delta method on the log means is a consistent extension of Yuan’s method: it keeps his estimate, applies to serial-sacrifice and batch designs alike through the covariance of the means, includes the covariance of \(\hat{\lambda}_z\) with \(C_{last}\) and the AUC to \(t_{last}\), and needs nothing beyond the quantities already used for the sparse AUC to \(t_{last}\). In the simulation it brought the coverage of the AUC to infinity much closer to nominal when the extrapolated area was large (90.2 to 93.6%, against 81.3 to 89.1% for Yuan), though still a few points short, and it cost little when the extrapolated area was small. For this reason, it is PKNCA’s default (sparse_lambda_z_se = "delta"). The same delta method carries the uncertainty of the back-extrapolated \(C_0\) after an IV bolus, and the IV bolus scenarios covered as well as the oral ones.

Estimating \(\lambda_z\) from the individual samples with a (mixed-effects) Tobit regression gave a more stable point estimate when samples were BLQ, but adding its variance to Yuan’s did not give adequate coverage, and its mixed-effects fit was itself occasionally unstable in the batch design with BLQ samples and slow elimination. A bootstrap over animals, which includes the covariance of that \(\lambda_z\) with the means, was too narrow with few animals per time, as every bootstrap was. The Tobit fit of the mean profile, which uses the LLOQs, gave nearly the same results as the log-linear fit. Neither Tobit fit removed the bias of the sparse estimate with BLQ samples, which comes mostly from counting BLQ samples as zero in the late means and from the mean profile ending earlier; estimating those means with the censoring taken into account would be the place to reduce it.

Limitations remain. The delta method conditions on the selected half-life points, so the uncertainty of the automatic point selection is not included; it is a first-order approximation, which is least accurate when \(\lambda_z\) is imprecise, as for the AUMC with slow elimination; and the arithmetic mean profile of animals with different elimination rates is not monoexponential, so the extrapolation \(C_{last}/\lambda_z\) is itself biased for the population mean curve, which is part of why the AUMC to infinity is underestimated. Where the extrapolated area is large, sampling longer remains the better remedy.

References

Bailer AJ. Testing for the equality of area under the curves when using destructive measurement techniques. Journal of Pharmacokinetics and Biopharmaceutics. 1988;16(3):303-309.

Holder DJ. Comments on Nedelman and Jia’s extension of Satterthwaite’s approximation applied to pharmacokinetics. Journal of Biopharmaceutical Statistics. 2001;11(1-2):75-79. doi:10.1081/BIP-100104199

Nedelman JR, Gibiansky E, Lau DTW. Applying Bailer’s method for AUC confidence intervals to sparse sampling. Pharmaceutical Research. 1995;12(1):124-128.

Nedelman JR, Jia X. An extension of Satterthwaite’s approximation applied to pharmacokinetics. Journal of Biopharmaceutical Statistics. 1998;8(2):317-328. doi:10.1080/10543409808835241

Shen M, Machado SG. Bioequivalence evaluation of sparse sampling pharmacokinetics data using bootstrap resampling method. Journal of Biopharmaceutical Statistics. 2017;27(2):257-264. doi:10.1080/10543406.2016.1265543

Takemoto S, Yamaoka K, Nishikawa M, Takakura Y. Histogram analysis of pharmacokinetic parameters by bootstrap resampling from one-point sampling data in animal experiments. Drug Metabolism and Pharmacokinetics. 2006;21(6):458-464. doi:10.2133/dmpk.21.458

Yuan J. Estimation of variance for AUC in animal studies. Journal of Pharmaceutical Sciences. 1993;82(7):761-763. doi:10.1002/jps.2600820718