-
Notifications
You must be signed in to change notification settings - Fork 23
/
half.life.R
416 lines (413 loc) · 15.4 KB
/
half.life.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
#' Compute the half-life and associated parameters
#'
#' The terminal elimination half-life is estimated from the final points in the
#' concentration-time curve using semi-log regression (`log(conc)~time`)
#' with automated selection of the points for calculation (unless
#' `manually.selected.points` is `TRUE`).
#'
#' See the "Half-Life Calculation" vignette for more details on the calculation
#' methods used.
#'
#' @details If `manually.selected.points` is `FALSE` (default), the
#' half-life is calculated by computing the best fit line for all points at or
#' after tmax (based on the value of `allow.tmax.in.half.life`). The best
#' half-life is chosen by the following rules in order:
#'
#' \itemize{
#' \item{At least `min.hl.points` points included}
#' \item{A `lambda.z` > 0 and at the same time the best adjusted r-squared
#' (within `adj.r.squared.factor`)}
#' \item{The one with the most points included}
#' }
#'
#' If `manually.selected.points` is `TRUE`, the `conc`
#' and `time` data are used as-is without any form of selection for
#' the best-fit half-life.
#'
#' @inheritParams assert_conc_time
#' @inheritParams choose_interval_method
#' @inheritParams PKNCA.choose.option
#' @param tmax Time of maximum concentration (will be calculated and
#' included in the return data frame if not given)
#' @param tlast Time of last concentration above the limit of
#' quantification (will be calculated and included in the return data
#' frame if not given)
#' @param manually.selected.points Have the input points (`conc` and
#' `time`) been manually selected? The impact of setting this to
#' `TRUE` is that no selection for the best points will be done. When
#' `TRUE`, this option causes the options of `adj.r.squared.factor`,
#' `min.hl.points`, and `allow.tmax.in.half.life` to be ignored.
#' @param min.hl.points The minimum number of points that must be
#' included to calculate the half-life
#' @param adj.r.squared.factor The allowance in adjusted r-squared for
#' adding another point.
#' @param conc.blq See [clean.conc.blq()]
#' @param conc.na See [clean.conc.na()]
#' @param check Run [assert_conc_time()],
#' [clean.conc.blq()], and [clean.conc.na()]?
#' @param first.tmax See [pk.calc.tmax()].
#' @param allow.tmax.in.half.life Allow the concentration point for tmax
#' to be included in the half-life slope calculation.
#' @return A data frame with one row and columns for
#' \describe{
#' \item{tmax}{Time of maximum observed concentration (only included
#' if not given as an input)}
#' \item{tlast}{Time of last observed concentration above the LOQ (only
#' included if not given as an input)}
#' \item{r.squared}{coefficient of determination}
#' \item{adj.r.squared}{adjusted coefficient of determination}
#' \item{lambda.z}{elimination rate}
#' \item{lambda.z.time.first}{first time for half-life calculation}
#' \item{lambda.z.n.points}{number of points in half-life calculation}
#' \item{clast.pred}{Concentration at tlast as predicted by the half-life
#' line}
#' \item{half.life}{half-life}
#' \item{span.ratio}{span ratio [ratio of half-life to time used for
#' half-life calculation}
#' }
#' @references
#'
#' Gabrielsson J, Weiner D. "Section 2.8.4 Strategies for estimation of
#' lambda-z." Pharmacokinetic & Pharmacodynamic Data Analysis: Concepts
#' and Applications, 4th Edition. Stockholm, Sweden: Swedish
#' Pharmaceutical Press, 2000. 167-9.
#' @family NCA parameter calculations
#' @export
pk.calc.half.life <- function(conc, time, tmax, tlast,
manually.selected.points=FALSE,
options=list(),
min.hl.points=NULL,
adj.r.squared.factor=NULL,
conc.blq=NULL,
conc.na=NULL,
first.tmax=NULL,
allow.tmax.in.half.life=NULL,
check=TRUE) {
# Check inputs
min.hl.points <-
PKNCA.choose.option(
name="min.hl.points", value=min.hl.points, options=options
)
adj.r.squared.factor <-
PKNCA.choose.option(
name="adj.r.squared.factor", value=adj.r.squared.factor, options=options
)
conc.blq <-
PKNCA.choose.option(
name="conc.blq", value=conc.blq, options=options
)
conc.na <-
PKNCA.choose.option(
name="conc.na", value=conc.na, options=options
)
first.tmax <-
PKNCA.choose.option(
name="first.tmax", value=first.tmax, options=options
)
allow.tmax.in.half.life <-
PKNCA.choose.option(
name="allow.tmax.in.half.life", value=allow.tmax.in.half.life, options=options
)
if (check) {
assert_conc_time(conc = conc, time = time)
data <- clean.conc.blq(conc, time, conc.blq=conc.blq, conc.na=conc.na)
} else {
data <- data.frame(conc, time)
}
# if (inherits(data$conc, "units")) {
# conc_units <- units(data$conc)
# } else {
conc_units <- NULL
#}
data$log_conc <- log(data$conc)
# as.numeric() to handle units objects
data <- data[as.numeric(data$conc) > 0,]
# Prepare the return values
ret <- data.frame(
# Terminal elimination slope
lambda.z=NA,
# R-squared of terminal elimination slope
r.squared=NA,
# Adjusted r-squared of terminal elimination slope
adj.r.squared=NA,
# First time point used in the slope estimation (for plotting later)
lambda.z.time.first=NA,
# Number of points in the half-life estimate
lambda.z.n.points=NA,
# Concentration at Tlast predicted by the half-life
clast.pred=NA,
# Half-life
half.life=NA,
# T1/2 span ratio
span.ratio=NA)
ret_replacements <-
c("lambda.z", "r.squared", "adj.r.squared", "lambda.z.time.first",
"lambda.z.n.points", "clast.pred", "half.life", "span.ratio")
if (missing(tmax)) {
ret$tmax <-
pk.calc.tmax(
data$conc, data$time, first.tmax=first.tmax, check=FALSE
)
} else {
ret$tmax <- tmax
}
if (missing(tlast)) {
ret$tlast <- pk.calc.tlast(data$conc, data$time, check=FALSE)
} else {
ret$tlast <- tlast
}
# Data frame to use for computation of half-life
if (allow.tmax.in.half.life) {
# as.numeric is for units handling
dfK <- data[as.numeric(data$time) >= as.numeric(ret$tmax), ]
} else {
# as.numeric is for units handling
dfK <- data[as.numeric(data$time) > as.numeric(ret$tmax), ]
}
if (manually.selected.points) {
if (nrow(data) > 0) {
fit <- fit_half_life(data=data, tlast=ret$tlast, conc_units=conc_units)
ret[,ret_replacements] <- fit[,ret_replacements]
} else {
warning("No data to manually fit for half-life (all concentrations may be 0 or excluded)")
ret <-
structure(
ret,
exclude="No data to manually fit for half-life (all concentrations may be 0 or excluded)"
)
}
} else if (nrow(dfK) >= min.hl.points) {
# If we have enough data to estimate a slope, then
half_lives_for_selection <-
data.frame(
r.squared=-Inf,
adj.r.squared=-Inf,
clast.pred=NA_real_,
lambda.z=-Inf,
lambda.z.n.points=NA_integer_,
lambda.z.time.first=dfK$time,
log_conc=dfK$log_conc,
span.ratio=NA_real_,
half.life=NA_real_
)
half_lives_for_selection <-
half_lives_for_selection[order(-half_lives_for_selection$lambda.z.time.first), ]
for(i in min.hl.points:nrow(half_lives_for_selection)) {
# Fit the terminal slopes until the adjusted r-squared value
# is not improving (or it only gets worse by a small factor).
fit <-
fit_half_life(
data=
data.frame(
# pass in the conc so that we can use its units, if applicable
log_conc=half_lives_for_selection$log_conc[1:i],
time=half_lives_for_selection$lambda.z.time.first[1:i]
),
tlast=ret$tlast,
conc_units=conc_units
)
half_lives_for_selection[i,names(fit)] <- fit
}
# Find the best model
mask_best <-
half_lives_for_selection$lambda.z > 0 &
if (min.hl.points == 2 & nrow(half_lives_for_selection) == 2) {
rlang::warn(
message = "2 points used for half-life calculation",
class = "pknca_halflife_2points"
)
TRUE
} else {
half_lives_for_selection$adj.r.squared >
(max(half_lives_for_selection$adj.r.squared, na.rm=TRUE) - adj.r.squared.factor)
}
# Missing values are not the best
mask_best[is.na(mask_best)] <- FALSE
if (sum(mask_best) > 1) {
# If more than one models qualify, choose the one with the
# most points used.
mask_best <-
(mask_best &
half_lives_for_selection$lambda.z.n.points == max(half_lives_for_selection$lambda.z.n.points[mask_best]))
}
# If the half-life fit, set all associated parameters
if (any(mask_best)) {
# Put in all the computed values
ret[,ret_replacements] <- half_lives_for_selection[mask_best, ret_replacements]
}
} else {
attr(ret, "exclude") <-
sprintf(
"Too few points for half-life calculation (min.hl.points=%g with only %g points)",
min.hl.points, nrow(dfK)
)
rlang::warn(
message = attr(ret, "exclude"),
class = "pknca_halflife_too_few_points"
)
}
# Drop the inputs of tmax and tlast, if given.
if (!missing(tmax))
ret$tmax <- NULL
if (!missing(tlast))
ret$tlast <- NULL
ret
}
#' Perform the half-life fit given the data. The function simply fits
#' the data without any validation. No selection of points or any other
#' components are done.
#'
#' @param data The data to fit. Must have two columns named "log_conc"
#' and "time"
#' @param tlast The time of last observed concentration above the limit
#' of quantification.
#' @param conc_units NULL or the units to set for concentration measures
#' @return A data.frame with one row and columns named "r.squared",
#' "adj.r.squared", "PROB", "lambda.z", "clast.pred",
#' "lambda.z.n.points", "half.life", "span.ratio"
#' @seealso [pk.calc.half.life()]
fit_half_life <- function(data, tlast, conc_units) {
fit <- stats::.lm.fit(x=cbind(1, data$time), y=data$log_conc)
# unit handling
# if (inherits(tlast, "units")) {
# time_units <- units(tlast)
# } else if (inherits(tlast, "mixed_units")) {
# time_units <- units(units::as_units(tlast))
# } else {
time_units <- NULL
# }
# if (!is.null(time_units)) {
# inverse_time_units <- time_units
# inverse_time_units$numerator <- time_units$denominator
# inverse_time_units$denominator <- time_units$numerator
# } else {
inverse_time_units <- NULL
# }
# as.numeric is so that it works for units objects
r_squared <- 1 - as.numeric(sum(fit$residuals^2))/as.numeric(sum((data$log_conc - mean(data$log_conc))^2))
clast_pred <- exp(sum(fit$coefficients*c(1, as.numeric(tlast))))
# if (!is.null(conc_units)) {
# clast_pred <- units::set_units(clast_pred, conc_units, mode="standard")
# }
lambda_z <- -fit$coefficients[2]
# if (!is.null(inverse_time_units)) {
# lambda_z <- units::set_units(lambda_z, inverse_time_units, mode="standard")
# }
ret <-
data.frame(
r.squared=r_squared,
adj.r.squared=adj.r.squared(r_squared, nrow(data)),
lambda.z=lambda_z,
clast.pred=clast_pred,
lambda.z.time.first=min(data$time, na.rm=TRUE),
lambda.z.n.points=nrow(data)
)
ret$half.life <- log(2)/ret$lambda.z
ret$span.ratio <- (max(data$time) - min(data$time))/ret$half.life
ret
}
# Add the column to the interval specification
add.interval.col("half.life",
FUN="pk.calc.half.life",
values=c(FALSE, TRUE),
unit_type="time",
pretty_name="Half-life",
desc="The (terminal) half-life",
depends=c("tmax", "tlast"))
PKNCA.set.summary(
name="half.life",
description="arithmetic mean and standard deviation",
point=business.mean,
spread=business.sd
)
add.interval.col("r.squared",
FUN=NA,
values=c(FALSE, TRUE),
unit_type="unitless",
pretty_name="$r^2$",
desc="The r^2 value of the half-life calculation",
depends="half.life")
PKNCA.set.summary(
name="r.squared",
description="arithmetic mean and standard deviation",
point=business.mean,
spread=business.sd
)
add.interval.col("adj.r.squared",
FUN=NA,
values=c(FALSE, TRUE),
unit_type="unitless",
pretty_name="$r^2_{adj}$",
desc="The adjusted r^2 value of the half-life calculation",
depends="half.life")
PKNCA.set.summary(
name="adj.r.squared",
description="arithmetic mean and standard deviation",
point=business.mean,
spread=business.sd
)
add.interval.col("lambda.z",
FUN=NA,
values=c(FALSE, TRUE),
unit_type="inverse_time",
pretty_name="$\\lambda_z$",
desc="The elimination rate of the terminal half-life",
depends="half.life")
PKNCA.set.summary(
name="lambda.z",
description="geometric mean and geometric coefficient of variation",
point=business.geomean,
spread=business.geocv
)
add.interval.col("lambda.z.time.first",
FUN=NA,
values=c(FALSE, TRUE),
unit_type="time",
pretty_name="First time for $\\lambda_z$",
desc="The first time point used for the calculation of half-life",
depends="half.life")
PKNCA.set.summary(
name="lambda.z.time.first",
description="median and range",
point=business.median,
spread=business.range
)
add.interval.col("lambda.z.n.points",
FUN=NA,
values=c(FALSE, TRUE),
unit_type="count",
pretty_name="Number of points used for lambda_z",
desc="The number of points used for the calculation of half-life",
depends="half.life")
PKNCA.set.summary(
name="lambda.z.n.points",
description="median and range",
point=business.median,
spread=business.range
)
add.interval.col("clast.pred",
FUN=NA,
values=c(FALSE, TRUE),
unit_type="conc",
pretty_name="Clast,pred",
desc="The concentration at Tlast as predicted by the half-life",
depends="half.life")
PKNCA.set.summary(
name="clast.pred",
description="geometric mean and geometric coefficient of variation",
point=business.geomean,
spread=business.geocv
)
add.interval.col("span.ratio",
FUN=NA,
values=c(FALSE, TRUE),
unit_type="fraction",
pretty_name="Span ratio",
desc="The ratio of the half-life to the duration used for half-life calculation",
depends="half.life")
PKNCA.set.summary(
name="span.ratio",
description="geometric mean and geometric coefficient of variation",
point=business.geomean,
spread=business.geocv
)