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
|
#' Functions to transform and backtransform kinetic parameters for fitting
#'
#' The transformations are intended to map parameters that should only take on
#' restricted values to the full scale of real numbers. For kinetic rate
#' constants and other paramters that can only take on positive values, a
#' simple log transformation is used. For compositional parameters, such as the
#' formations fractions that should always sum up to 1 and can not be negative,
#' the \code{\link{ilr}} transformation is used.
#'
#' The transformation of sets of formation fractions is fragile, as it supposes
#' the same ordering of the components in forward and backward transformation.
#' This is no problem for the internal use in \code{\link{mkinfit}}.
#'
#' @aliases transform_odeparms backtransform_odeparms
#' @param parms Parameters of kinetic models as used in the differential
#' equations.
#' @param transparms Transformed parameters of kinetic models as used in the
#' fitting procedure.
#' @param mkinmod The kinetic model of class \code{\link{mkinmod}}, containing
#' the names of the model variables that are needed for grouping the
#' formation fractions before \code{\link{ilr}} transformation, the parameter
#' names and the information if the pathway to sink is included in the model.
#' @param transform_rates Boolean specifying if kinetic rate constants should
#' be transformed in the model specification used in the fitting for better
#' compliance with the assumption of normal distribution of the estimator. If
#' TRUE, also alpha and beta parameters of the FOMC model are
#' log-transformed, as well as k1 and k2 rate constants for the DFOP and HS
#' models and the break point tb of the HS model.
#' @param transform_fractions Boolean specifying if formation fractions
#' constants should be transformed in the model specification used in the
#' fitting for better compliance with the assumption of normal distribution
#' of the estimator. The default (TRUE) is to do transformations. The g
#' parameter of the DFOP and HS models are also transformed, as they can also
#' be seen as compositional data. The transformation used for these
#' transformations is the \code{\link{ilr}} transformation.
#' @return A vector of transformed or backtransformed parameters with the same
#' names as the original parameters.
#' @author Johannes Ranke
#' @examples
#'
#' SFO_SFO <- mkinmod(
#' parent = list(type = "SFO", to = "m1", sink = TRUE),
#' m1 = list(type = "SFO"))
#' # Fit the model to the FOCUS example dataset D using defaults
#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE)
#' fit.s <- summary(fit)
#' # Transformed and backtransformed parameters
#' print(fit.s$par, 3)
#' print(fit.s$bpar, 3)
#'
#' \dontrun{
#' # Compare to the version without transforming rate parameters
#' fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE, quiet = TRUE)
#' fit.2.s <- summary(fit.2)
#' print(fit.2.s$par, 3)
#' print(fit.2.s$bpar, 3)
#' }
#'
#' initials <- fit$start$value
#' names(initials) <- rownames(fit$start)
#' transformed <- fit$start_transformed$value
#' names(transformed) <- rownames(fit$start_transformed)
#' transform_odeparms(initials, SFO_SFO)
#' backtransform_odeparms(transformed, SFO_SFO)
#'
#' \dontrun{
#' # The case of formation fractions
#' SFO_SFO.ff <- mkinmod(
#' parent = list(type = "SFO", to = "m1", sink = TRUE),
#' m1 = list(type = "SFO"),
#' use_of_ff = "max")
#'
#' fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE)
#' fit.ff.s <- summary(fit.ff)
#' print(fit.ff.s$par, 3)
#' print(fit.ff.s$bpar, 3)
#' initials <- c("f_parent_to_m1" = 0.5)
#' transformed <- transform_odeparms(initials, SFO_SFO.ff)
#' backtransform_odeparms(transformed, SFO_SFO.ff)
#'
#' # And without sink
#' SFO_SFO.ff.2 <- mkinmod(
#' parent = list(type = "SFO", to = "m1", sink = FALSE),
#' m1 = list(type = "SFO"),
#' use_of_ff = "max")
#'
#'
#' fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_2006_D, quiet = TRUE)
#' fit.ff.2.s <- summary(fit.ff.2)
#' print(fit.ff.2.s$par, 3)
#' print(fit.ff.2.s$bpar, 3)
#' }
#'
#' @export transform_odeparms
transform_odeparms <- function(parms, mkinmod,
transform_rates = TRUE,
transform_fractions = TRUE)
{
# We need the model specification for the names of the model
# variables and the information on the sink
spec = mkinmod$spec
# Set up container for transformed parameters
transparms <- numeric(0)
# Do not transform initial values for state variables
state.ini.optim <- parms[grep("_0$", names(parms))]
transparms[names(state.ini.optim)] <- state.ini.optim
# Log transformation for rate constants if requested
k <- parms[grep("^k_", names(parms))]
k__iore <- parms[grep("^k__iore_", names(parms))]
k <- c(k, k__iore)
if (length(k) > 0) {
if(transform_rates) {
transparms[paste0("log_", names(k))] <- log(k)
} else transparms[names(k)] <- k
}
# Do not transform exponents in IORE models
N <- parms[grep("^N", names(parms))]
transparms[names(N)] <- N
# Go through state variables and apply isometric logratio transformation to
# formation fractions if requested
mod_vars = names(spec)
for (box in mod_vars) {
f <- parms[grep(paste("^f", box, sep = "_"), names(parms))]
if (length(f) > 0) {
if(transform_fractions) {
if (spec[[box]]$sink) {
trans_f <- ilr(c(f, 1 - sum(f)))
trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
transparms[trans_f_names] <- trans_f
} else {
if (length(f) > 1) {
trans_f <- ilr(f)
trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
transparms[trans_f_names] <- trans_f
}
}
} else {
transparms[names(f)] <- f
}
}
}
# Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2
# and HS parameter tb as well as logistic model parameters kmax, k0 and r if
# transformation of rates is requested
for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
if (!is.na(parms[pname])) {
if (transform_rates) {
transparms[paste0("log_", pname)] <- log(parms[pname])
} else {
transparms[pname] <- parms[pname]
}
}
}
# DFOP parameter g is treated as a fraction
if (!is.na(parms["g"])) {
g <- parms["g"]
if (transform_fractions) {
transparms["g_ilr"] <- ilr(c(g, 1 - g))
} else {
transparms["g"] <- g
}
}
return(transparms)
}
#' @describeIn transform_odeparms Backtransform the set of transformed parameters
#' @export backtransform_odeparms
backtransform_odeparms <- function(transparms, mkinmod,
transform_rates = TRUE,
transform_fractions = TRUE)
{
# We need the model specification for the names of the model
# variables and the information on the sink
spec = mkinmod$spec
# Set up container for backtransformed parameters
parms <- numeric(0)
# Do not transform initial values for state variables
state.ini.optim <- transparms[grep("_0$", names(transparms))]
parms[names(state.ini.optim)] <- state.ini.optim
# Exponential transformation for rate constants
if(transform_rates) {
trans_k <- transparms[grep("^log_k_", names(transparms))]
trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))]
trans_k = c(trans_k, trans_k__iore)
if (length(trans_k) > 0) {
k_names <- gsub("^log_k", "k", names(trans_k))
parms[k_names] <- exp(trans_k)
}
} else {
trans_k <- transparms[grep("^k_", names(transparms))]
parms[names(trans_k)] <- trans_k
trans_k__iore <- transparms[grep("^k__iore_", names(transparms))]
parms[names(trans_k__iore)] <- trans_k__iore
}
# Do not transform exponents in IORE models
N <- transparms[grep("^N", names(transparms))]
parms[names(N)] <- N
# Go through state variables and apply inverse isometric logratio transformation
mod_vars = names(spec)
for (box in mod_vars) {
# Get the names as used in the model
f_names = grep(paste("^f", box, sep = "_"), mkinmod$parms, value = TRUE)
# Get the formation fraction parameters
trans_f = transparms[grep(paste("^f", box, sep = "_"), names(transparms))]
if (length(trans_f) > 0) {
if(transform_fractions) {
f <- invilr(trans_f)
if (spec[[box]]$sink) {
parms[f_names] <- f[1:length(f)-1]
} else {
parms[f_names] <- f
}
} else {
parms[names(trans_f)] <- trans_f
}
}
}
# Transform parameters also for FOMC, DFOP, HS and logistic models
for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
if (transform_rates) {
pname_trans = paste0("log_", pname)
if (!is.na(transparms[pname_trans])) {
parms[pname] <- exp(transparms[pname_trans])
}
} else {
if (!is.na(transparms[pname])) {
parms[pname] <- transparms[pname]
}
}
}
# DFOP parameter g is treated as a fraction
if (!is.na(transparms["g_ilr"])) {
g_ilr <- transparms["g_ilr"]
parms["g"] <- invilr(g_ilr)[1]
}
if (!is.na(transparms["g"])) {
parms["g"] <- transparms["g"]
}
return(parms)
}
# vim: set ts=2 sw=2 expandtab:
|
