6 Half-Life and Terminal Phase

6.1 How PKNCA estimates λz

PKNCA uses curve stripping: it systematically evaluates every terminal window ending at the last measurable point (last 3 points, last 4, and so on up to the full tail), fits log-linear regression to each subset, and selects the best fit based on adjusted R².

Before any window is fit, BLQ and zero concentrations are dropped, and so are points at or before the end of the last dose (including an infusion’s duration). Only points after Tmax are candidates unless allow.tmax.in.half.life = TRUE.

The selection criterion uses a tolerance window: among all subsets that yield a valid fit (λz > 0), PKNCA picks the one with the most points whose adjusted R² is within adj.r.squared.factor of the best adjusted R²:

selected n = the largest n whose adj.R² is within adj.r.squared.factor of the best adj.R²

The small tolerance (adj.r.squared.factor, default 0.0001) means: prefer longer windows (more points → more accurate λz) unless the longer window’s adj.R² is noticeably worse than the best. This is a model-selection tolerance, not a per-point penalty added to R². Selection can use the unadjusted R² instead, and the candidate windows can be limited by number of points and start time; see Controlling automatic point selection.

flowchart TD
    A["Drop BLQ/zero concentrations and<br/>points at or before the end of the last dose"] --> B["Candidates: points after Tmax<br/>(and at or after min.hl.start.time)"]
    B --> C["Fit log-linear regression to the last 3, last 4, ...<br/>points (min.hl.points to max.hl.points)"]
    C --> D{"Valid fit?<br/>λz > 0"}
    D -- yes --> E["Compute adj.R² (or R²)"]
    D -- no --> F["Discard subset"]
    E --> G["Select the longest window within<br/>adj.r.squared.factor (or r.squared.factor) of the best"]
    G --> H["Return λz, half.life, R², n.points"]

The algorithm is described step-by-step in the Half-Life Calculation vignette.

6.2 Basic half-life output

hl_interval <- data.frame(
  start             = 0,
  end               = Inf,
  half.life         = TRUE,
  lambda.z          = TRUE,
  lambda.z.n.points = TRUE,
  lambda.z.time.first = TRUE,
  r.squared         = TRUE,
  adj.r.squared     = TRUE,
  span.ratio        = TRUE
)

o_nca <- pk.nca(PKNCAdata(o_conc, o_dose, intervals = hl_interval))

