-
Notifications
You must be signed in to change notification settings - Fork 4
/
endpoints.R
255 lines (238 loc) · 9.62 KB
/
endpoints.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
#' Function to calculate endpoints for further use from kinetic models fitted
#' with mkinfit
#'
#' This function calculates DT50 and DT90 values as well as formation fractions
#' from kinetic models fitted with mkinfit. If the SFORB model was specified
#' for one of the parents or metabolites, the Eigenvalues are returned. These
#' are equivalent to the rate constants of the DFOP model, but with the
#' advantage that the SFORB model can also be used for metabolites.
#'
#' Additional DT50 values are calculated from the FOMC DT90 and k1 and k2 from
#' HS and DFOP, as well as from Eigenvalues b1 and b2 of any SFORB models
#'
#' @param fit An object of class [mkinfit], [nlme.mmkin] or [saem.mmkin], or
#' another object that has list components mkinmod containing an [mkinmod]
#' degradation model, and two numeric vectors, bparms.optim and bparms.fixed,
#' that contain parameter values for that model.
#' @param covariates Numeric vector with covariate values for all variables in
#' any covariate models in the object. If given, it overrides 'covariate_quantile'.
#' @param covariate_quantile This argument only has an effect if the fitted
#' object has covariate models. If so, the default is to show endpoints
#' for the median of the covariate values (50th percentile).
#' @importFrom stats optimize
#' @return A list with a matrix of dissipation times named distimes, and, if
#' applicable, a vector of formation fractions named ff and, if the SFORB model
#' was in use, a vector of eigenvalues of these SFORB models, equivalent to
#' DFOP rate constants
#' @note The function is used internally by [summary.mkinfit],
#' [summary.nlme.mmkin] and [summary.saem.mmkin].
#' @author Johannes Ranke
#' @examples
#'
#' fit <- mkinfit("FOMC", FOCUS_2006_C, quiet = TRUE)
#' endpoints(fit)
#' \dontrun{
#' fit_2 <- mkinfit("DFOP", FOCUS_2006_C, quiet = TRUE)
#' endpoints(fit_2)
#' fit_3 <- mkinfit("SFORB", FOCUS_2006_C, quiet = TRUE)
#' endpoints(fit_3)
#' }
#'
#' @export
endpoints <- function(fit, covariates = NULL, covariate_quantile = 0.5) {
mkinmod <- fit$mkinmod
obs_vars <- names(mkinmod$spec)
if (!is.null(fit$covariate_models)) {
if (is.null(covariates)) {
covariates = as.data.frame(
apply(fit$covariates, 2, quantile,
covariate_quantile, simplify = FALSE))
} else {
covariate_m <- matrix(covariates, byrow = TRUE)
colnames(covariate_m) <- names(covariates)
rownames(covariate_m) <- "User"
covariates <- as.data.frame(covariate_m)
}
degparms_trans <- parms(fit, covariates = covariates)[, 1]
if (inherits(fit, "saem.mmkin") & (fit$transformations == "saemix")) {
degparms <- degparms_trans
} else {
degparms <- backtransform_odeparms(degparms_trans,
fit$mkinmod,
transform_rates = fit$transform_rates,
transform_fractions = fit$transform_fractions)
}
} else {
degparms <- c(fit$bparms.optim, fit$bparms.fixed)
}
# Set up object to return
ep <- list()
ep$covariates <- covariates
ep$ff <- vector()
ep$SFORB <- vector()
ep$distimes <- data.frame(
DT50 = rep(NA, length(obs_vars)),
DT90 = rep(NA, length(obs_vars)),
row.names = obs_vars)
for (obs_var in obs_vars) {
type = names(mkinmod$map[[obs_var]])[1]
# Get formation fractions if directly fitted, and calculate remaining fraction to sink
f_names = grep(paste("^f", obs_var, sep = "_"), names(degparms), value=TRUE)
if (length(f_names) > 0) {
f_values = degparms[f_names]
f_to_sink = 1 - sum(f_values)
names(f_to_sink) = ifelse(type == "SFORB",
paste(obs_var, "free", "sink", sep = "_"),
paste(obs_var, "sink", sep = "_"))
for (f_name in f_names) {
ep$ff[[sub("f_", "", sub("_to_", "_", f_name))]] = f_values[[f_name]]
}
ep$ff = append(ep$ff, f_to_sink)
}
# Get the rest
if (type == "SFO") {
k_names = grep(paste("^k", obs_var, sep="_"), names(degparms), value=TRUE)
k_tot = sum(degparms[k_names])
DT50 = log(2)/k_tot
DT90 = log(10)/k_tot
if (mkinmod$use_of_ff == "min" && length(obs_vars) > 1) {
for (k_name in k_names)
{
ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
}
}
}
if (type == "FOMC") {
alpha = degparms["alpha"]
beta = degparms["beta"]
DT50 = beta * (2^(1/alpha) - 1)
DT90 = beta * (10^(1/alpha) - 1)
DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
ep$distimes[obs_var, c("DT50back")] = DT50_back
}
if (type == "IORE") {
k_names = grep(paste("^k__iore", obs_var, sep="_"), names(degparms), value=TRUE)
k_tot = sum(degparms[k_names])
# From the NAFTA kinetics guidance, p. 5
n = degparms[paste("N", obs_var, sep = "_")]
k = k_tot
# Use the initial concentration of the parent compound
source_name = mkinmod$map[[1]][[1]]
c0 = degparms[paste(source_name, "0", sep = "_")]
alpha = 1 / (n - 1)
beta = (c0^(1 - n))/(k * (n - 1))
DT50 = beta * (2^(1/alpha) - 1)
DT90 = beta * (10^(1/alpha) - 1)
DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
ep$distimes[obs_var, c("DT50back")] = DT50_back
if (mkinmod$use_of_ff == "min") {
for (k_name in k_names)
{
ep$ff[[sub("k_", "", k_name)]] = degparms[[k_name]] / k_tot
}
}
}
if (type == "DFOP") {
k1 = degparms["k1"]
k2 = degparms["k2"]
g = degparms["g"]
f <- function(log_t, x) {
t <- exp(log_t)
fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
(fraction - (1 - x/100))^2
}
DT50_k1 = log(2)/k1
DT50_k2 = log(2)/k2
DT90_k1 = log(10)/k1
DT90_k2 = log(10)/k2
DT50 <- try(exp(optimize(f, c(log(DT50_k1), log(DT50_k2)), x=50)$minimum),
silent = TRUE)
DT90 <- try(exp(optimize(f, c(log(DT90_k1), log(DT90_k2)), x=90)$minimum),
silent = TRUE)
if (inherits(DT50, "try-error")) DT50 = NA
if (inherits(DT90, "try-error")) DT90 = NA
DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
ep$distimes[obs_var, c("DT50back")] = DT50_back
ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
}
if (type == "HS") {
k1 = degparms["k1"]
k2 = degparms["k2"]
tb = degparms["tb"]
DTx <- function(x) {
DTx.a <- (log(100/(100 - x)))/k1
DTx.b <- tb + (log(100/(100 - x)) - k1 * tb)/k2
if (DTx.a < tb) DTx <- DTx.a
else DTx <- DTx.b
return(DTx)
}
DT50 <- DTx(50)
DT90 <- DTx(90)
DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
DT50_k1 = log(2)/k1
DT50_k2 = log(2)/k2
ep$distimes[obs_var, c("DT50back")] = DT50_back
ep$distimes[obs_var, c("DT50_k1")] = DT50_k1
ep$distimes[obs_var, c("DT50_k2")] = DT50_k2
}
if (type == "SFORB") {
# FOCUS kinetics (2006), p. 60 f
k_out_names = grep(paste("^k", obs_var, "free", sep="_"), names(degparms), value=TRUE)
k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
k_1output = sum(degparms[k_out_names])
k_12 = degparms[paste("k", obs_var, "free", "bound", sep="_")]
k_21 = degparms[paste("k", obs_var, "bound", "free", sep="_")]
sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 - k_1output * k_21)
b1 = 0.5 * (k_12 + k_21 + k_1output) + sqrt_exp
b2 = 0.5 * (k_12 + k_21 + k_1output) - sqrt_exp
g = (k_12 + k_21 - b1)/(b2 - b1)
DT50_b1 = log(2)/b1
DT50_b2 = log(2)/b2
DT90_b1 = log(10)/b1
DT90_b2 = log(10)/b2
SFORB_fraction = function(t) {
g * exp(-b1 * t) + (1 - g) * exp(-b2 * t)
}
f_50 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.5)^2
log_DT50 <- try(optimize(f_50, c(log(DT50_b1), log(DT50_b2)))$minimum,
silent = TRUE)
f_90 <- function(log_t) (SFORB_fraction(exp(log_t)) - 0.1)^2
log_DT90 <- try(optimize(f_90, c(log(DT90_b1), log(DT90_b2)))$minimum,
silent = TRUE)
DT50 = if (inherits(log_DT50, "try-error")) NA
else exp(log_DT50)
DT90 = if (inherits(log_DT90, "try-error")) NA
else exp(log_DT90)
DT50_back = DT90 / (log(10)/log(2)) # Backcalculated DT50 as recommended in FOCUS 2011
for (k_out_name in k_out_names)
{
ep$ff[[sub("k_", "", k_out_name)]] = degparms[[k_out_name]] / k_1output
}
# Return the eigenvalues for comparison with DFOP rate constants
ep$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
ep$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
# Return g for comparison with DFOP
ep$SFORB[[paste(obs_var, "g", sep="_")]] = g
ep$distimes[obs_var, c("DT50back")] = DT50_back
ep$distimes[obs_var, c(paste("DT50", obs_var, "b1", sep = "_"))] = DT50_b1
ep$distimes[obs_var, c(paste("DT50", obs_var, "b2", sep = "_"))] = DT50_b2
}
if (type == "logistic") {
# FOCUS kinetics (2014) p. 67
kmax = degparms["kmax"]
k0 = degparms["k0"]
r = degparms["r"]
DT50 = (1/r) * log(1 - ((kmax/k0) * (1 - 2^(r/kmax))))
DT90 = (1/r) * log(1 - ((kmax/k0) * (1 - 10^(r/kmax))))
DT50_k0 = log(2)/k0
DT50_kmax = log(2)/kmax
ep$distimes[obs_var, c("DT50_k0")] = DT50_k0
ep$distimes[obs_var, c("DT50_kmax")] = DT50_kmax
}
ep$distimes[obs_var, c("DT50", "DT90")] = c(DT50, DT90)
}
if (length(ep$ff) == 0) ep$ff <- NULL
if (length(ep$SFORB) == 0) ep$SFORB <- NULL
return(ep)
}