as.data.frame(o_nca) |>
  filter(PPTESTCD %in% c("half.life", "lambda.z", "lambda.z.n.points",
                          "lambda.z.time.first", "adj.r.squared", "span.ratio")) |>
  select(subject, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  arrange(subject)

# A tibble: 12 × 7
   subject lambda.z adj.r.squared lambda.z.time.first lambda.z.n.points
   <ord>      <dbl>         <dbl>               <dbl>             <dbl>
 1 6         0.0878         0.998                2.03                 7
 2 7         0.0883         0.998                6.98                 4
 3 8         0.0815         0.989                3.53                 6
 4 11        0.0955         1.000                9.03                 3
 5 3         0.102          0.999                9                    3
 6 2         0.104          0.996                7.03                 4
 7 4         0.0993         0.998                9.02                 3
 8 9         0.0825         0.999                8.8                  3
 9 12        0.110          0.999                9.03                 3
10 10        0.0750         0.999                9.38                 3
11 1         0.0485         1.000                9.05                 3
12 5         0.0866         0.998                7.02                 4
# ℹ 2 more variables: half.life <dbl>, span.ratio <dbl>

What each output means:

Parameter	Meaning
`lambda.z`	Terminal elimination rate constant (h⁻¹)
`half.life`	ln(2) / λz
`lambda.z.n.points`	Number of points used in the selected regression
`lambda.z.time.first`	First timepoint included in the regression window
`r.squared` / `adj.r.squared`	Fit quality of the selected regression
`span.ratio`	(time of last − time of first regression point) / half.life — post-hoc quality flag

6.3 Quality filters

PKNCA provides limits enforced during automatic point selection, several post-hoc exclusion thresholds applied after fitting with exclude(), and one fitting option:

Important distinction: min.hl.points, max.hl.points, and min.hl.start.time are enforced inside pk.calc.half.life() while it selects points — with too few candidate points, half.life is NA immediately and the exclude column gives the reason. min.hl.r.squared and min.span.ratio are not consulted during fitting; they are thresholds for exclude() calls after pk.nca().

6.3.1 Minimum number of points (enforced during fitting)

# Default: 3 points minimum; with too few points, half.life is NA during pk.nca()
o_data_4pts <- PKNCAdata(o_conc, o_dose, intervals = hl_interval,
                          options = list(min.hl.points = 4))
o_nca_4pts  <- pk.nca(o_data_4pts)

as.data.frame(o_nca_4pts) |>
  filter(PPTESTCD %in% c("half.life", "lambda.z.n.points")) |>
  select(subject, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  arrange(subject)

# A tibble: 12 × 3
   subject lambda.z.n.points half.life
   <ord>               <dbl>     <dbl>
 1 6                       7      7.89
 2 7                       4      7.85
 3 8                       6      8.51
 4 11                      4      7.22
 5 3                       6      7.36
 6 2                       4      6.66
 7 4                       4      7.32
 8 9                       4      8.70
 9 12                      5      6.67
10 10                      4      9.46
11 1                       5     14.4 
12 5                       4      8.00

6.3.2 R² threshold (post-hoc exclusion via `exclude()`)

min.hl.r.squared sets the minimum unadjusted R² threshold. It has no effect during pk.nca(); use it with exclude_nca_min.hl.r.squared() after fitting. Theoph’s log-linear fits are all excellent — every subject’s R² is above 0.99 — so a conventional threshold like 0.95 flags nothing here; this demonstration uses a stricter 0.995 to make the effect visible:

# Compute half-life for all subjects
o_nca_r2 <- pk.nca(PKNCAdata(o_conc, o_dose, intervals = hl_interval))

# Then flag results with unadjusted R² < 0.995
o_nca_r2_excl <- exclude(
  o_nca_r2,
  reason = "r.squared < 0.995",
  FUN    = exclude_nca_min.hl.r.squared(0.995)
)

as.data.frame(o_nca_r2_excl) |>
  filter(PPTESTCD %in% c("half.life", "r.squared")) |>
  select(subject, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  arrange(subject)

# A tibble: 12 × 3
   subject r.squared half.life
   <ord>       <dbl>     <dbl>
 1 6           0.998      7.89
 2 7           0.999      7.85
 3 8           0.991      8.51
 4 11          1.000      7.26
 5 3           0.999      6.77
 6 2           0.997      6.66
 7 4           0.999      6.98
 8 9           0.999      8.41
 9 12          0.999      6.29
10 10          1.000      9.25
11 1           1.000     14.3 
12 5           0.999      8.00

The flag itself is recorded in the exclude column. Subject 8 (R² = 0.991) is the only one below the threshold:

as.data.frame(o_nca_r2_excl) |>
  filter(!is.na(exclude)) |>
  distinct(subject, exclude)

# A tibble: 1 × 2
  subject exclude          
  <ord>   <chr>            
1 8       r.squared < 0.995

For the stricter adjusted R² filter: exclude_nca_min.hl.adj.r.squared(min).

6.3.3 Span ratio threshold (post-hoc exclusion via `exclude()`)

min.span.ratio has no effect during pk.nca(); use it with exclude_nca_span.ratio():

# Flag half-lives with span.ratio < 2
o_nca_span_excl <- exclude(
  o_nca_r2,
  reason = "span.ratio < 2",
  FUN    = exclude_nca_span.ratio(2)
)

as.data.frame(o_nca_span_excl) |>
  filter(PPTESTCD %in% c("half.life", "span.ratio")) |>
  select(subject, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  arrange(subject)

# A tibble: 12 × 3
   subject half.life span.ratio
   <ord>       <dbl>      <dbl>
 1 6            7.89       2.76
 2 7            7.85       2.20
 3 8            8.51       2.42
 4 11           7.26       2.07
 5 3            6.77       2.24
 6 2            6.66       2.59
 7 4            6.98       2.24
 8 9            8.41       1.86
 9 12           6.29       2.41
10 10           9.25       1.55
11 1           14.3        1.07
12 5            8.00       2.17

Three subjects fall below a span ratio of 2 and are flagged in the exclude column:

as.data.frame(o_nca_span_excl) |>
  filter(!is.na(exclude)) |>
  distinct(subject, exclude)

# A tibble: 3 × 2
  subject exclude       
  <ord>   <chr>         
1 9       span.ratio < 2
2 10      span.ratio < 2
3 1       span.ratio < 2

6.3.4 Allowing or excluding Tmax

By default Tmax is excluded from the λz regression (drug may still be distributing).

o_data_tmax <- PKNCAdata(o_conc, o_dose, intervals = hl_interval,
                          options = list(allow.tmax.in.half.life = TRUE))
o_nca_tmax  <- pk.nca(o_data_tmax)

as.data.frame(o_nca_tmax) |>
  filter(PPTESTCD %in% c("half.life", "lambda.z.n.points", "lambda.z.time.first")) |>
  select(subject, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  arrange(subject)

# A tibble: 12 × 4
   subject lambda.z.time.first lambda.z.n.points half.life
   <ord>                 <dbl>             <dbl>     <dbl>
 1 6                      2.03                 7      7.89
 2 7                      6.98                 4      7.85
 3 8                      2.02                 7      8.47
 4 11                     9.03                 3      7.26
 5 3                      9                    3      6.77
 6 2                      7.03                 4      6.66
 7 4                      9.02                 3      6.98
 8 9                      8.8                  3      8.41
 9 12                     9.03                 3      6.29
10 10                     9.38                 3      9.25
11 1                      9.05                 3     14.3 
12 5                      7.02                 4      8.00

6.3.5 Controlling automatic point selection

Three options (≥ 0.13.0) change which windows curve stripping considers or how it ranks them. Each is a PKNCA.options() option (set per analysis with PKNCAdata(options = ...)) and an argument of pk.calc.half.life() with the same name:

Option	Default	Effect
`max.hl.points`	`Inf` (no limit)	The most points a window may use; must be more than `min.hl.points`
`min.hl.start.time`	`0` (no restriction)	The earliest time a window may start; within `pk.nca()` it is the time since the start of the interval
`r.squared.factor`	`NA` (not in use)	Select on the unadjusted R² with this tolerance instead of on the adjusted R²

max.hl.points and min.hl.start.time mirror the “Lambda Z Rules for Best Fit” of Phoenix WinNonlin, and they apply to both the log-linear and the Tobit methods. adj.r.squared.factor and r.squared.factor are alternatives: setting either one takes the other out of use, so exactly one is in use at a time, and setting one to NA selects on the other. Unlike the adjusted R², the unadjusted R² does not reward more points, so it usually selects fewer of them.

hl_select_interval <- data.frame(
  start               = 0,
  end                 = Inf,
  half.life           = TRUE,
  lambda.z.n.points   = TRUE,
  lambda.z.time.first = TRUE
)

selection_settings <- list(
  "default"               = list(),
  "max.hl.points=4"       = list(max.hl.points = 4),
  "min.hl.start.time=9"   = list(min.hl.start.time = 9),
  "r.squared.factor=1e-4" = list(r.squared.factor = 0.0001)
)

hl_select <- lapply(names(selection_settings), function(setting) {
  o_data_setting <- PKNCAdata(o_conc, o_dose, intervals = hl_select_interval,
                              options = selection_settings[[setting]])
  as.data.frame(pk.nca(o_data_setting)) |>
    mutate(setting = setting)
}) |>
  bind_rows()

Warning: Too few points for half-life calculation (min.hl.points=3 with only 2
points at or after min.hl.start.time=9)

hl_select |>
  filter(subject %in% c("6", "8", "9", "11"),
         PPTESTCD %in% c("lambda.z.time.first", "lambda.z.n.points", "half.life")) |>
  select(subject, setting, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  arrange(subject, setting)

# A tibble: 16 × 5
   subject setting               lambda.z.time.first lambda.z.n.points half.life
   <ord>   <chr>                               <dbl>             <dbl>     <dbl>
 1 6       default                              2.03                 7      7.89
 2 6       max.hl.points=4                      9.22                 3      7.57
 3 6       min.hl.start.time=9                  9.22                 3      7.57
 4 6       r.squared.factor=1e-4                9.22                 3      7.57
 5 8       default                              3.53                 6      8.51
 6 8       max.hl.points=4                      7.15                 4      8.59
 7 8       min.hl.start.time=9                  9.07                 3      8.42
 8 8       r.squared.factor=1e-4                3.53                 6      8.51
 9 11      default                              9.03                 3      7.26
10 11      max.hl.points=4                      9.03                 3      7.26
11 11      min.hl.start.time=9                  9.03                 3      7.26
12 11      r.squared.factor=1e-4                7.03                 4      7.22
13 9       default                              8.8                  3      8.41
14 9       max.hl.points=4                      8.8                  3      8.41
15 9       min.hl.start.time=9                 NA                   NA     NA   
16 9       r.squared.factor=1e-4                8.8                  3      8.41

For the four subjects shown:

max.hl.points = 4 shortens subject 6’s window from 7 points to 3 (the better of the 3- and 4-point windows) and subject 8’s from 6 to 4; subjects 9 and 11 already used 3 points.
min.hl.start.time = 9 moves the start of subject 6’s and subject 8’s windows to their first samples after 9 h. Subject 9 has only two measurable points at or after 9 h, so its half-life is not calculated; that is the warning above, and the same text is recorded in the exclude column:

hl_select |>
  filter(setting == "min.hl.start.time=9", subject == "9",
         PPTESTCD == "half.life") |>
  pull(exclude)

[1] "Too few points for half-life calculation (min.hl.points=3 with only 2 points at or after min.hl.start.time=9)"

r.squared.factor = 1e-4 selects a shorter window for subject 6 (3 points instead of 7) but a longer one for subject 11 (4 instead of 3), and leaves subjects 8 and 9 unchanged.

When no window survives selection — for example, when a well-fitting window with λz ≤ 0 sets the R² tolerance so that no window with λz > 0 is within it — half.life is NA with an exclusion reason beginning “No valid terminal phase” (≥ 0.13.0).

6.4 Manually controlling which points are used

See also manual point selection in the Half-Life Calculation vignette.

6.4.1 Excluding specific points (`exclude_half.life`)

Both point-selection columns are optional columns of the concentration data, named with the exclude_half.life and include_half.life arguments of PKNCAconc(). Since PKNCA 0.13.0 they must be logical (TRUE/FALSE/NA), naming a column that does not exist is an error, and only one of the two may be in use for an interval (using both is an error). A column is “in use” for an interval when any of its values there is not NA — even if every value is FALSE — so initialize the column to NA wherever it should not apply. The Half-Life Calculation vignette has a decision tree for choosing among automatic selection, exclude_half.life, and include_half.life.

Mark individual observations that should never be included in the λz regression — for example, a suspected outlier or a secondary peak. The remaining points go through the usual curve stripping.

d_conc_excl <- d_conc |>
  mutate(excl_hl = ifelse(subject == "1" & time > 16, TRUE, NA))

o_conc_excl <- PKNCAconc(d_conc_excl, conc ~ time | subject,
                          exclude_half.life = "excl_hl")

o_nca_excl <- pk.nca(PKNCAdata(o_conc_excl, o_dose, intervals = hl_interval))

# Compare half.life and n.points for subject 1
bind_rows(
  as.data.frame(o_nca) |>
    filter(subject == "1", PPTESTCD %in% c("half.life", "lambda.z.n.points")) |>
    mutate(version = "default"),
  as.data.frame(o_nca_excl) |>
    filter(subject == "1", PPTESTCD %in% c("half.life", "lambda.z.n.points")) |>
    mutate(version = "late points excluded")
) |>
  select(version, PPTESTCD, PPORRES) |>
  arrange(version, PPTESTCD)

# A tibble: 4 × 3
  version              PPTESTCD          PPORRES
  <chr>                <chr>               <dbl>
1 default              half.life            14.3
2 default              lambda.z.n.points     3  
3 late points excluded half.life            15.3
4 late points excluded lambda.z.n.points     3

6.4.2 Specifying the exact points (`include_half.life`)

Mark the exact points to fit: include_half.life bypasses curve stripping entirely and fits only the marked points. Automatic selection — including min.hl.points, max.hl.points, min.hl.start.time, adj.r.squared.factor, r.squared.factor, and allow.tmax.in.half.life — is not applied, so even a two-point fit is accepted. Useful when the regression window is pre-specified, e.g. by an analysis plan.

d_conc_incl <- d_conc |>
  mutate(incl_hl = ifelse(subject == "1" & time %in% c(9.05, 12.12), TRUE, NA))

o_conc_incl <- PKNCAconc(d_conc_incl, conc ~ time | subject,
                          include_half.life = "incl_hl")

o_nca_incl <- pk.nca(PKNCAdata(o_conc_incl, o_dose, intervals = hl_interval))

Warning: n must be > 2 for adj.r.squared

as.data.frame(o_nca_incl) |>
  filter(subject == "1",
         PPTESTCD %in% c("half.life", "lambda.z.n.points", "lambda.z.time.first")) |>
  select(PPTESTCD, PPORRES)

# A tibble: 3 × 2
  PPTESTCD            PPORRES
  <chr>                 <dbl>
1 lambda.z.time.first    9.05
2 lambda.z.n.points      2   
3 half.life             14.3

The two-point fit is accepted, with the warning “n must be > 2 for adj.r.squared”: the adjusted R² is undefined for two points, so it is NA for this fit while half.life is reported.

Note: naming a point in include_half.life does not force it into the fit. BLQ and zero concentrations, and points at or before the end of the last dose (including an infusion’s duration), are dropped first, exactly as for automatic selection.

6.5 Tobit regression for BLQ-heavy data (≥ 0.13.0)

6.5.1 Why Tobit?

In studies where drug concentrations in the terminal elimination phase are frequently below the lower limit of quantification (BLQ), ordinary least-squares (OLS) log-linear regression produces biased estimates of λz. The bias arises because OLS either:

drops BLQ observations entirely — shrinking the regression window and losing information about the true terminal slope, or
substitutes BLQ = 0 — which is undefined on the log scale and forces arbitrary imputation.

Tobit regression handles this correctly by treating each BLQ observation as a left-censored data point. Rather than requiring a specific concentration value, it only asserts C ≤ LLOQ. The regression uses the full profile, including censored observations, without fabricating values.

6.5.2 The key statistical difference from OLS

OLS minimizes the sum of squared residuals across observed concentrations. Tobit instead maximizes a likelihood that has two components:

Observed concentrations (above LLOQ): contribute via the normal log-likelihood of the log-transformed concentration given the fitted log-linear line.
BLQ observations: contribute the cumulative probability that the true (unobserved) concentration falls at or below the LLOQ, i.e. Φ((log(LLOQ) − predicted) / σ), where Φ is the standard normal CDF and σ the residual standard deviation.

This means BLQ timepoints actively inform the slope estimate — they pull the λz estimate toward the true terminal slope — rather than being discarded.

6.5.3 When to consider Tobit

Tobit regression is most beneficial in the following situations:

A substantial fraction of terminal-phase observations — as a rough rule of thumb (not a documented PKNCA threshold), a fifth or more — are BLQ. With few or no BLQ points, OLS and Tobit give similar estimates.
The drug has low tissue concentrations or a short distribution half-life relative to the sampling schedule, making late timepoints frequently BLQ.
Pediatric or geriatric studies with limited sampling windows, where sparse late-time samples are often at or near the LLOQ.
The terminal slope is important for downstream parameters (AUC extrapolation, mean residence time) and precision matters for regulatory submissions.

6.5.4 Using Tobit directly with `pk.calc.half.life()` (≥ 0.13.0)

Tobit can be called directly through pk.calc.half.life() or used in the full pk.nca() pipeline (next section). For a direct call, supply the full concentration-time vector and the LLOQ (a scalar, or a vector with one value per timepoint):

# Simulated profile — concentrations at 16h and 24h are below the LLOQ of 0.1 mg/L
conc_tobit <- c(0, 2.5, 4.8, 4.2, 2.9, 1.4, 0.6, 0.05, 0.01)
time_tobit <- c(0, 0.5,   1,   2,   4,   8,  12,  16,  24)

# OLS reference — the BLQ observations (16 and 24 h) are dropped
is_blq <- conc_tobit > 0 & conc_tobit < 0.1
hl_ols   <- pk.calc.half.life(conc = conc_tobit[!is_blq], time = time_tobit[!is_blq],
                               tmax = 1, tlast = 12, hl_method = "log-linear")

# Tobit — includes the BLQ observations as left-censored at the scalar lloq
hl_tobit <- pk.calc.half.life(conc = conc_tobit, time = time_tobit, tmax = 1,
                               tlast = 24, lloq = 0.1, hl_method = "tobit")

# Compare
data.frame(
  method    = c("log-linear", "tobit"),
  lambda.z  = c(hl_ols$lambda.z,  hl_tobit$lambda.z),
  half.life = c(hl_ols$half.life, hl_tobit$half.life),
  n.points  = c(hl_ols$lambda.z.n.points, hl_tobit$lambda.z.n.points)
)

      method  lambda.z half.life n.points
1 log-linear 0.1936354  3.579651        4
2      tobit 0.2658595  2.607194        6

With the two BLQ samples dropped, the log-linear fit uses the 4 remaining points after Tmax and gives a half-life of about 3.58 h. Tobit keeps the 16- and 24-hour samples as censored at the LLOQ, so its lambda.z.n.points of 6 counts them, and the steeper slope they imply gives a half-life of about 2.61 h.

6.5.5 Tobit through the `pk.nca()` pipeline (≥ 0.13.0)

To run Tobit across a full dataset, declare the LLOQ on the concentration object with the lloq argument of PKNCAconc() (a column name or a scalar) and select the method with the hl_method option. Every non-missing lloq value must be a finite number greater than zero; anything else is an error (pknca_error_conc_invalid_lloq). Missing values are allowed, and an entirely NA column is accepted.

# Subject 1 with an assay LLOQ of 0.5 mg/L
d_conc_subset <- d_conc |> filter(subject == "1")
d_dose_subset <- d_dose |> filter(subject == "1")

o_conc_tobit <- PKNCAconc(d_conc_subset, conc ~ time | subject, lloq = 0.5)
o_dose_tobit <- PKNCAdose(d_dose_subset, dose ~ time | subject, route = "extravascular")

o_data_tobit <- PKNCAdata(
  o_conc_tobit, o_dose_tobit,
  intervals = data.frame(start = 0, end = Inf, half.life = TRUE),
  options   = list(hl_method = "tobit")
)
o_nca_tobit <- pk.nca(o_data_tobit)

as.data.frame(o_nca_tobit) |>
  filter(PPTESTCD %in% c("half.life", "lambda.z.n.points_blq",
                          "tobit_residual", "adj_tobit_residual")) |>
  select(PPTESTCD, PPORRES)

# A tibble: 4 × 2
  PPTESTCD                 PPORRES
  <chr>                      <dbl>
1 lambda.z.n.points_blq  0        
2 half.life             14.3      
3 tobit_residual         0.000167 
4 adj_tobit_residual     0.0000834

None of subject 1’s terminal-phase samples fall below this LLOQ (lambda.z.n.points_blq is 0), so the Tobit estimate matches the OLS fit here — the example shows the wiring, and the Tobit-specific outputs appear alongside the usual half-life results.

Tobit adds new output parameters to the results:

Parameter	Description
`lambda.z.n.points_blq`	Number of BLQ points included in the Tobit regression
`tobit_residual`	Estimated residual standard deviation of the Tobit fit (log-concentration scale)
`adj_tobit_residual`	Adjusted Tobit residual (analogous to adjusted R²)

6.5.6 Tobit-related options

Option	Default	Effect
`hl_method`	`"log-linear"`	`"log-linear"` (OLS) or `"tobit"` (censored-likelihood regression)
`tobit_n_points_penalty`	`0`	Exponent in the window-selection criterion (residual × n^penalty); positive values favor shorter windows
`tobit_optim_control`	`list()`	Control list forwarded to `stats::optim()` for the Tobit optimization

Window selection also differs from the log-linear fits: instead of the adjusted-R² rule, PKNCA fits each candidate terminal window and keeps the one that minimizes tobit_residual × n_points^tobit_n_points_penalty, with ties going to the larger window. min.hl.points and max.hl.points count only the above-LLOQ points of a Tobit window, and min.hl.start.time applies as it does for the log-linear method.

6.5.7 Tobit with sparse PK (≥ 0.13.0)

The Tobit half-life also works with sparse sampling. For sparse data the half-life is calculated on the arithmetic-mean profile, and that profile carries an LLOQ at each time: the median of the lloq values of its samples, ignoring missing ones. A time whose mean is set to zero because more than half of its samples are BLQ is therefore censored at that LLOQ rather than dropped. See the sparse NCA chapter for sparse analyses.

See the vignette Half-life calculation with Tobit regression for a detailed simulation comparison between OLS and Tobit estimates.

6.6 `lambda.z.corrxy` — regression correlation

lambda.z.corrxy (≥ 0.13.0) returns the Pearson correlation between time and log-concentration for the points used in the λz regression (numeric NA when it cannot be computed). It is a unitless diagnostic — values close to −1 indicate a strong, clean terminal slope. Like r.squared and adj.r.squared, it is returned only by the log-linear method.

hl_corrxy_interval <- data.frame(
  start              = 0,
  end                = Inf,
  half.life          = TRUE,
  lambda.z           = TRUE,
  lambda.z.corrxy    = TRUE,
  r.squared          = TRUE,
  adj.r.squared      = TRUE
)

o_nca_corrxy <- pk.nca(PKNCAdata(o_conc, o_dose, intervals = hl_corrxy_interval))

as.data.frame(o_nca_corrxy) |>
  filter(PPTESTCD %in% c("half.life", "lambda.z.corrxy", "r.squared", "adj.r.squared")) |>
  select(subject, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |>
  arrange(subject)

# A tibble: 12 × 5
   subject r.squared adj.r.squared lambda.z.corrxy half.life
   <ord>       <dbl>         <dbl>           <dbl>     <dbl>
 1 6           0.998         0.998          -0.999      7.89
 2 7           0.999         0.998          -0.999      7.85
 3 8           0.991         0.989          -0.995      8.51
 4 11          1.000         1.000          -1.000      7.26
 5 3           0.999         0.999          -1.000      6.77
 6 2           0.997         0.996          -0.999      6.66
 7 4           0.999         0.998          -0.999      6.98
 8 9           0.999         0.999          -1.000      8.41
 9 12          0.999         0.999          -1.000      6.29
10 10          1.000         0.999          -1.000      9.25
11 1           1.000         1.000          -1.000     14.3 
12 5           0.999         0.998          -0.999      8.00

6.7 Visualizing the chosen regression window

Three functions extract what is needed to plot the fit, each from a PKNCAresults object or (by running only the half-life calculation) a PKNCAdata object:

get_halflife_points() gives one logical value per row of the concentration data: TRUE for points used in the fit, FALSE for points not used, and NA where no half-life was calculated for the interval.
get_halflife_fit() (≥ 0.13.0) gives the intercept and slope of the line log(conc) = intercept + slope*time and the time_first/time_last of the points used, one row per group and interval. Times are on the scale of the concentration data, not relative to the interval start, and the line is NA when the half-life was not calculated or was excluded.
get_halflife_curve() (≥ 0.13.0) gives concentrations along that line, at the times in tout or at n equally spaced times (default 50) spanning the points used. It extrapolates after the last point used by default (extrapolate_later = TRUE) and not before the first (extrapolate_earlier = FALSE); times outside the fit that are not extrapolated give NA.

get_halflife_fit(o_nca) |>
  filter(subject %in% c("1", "6"))

  subject start end intercept       slope time_first time_last
1       6     0 Inf  2.033404 -0.08779574       2.03     23.85
2       1     0 Inf  2.368785 -0.04845700       9.05     24.37

# Mark the points used for each subject's fit
d_conc_hl <- d_conc |>
  mutate(used_in_lambda_z = get_halflife_points(o_nca))

# The fitted line on a grid from 0 to 36 h: NA before the first point used,
# extrapolated after the last one
d_curve <- get_halflife_curve(o_nca, tout = seq(0, 36, by = 0.5),
                              extrapolate_earlier = FALSE)

ggplot(d_conc_hl |> filter(subject %in% c("1", "6"), conc > 0),
       aes(x = time, y = conc)) +
  geom_line(colour = "grey60") +
  geom_point(aes(colour = used_in_lambda_z), size = 3) +
  geom_line(data = d_curve |> filter(subject %in% c("1", "6"), !is.na(conc)),
            colour = "firebrick", linetype = "dashed") +
  scale_colour_manual(values = c("FALSE" = "grey60", "TRUE" = "firebrick"),
                      labels = c("FALSE" = "not used in λz", "TRUE" = "used in λz"),
                      name = NULL) +
  scale_y_log10() +
  facet_wrap(~subject, labeller = label_both) +
  labs(title = "Terminal regression window and λz fit (log scale)",
       x = "Time (h)", y = "Concentration (mg/L) [log scale]") +
  theme_minimal()

Subject 1’s fit uses its last 3 points (from 9.05 h), while subject 6’s uses 7 points (from 2.03 h). The dashed line starts at the first point used because extrapolate_earlier = FALSE, and it continues to 36 h because later times are extrapolated by default. For many subjects at once, the same functions feed a faceted plot or a geom_line() overlay in a QC report.

pkgdown reference: PKNCAconc() · PKNCAdose() · PKNCAdata() · pk.nca() · pk.calc.half.life() · PKNCA.options() · get_halflife_points() · get_halflife_fit() · get_halflife_curve() · exclude() · exclude_nca_min.hl.r.squared() · exclude_nca_span.ratio() · exclude_nca_min.hl.adj.r.squared()

--- title: "Half-Life and Terminal Phase" --- ```{r setup, include=FALSE} library(PKNCA) library(dplyr) library(ggplot2) conflicted::conflicts_prefer(dplyr::filter, dplyr::select, .quiet = TRUE) 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], .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") ``` ## How PKNCA estimates λz PKNCA uses **curve stripping**: it systematically evaluates every terminal window ending at the last measurable point (last 3 points, last 4, and so on up to the full tail), fits log-linear regression to each subset, and selects the best fit based on adjusted R². Before any window is fit, BLQ and zero concentrations are dropped, and so are points at or before the end of the last dose (including an infusion's duration). Only points after Tmax are candidates unless `allow.tmax.in.half.life = TRUE`. The selection criterion uses a **tolerance window**: among all subsets that yield a valid fit (λz > 0), PKNCA picks the one with the **most points** whose adjusted R² is within `adj.r.squared.factor` of the best adjusted R²: > selected n = the largest n whose adj.R² is within `adj.r.squared.factor` of the best adj.R² The small tolerance (`adj.r.squared.factor`, default 0.0001) means: prefer longer windows (more points → more accurate λz) unless the longer window's adj.R² is noticeably worse than the best. This is a **model-selection tolerance**, not a per-point penalty added to R². Selection can use the unadjusted R² instead, and the candidate windows can be limited by number of points and start time; see [Controlling automatic point selection](#controlling-automatic-point-selection). ```{mermaid} flowchart TD A["Drop BLQ/zero concentrations and points at or before the end of the last dose"] --> B["Candidates: points after Tmax (and at or after min.hl.start.time)"] B --> C["Fit log-linear regression to the last 3, last 4, ... points (min.hl.points to max.hl.points)"] C --> D{"Valid fit? λz > 0"} D -- yes --> E["Compute adj.R² (or R²)"] D -- no --> F["Discard subset"] E --> G["Select the longest window within adj.r.squared.factor (or r.squared.factor) of the best"] G --> H["Return λz, half.life, R², n.points"] ``` The algorithm is described step-by-step in the [Half-Life Calculation vignette](https://humanpred.github.io/pknca/articles/v06-half-life-calculation.html). --- ## Basic half-life output ```{r} hl_interval <- data.frame( start = 0, end = Inf, half.life = TRUE, lambda.z = TRUE, lambda.z.n.points = TRUE, lambda.z.time.first = TRUE, r.squared = TRUE, adj.r.squared = TRUE, span.ratio = TRUE ) o_nca <- pk.nca(PKNCAdata(o_conc, o_dose, intervals = hl_interval)) as.data.frame(o_nca) |> filter(PPTESTCD %in% c("half.life", "lambda.z", "lambda.z.n.points", "lambda.z.time.first", "adj.r.squared", "span.ratio")) |> select(subject, PPTESTCD, PPORRES) |> tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |> arrange(subject) ``` **What each output means:** | Parameter | Meaning | |---|---| | `lambda.z` | Terminal elimination rate constant (h⁻¹) | | `half.life` | ln(2) / λz | | `lambda.z.n.points` | Number of points used in the selected regression | | `lambda.z.time.first` | First timepoint included in the regression window | | `r.squared` / `adj.r.squared` | Fit quality of the selected regression | | `span.ratio` | (time of last − time of first regression point) / half.life — post-hoc quality flag | --- ## Quality filters PKNCA provides limits enforced **during automatic point selection**, several post-hoc exclusion thresholds applied **after fitting** with `exclude()`, and one fitting option: > **Important distinction:** `min.hl.points`, `max.hl.points`, and `min.hl.start.time` are enforced inside `pk.calc.half.life()` while it selects points — with too few candidate points, `half.life` is `NA` immediately and the `exclude` column gives the reason. `min.hl.r.squared` and `min.span.ratio` are **not** consulted during fitting; they are thresholds for `exclude()` calls after `pk.nca()`. ### Minimum number of points (enforced during fitting) ```{r} # Default: 3 points minimum; with too few points, half.life is NA during pk.nca() o_data_4pts <- PKNCAdata(o_conc, o_dose, intervals = hl_interval, options = list(min.hl.points = 4)) o_nca_4pts <- pk.nca(o_data_4pts) as.data.frame(o_nca_4pts) |> filter(PPTESTCD %in% c("half.life", "lambda.z.n.points")) |> select(subject, PPTESTCD, PPORRES) |> tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |> arrange(subject) ``` ### R² threshold (post-hoc exclusion via `exclude()`) `min.hl.r.squared` sets the minimum **unadjusted R²** threshold. It has no effect during `pk.nca()`; use it with `exclude_nca_min.hl.r.squared()` after fitting. Theoph's log-linear fits are all excellent — every subject's R² is above 0.99 — so a conventional threshold like 0.95 flags nothing here; this demonstration uses a stricter 0.995 to make the effect visible: ```{r} # Compute half-life for all subjects o_nca_r2 <- pk.nca(PKNCAdata(o_conc, o_dose, intervals = hl_interval)) # Then flag results with unadjusted R² < 0.995 o_nca_r2_excl <- exclude( o_nca_r2, reason = "r.squared < 0.995", FUN = exclude_nca_min.hl.r.squared(0.995) ) as.data.frame(o_nca_r2_excl) |> filter(PPTESTCD %in% c("half.life", "r.squared")) |> select(subject, PPTESTCD, PPORRES) |> tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |> arrange(subject) ``` The flag itself is recorded in the `exclude` column. Subject 8 (R² = 0.991) is the only one below the threshold: ```{r} as.data.frame(o_nca_r2_excl) |> filter(!is.na(exclude)) |> distinct(subject, exclude) ``` For the stricter **adjusted R²** filter: `exclude_nca_min.hl.adj.r.squared(min)`. ### Span ratio threshold (post-hoc exclusion via `exclude()`) `min.span.ratio` has no effect during `pk.nca()`; use it with `exclude_nca_span.ratio()`: ```{r} # Flag half-lives with span.ratio < 2 o_nca_span_excl <- exclude( o_nca_r2, reason = "span.ratio < 2", FUN = exclude_nca_span.ratio(2) ) as.data.frame(o_nca_span_excl) |> filter(PPTESTCD %in% c("half.life", "span.ratio")) |> select(subject, PPTESTCD, PPORRES) |> tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |> arrange(subject) ``` Three subjects fall below a span ratio of 2 and are flagged in the `exclude` column: ```{r} as.data.frame(o_nca_span_excl) |> filter(!is.na(exclude)) |> distinct(subject, exclude) ``` ### Allowing or excluding Tmax By default Tmax is excluded from the λz regression (drug may still be distributing). ```{r} o_data_tmax <- PKNCAdata(o_conc, o_dose, intervals = hl_interval, options = list(allow.tmax.in.half.life = TRUE)) o_nca_tmax <- pk.nca(o_data_tmax) as.data.frame(o_nca_tmax) |> filter(PPTESTCD %in% c("half.life", "lambda.z.n.points", "lambda.z.time.first")) |> select(subject, PPTESTCD, PPORRES) |> tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |> arrange(subject) ``` ### Controlling automatic point selection {#controlling-automatic-point-selection} Three options (≥ 0.13.0) change which windows curve stripping considers or how it ranks them. Each is a `PKNCA.options()` option (set per analysis with `PKNCAdata(options = ...)`) and an argument of `pk.calc.half.life()` with the same name: | Option | Default | Effect | |---|---|---| | `max.hl.points` | `Inf` (no limit) | The most points a window may use; must be more than `min.hl.points` | | `min.hl.start.time` | `0` (no restriction) | The earliest time a window may start; within `pk.nca()` it is the time since the start of the interval | | `r.squared.factor` | `NA` (not in use) | Select on the unadjusted R² with this tolerance instead of on the adjusted R² | `max.hl.points` and `min.hl.start.time` mirror the "Lambda Z Rules for Best Fit" of Phoenix WinNonlin, and they apply to both the log-linear and the Tobit methods. `adj.r.squared.factor` and `r.squared.factor` are alternatives: setting either one takes the other out of use, so exactly one is in use at a time, and setting one to `NA` selects on the other. Unlike the adjusted R², the unadjusted R² does not reward more points, so it usually selects fewer of them. ```{r} hl_select_interval <- data.frame( start = 0, end = Inf, half.life = TRUE, lambda.z.n.points = TRUE, lambda.z.time.first = TRUE ) selection_settings <- list( "default" = list(), "max.hl.points=4" = list(max.hl.points = 4), "min.hl.start.time=9" = list(min.hl.start.time = 9), "r.squared.factor=1e-4" = list(r.squared.factor = 0.0001) ) hl_select <- lapply(names(selection_settings), function(setting) { o_data_setting <- PKNCAdata(o_conc, o_dose, intervals = hl_select_interval, options = selection_settings[[setting]]) as.data.frame(pk.nca(o_data_setting)) |> mutate(setting = setting) }) |> bind_rows() hl_select |> filter(subject %in% c("6", "8", "9", "11"), PPTESTCD %in% c("lambda.z.time.first", "lambda.z.n.points", "half.life")) |> select(subject, setting, PPTESTCD, PPORRES) |> tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |> arrange(subject, setting) ``` For the four subjects shown: - **`max.hl.points = 4`** shortens subject 6's window from 7 points to 3 (the better of the 3- and 4-point windows) and subject 8's from 6 to 4; subjects 9 and 11 already used 3 points. - **`min.hl.start.time = 9`** moves the start of subject 6's and subject 8's windows to their first samples after 9 h. Subject 9 has only two measurable points at or after 9 h, so its half-life is not calculated; that is the warning above, and the same text is recorded in the `exclude` column: ```{r} hl_select |> filter(setting == "min.hl.start.time=9", subject == "9", PPTESTCD == "half.life") |> pull(exclude) ``` - **`r.squared.factor = 1e-4`** selects a shorter window for subject 6 (3 points instead of 7) but a longer one for subject 11 (4 instead of 3), and leaves subjects 8 and 9 unchanged. When no window survives selection — for example, when a well-fitting window with λz ≤ 0 sets the R² tolerance so that no window with λz > 0 is within it — `half.life` is `NA` with an exclusion reason beginning "No valid terminal phase" (≥ 0.13.0). --- ## Manually controlling which points are used See also [manual point selection](https://humanpred.github.io/pknca/articles/v06-half-life-calculation.html#manual-point-selection) in the Half-Life Calculation vignette. ### Excluding specific points (`exclude_half.life`) Both point-selection columns are optional columns of the concentration data, named with the `exclude_half.life` and `include_half.life` arguments of `PKNCAconc()`. Since PKNCA 0.13.0 they must be logical (`TRUE`/`FALSE`/`NA`), naming a column that does not exist is an error, and only one of the two may be in use for an interval (using both is an error). A column is "in use" for an interval when any of its values there is not `NA` — even if every value is `FALSE` — so initialize the column to `NA` wherever it should not apply. The [Half-Life Calculation vignette](https://humanpred.github.io/pknca/articles/v06-half-life-calculation.html#choosing-a-point-selection-method) has a decision tree for choosing among automatic selection, `exclude_half.life`, and `include_half.life`. Mark individual observations that should never be included in the λz regression — for example, a suspected outlier or a secondary peak. The remaining points go through the usual curve stripping. ```{r} d_conc_excl <- d_conc |> mutate(excl_hl = ifelse(subject == "1" & time > 16, TRUE, NA)) o_conc_excl <- PKNCAconc(d_conc_excl, conc ~ time | subject, exclude_half.life = "excl_hl") o_nca_excl <- pk.nca(PKNCAdata(o_conc_excl, o_dose, intervals = hl_interval)) # Compare half.life and n.points for subject 1 bind_rows( as.data.frame(o_nca) |> filter(subject == "1", PPTESTCD %in% c("half.life", "lambda.z.n.points")) |> mutate(version = "default"), as.data.frame(o_nca_excl) |> filter(subject == "1", PPTESTCD %in% c("half.life", "lambda.z.n.points")) |> mutate(version = "late points excluded") ) |> select(version, PPTESTCD, PPORRES) |> arrange(version, PPTESTCD) ``` ### Specifying the exact points (`include_half.life`) Mark the exact points to fit: `include_half.life` **bypasses curve stripping entirely** and fits only the marked points. Automatic selection — including `min.hl.points`, `max.hl.points`, `min.hl.start.time`, `adj.r.squared.factor`, `r.squared.factor`, and `allow.tmax.in.half.life` — is not applied, so even a two-point fit is accepted. Useful when the regression window is pre-specified, e.g. by an analysis plan. ```{r} d_conc_incl <- d_conc |> mutate(incl_hl = ifelse(subject == "1" & time %in% c(9.05, 12.12), TRUE, NA)) o_conc_incl <- PKNCAconc(d_conc_incl, conc ~ time | subject, include_half.life = "incl_hl") o_nca_incl <- pk.nca(PKNCAdata(o_conc_incl, o_dose, intervals = hl_interval)) as.data.frame(o_nca_incl) |> filter(subject == "1", PPTESTCD %in% c("half.life", "lambda.z.n.points", "lambda.z.time.first")) |> select(PPTESTCD, PPORRES) ``` The two-point fit is accepted, with the warning "n must be > 2 for adj.r.squared": the adjusted R² is undefined for two points, so it is `NA` for this fit while `half.life` is reported. > **Note:** naming a point in `include_half.life` does not force it into the fit. BLQ and zero concentrations, and points at or before the end of the last dose (including an infusion's duration), are dropped first, exactly as for automatic selection. --- ## Tobit regression for BLQ-heavy data (≥ 0.13.0) ### Why Tobit? In studies where drug concentrations in the terminal elimination phase are frequently below the lower limit of quantification (BLQ), ordinary least-squares (OLS) log-linear regression produces biased estimates of λz. The bias arises because OLS either: - **drops** BLQ observations entirely — shrinking the regression window and losing information about the true terminal slope, or - **substitutes** BLQ = 0 — which is undefined on the log scale and forces arbitrary imputation. Tobit regression handles this correctly by treating each BLQ observation as a **left-censored** data point. Rather than requiring a specific concentration value, it only asserts `C ≤ LLOQ`. The regression uses the full profile, including censored observations, without fabricating values. ### The key statistical difference from OLS OLS minimizes the sum of squared residuals across observed concentrations. Tobit instead **maximizes a likelihood** that has two components: 1. **Observed concentrations** (above LLOQ): contribute via the normal log-likelihood of the log-transformed concentration given the fitted log-linear line. 2. **BLQ observations**: contribute the cumulative probability that the true (unobserved) concentration falls at or below the LLOQ, i.e. `Φ((log(LLOQ) − predicted) / σ)`, where Φ is the standard normal CDF and σ the residual standard deviation. This means BLQ timepoints actively inform the slope estimate — they pull the λz estimate toward the true terminal slope — rather than being discarded. ### When to consider Tobit Tobit regression is most beneficial in the following situations: - A substantial fraction of terminal-phase observations — as a rough rule of thumb (not a documented PKNCA threshold), a fifth or more — are BLQ. With few or no BLQ points, OLS and Tobit give similar estimates. - The drug has **low tissue concentrations** or a short distribution half-life relative to the sampling schedule, making late timepoints frequently BLQ. - **Pediatric or geriatric studies** with limited sampling windows, where sparse late-time samples are often at or near the LLOQ. - The terminal slope is important for downstream parameters (AUC extrapolation, mean residence time) and precision matters for regulatory submissions. ### Using Tobit directly with `pk.calc.half.life()` (≥ 0.13.0) Tobit can be called directly through `pk.calc.half.life()` or used in the full `pk.nca()` pipeline (next section). For a direct call, supply the full concentration-time vector and the LLOQ (a scalar, or a vector with one value per timepoint): ```{r} # Simulated profile — concentrations at 16h and 24h are below the LLOQ of 0.1 mg/L conc_tobit <- c(0, 2.5, 4.8, 4.2, 2.9, 1.4, 0.6, 0.05, 0.01) time_tobit <- c(0, 0.5, 1, 2, 4, 8, 12, 16, 24) # OLS reference — the BLQ observations (16 and 24 h) are dropped is_blq <- conc_tobit > 0 & conc_tobit < 0.1 hl_ols <- pk.calc.half.life(conc = conc_tobit[!is_blq], time = time_tobit[!is_blq], tmax = 1, tlast = 12, hl_method = "log-linear") # Tobit — includes the BLQ observations as left-censored at the scalar lloq hl_tobit <- pk.calc.half.life(conc = conc_tobit, time = time_tobit, tmax = 1, tlast = 24, lloq = 0.1, hl_method = "tobit") # Compare data.frame( method = c("log-linear", "tobit"), lambda.z = c(hl_ols$lambda.z, hl_tobit$lambda.z), half.life = c(hl_ols$half.life, hl_tobit$half.life), n.points = c(hl_ols$lambda.z.n.points, hl_tobit$lambda.z.n.points) ) ``` With the two BLQ samples dropped, the log-linear fit uses the 4 remaining points after Tmax and gives a half-life of about 3.58 h. Tobit keeps the 16- and 24-hour samples as censored at the LLOQ, so its `lambda.z.n.points` of 6 counts them, and the steeper slope they imply gives a half-life of about 2.61 h. ### Tobit through the `pk.nca()` pipeline (≥ 0.13.0) To run Tobit across a full dataset, declare the LLOQ on the concentration object with the `lloq` argument of `PKNCAconc()` (a column name or a scalar) and select the method with the `hl_method` option. Every non-missing `lloq` value must be a finite number greater than zero; anything else is an error (`pknca_error_conc_invalid_lloq`). Missing values are allowed, and an entirely `NA` column is accepted. ```{r} # Subject 1 with an assay LLOQ of 0.5 mg/L d_conc_subset <- d_conc |> filter(subject == "1") d_dose_subset <- d_dose |> filter(subject == "1") o_conc_tobit <- PKNCAconc(d_conc_subset, conc ~ time | subject, lloq = 0.5) o_dose_tobit <- PKNCAdose(d_dose_subset, dose ~ time | subject, route = "extravascular") o_data_tobit <- PKNCAdata( o_conc_tobit, o_dose_tobit, intervals = data.frame(start = 0, end = Inf, half.life = TRUE), options = list(hl_method = "tobit") ) o_nca_tobit <- pk.nca(o_data_tobit) as.data.frame(o_nca_tobit) |> filter(PPTESTCD %in% c("half.life", "lambda.z.n.points_blq", "tobit_residual", "adj_tobit_residual")) |> select(PPTESTCD, PPORRES) ``` None of subject 1's terminal-phase samples fall below this LLOQ (`lambda.z.n.points_blq` is 0), so the Tobit estimate matches the OLS fit here — the example shows the wiring, and the Tobit-specific outputs appear alongside the usual half-life results. Tobit adds new output parameters to the results: | Parameter | Description | |---|---| | `lambda.z.n.points_blq` | Number of BLQ points included in the Tobit regression | | `tobit_residual` | Estimated residual standard deviation of the Tobit fit (log-concentration scale) | | `adj_tobit_residual` | Adjusted Tobit residual (analogous to adjusted R²) | ### Tobit-related options | Option | Default | Effect | |---|---|---| | `hl_method` | `"log-linear"` | `"log-linear"` (OLS) or `"tobit"` (censored-likelihood regression) | | `tobit_n_points_penalty` | `0` | Exponent in the window-selection criterion (residual × n^penalty); positive values favor shorter windows | | `tobit_optim_control` | `list()` | Control list forwarded to `stats::optim()` for the Tobit optimization | Window selection also differs from the log-linear fits: instead of the adjusted-R² rule, PKNCA fits each candidate terminal window and keeps the one that minimizes `tobit_residual × n_points^tobit_n_points_penalty`, with ties going to the larger window. `min.hl.points` and `max.hl.points` count only the above-LLOQ points of a Tobit window, and `min.hl.start.time` applies as it does for the log-linear method. ### Tobit with sparse PK (≥ 0.13.0) The Tobit half-life also works with sparse sampling. For sparse data the half-life is calculated on the arithmetic-mean profile, and that profile carries an LLOQ at each time: the median of the `lloq` values of its samples, ignoring missing ones. A time whose mean is set to zero because more than half of its samples are BLQ is therefore censored at that LLOQ rather than dropped. See the [sparse NCA chapter](sparse.qmd) for sparse analyses. See the vignette [Half-life calculation with Tobit regression](https://humanpred.github.io/pknca/articles/v06-half-life-calculation-tobit.html) for a detailed simulation comparison between OLS and Tobit estimates. --- ## `lambda.z.corrxy` — regression correlation `lambda.z.corrxy` (≥ 0.13.0) returns the Pearson correlation between time and log-concentration for the points used in the λz regression (numeric `NA` when it cannot be computed). It is a unitless diagnostic — values close to −1 indicate a strong, clean terminal slope. Like `r.squared` and `adj.r.squared`, it is returned only by the log-linear method. ```{r} hl_corrxy_interval <- data.frame( start = 0, end = Inf, half.life = TRUE, lambda.z = TRUE, lambda.z.corrxy = TRUE, r.squared = TRUE, adj.r.squared = TRUE ) o_nca_corrxy <- pk.nca(PKNCAdata(o_conc, o_dose, intervals = hl_corrxy_interval)) as.data.frame(o_nca_corrxy) |> filter(PPTESTCD %in% c("half.life", "lambda.z.corrxy", "r.squared", "adj.r.squared")) |> select(subject, PPTESTCD, PPORRES) |> tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES) |> arrange(subject) ``` --- ## Visualizing the chosen regression window Three functions extract what is needed to plot the fit, each from a `PKNCAresults` object or (by running only the half-life calculation) a `PKNCAdata` object: - `get_halflife_points()` gives one logical value per row of the concentration data: `TRUE` for points used in the fit, `FALSE` for points not used, and `NA` where no half-life was calculated for the interval. - `get_halflife_fit()` (≥ 0.13.0) gives the `intercept` and `slope` of the line `log(conc) = intercept + slope*time` and the `time_first`/`time_last` of the points used, one row per group and interval. Times are on the scale of the concentration data, not relative to the interval start, and the line is `NA` when the half-life was not calculated or was excluded. - `get_halflife_curve()` (≥ 0.13.0) gives concentrations along that line, at the times in `tout` or at `n` equally spaced times (default 50) spanning the points used. It extrapolates after the last point used by default (`extrapolate_later = TRUE`) and not before the first (`extrapolate_earlier = FALSE`); times outside the fit that are not extrapolated give `NA`. ```{r} get_halflife_fit(o_nca) |> filter(subject %in% c("1", "6")) ``` ```{r} # Mark the points used for each subject's fit d_conc_hl <- d_conc |> mutate(used_in_lambda_z = get_halflife_points(o_nca)) # The fitted line on a grid from 0 to 36 h: NA before the first point used, # extrapolated after the last one d_curve <- get_halflife_curve(o_nca, tout = seq(0, 36, by = 0.5), extrapolate_earlier = FALSE) ggplot(d_conc_hl |> filter(subject %in% c("1", "6"), conc > 0), aes(x = time, y = conc)) + geom_line(colour = "grey60") + geom_point(aes(colour = used_in_lambda_z), size = 3) + geom_line(data = d_curve |> filter(subject %in% c("1", "6"), !is.na(conc)), colour = "firebrick", linetype = "dashed") + scale_colour_manual(values = c("FALSE" = "grey60", "TRUE" = "firebrick"), labels = c("FALSE" = "not used in λz", "TRUE" = "used in λz"), name = NULL) + scale_y_log10() + facet_wrap(~subject, labeller = label_both) + labs(title = "Terminal regression window and λz fit (log scale)", x = "Time (h)", y = "Concentration (mg/L) [log scale]") + theme_minimal() ``` Subject 1's fit uses its last 3 points (from 9.05 h), while subject 6's uses 7 points (from 2.03 h). The dashed line starts at the first point used because `extrapolate_earlier = FALSE`, and it continues to 36 h because later times are extrapolated by default. For many subjects at once, the same functions feed a faceted plot or a `geom_line()` overlay in a QC report. --- ::: {.callout-note icon=false appearance="minimal"} **pkgdown reference:** [PKNCAconc()](https://humanpred.github.io/pknca/reference/PKNCAconc.html) · [PKNCAdose()](https://humanpred.github.io/pknca/reference/PKNCAdose.html) · [PKNCAdata()](https://humanpred.github.io/pknca/reference/PKNCAdata.html) · [pk.nca()](https://humanpred.github.io/pknca/reference/pk.nca.html) · [pk.calc.half.life()](https://humanpred.github.io/pknca/reference/pk.calc.half.life.html) · [PKNCA.options()](https://humanpred.github.io/pknca/reference/PKNCA.options.html) · [get_halflife_points()](https://humanpred.github.io/pknca/reference/get_halflife_points.html) · [get_halflife_fit()](https://humanpred.github.io/pknca/reference/get_halflife_fit.html) · [get_halflife_curve()](https://humanpred.github.io/pknca/reference/get_halflife_curve.html) · [exclude()](https://humanpred.github.io/pknca/reference/exclude.html) · [exclude_nca_min.hl.r.squared()](https://humanpred.github.io/pknca/reference/exclude_nca.html) · [exclude_nca_span.ratio()](https://humanpred.github.io/pknca/reference/exclude_nca.html) · [exclude_nca_min.hl.adj.r.squared()](https://humanpred.github.io/pknca/reference/exclude_nca.html) :::