jranke / mkin
1
#' Functions to transform and backtransform kinetic parameters for fitting
2
#'
3
#' The transformations are intended to map parameters that should only take on
4
#' restricted values to the full scale of real numbers. For kinetic rate
5
#' constants and other parameters that can only take on positive values, a
6
#' simple log transformation is used. For compositional parameters, such as the
7
#' formations fractions that should always sum up to 1 and can not be negative,
8
#' the [ilr] transformation is used.
9
#'
10
#' The transformation of sets of formation fractions is fragile, as it supposes
11
#' the same ordering of the components in forward and backward transformation.
12
#' This is no problem for the internal use in [mkinfit].
13
#'
14
#' @param parms Parameters of kinetic models as used in the differential
15
#'   equations.
16
#' @param transparms Transformed parameters of kinetic models as used in the
17
#'   fitting procedure.
18
#' @param mkinmod The kinetic model of class [mkinmod], containing
19
#'   the names of the model variables that are needed for grouping the
20
#'   formation fractions before [ilr] transformation, the parameter
21
#'   names and the information if the pathway to sink is included in the model.
22
#' @param transform_rates Boolean specifying if kinetic rate constants should
23
#'   be transformed in the model specification used in the fitting for better
24
#'   compliance with the assumption of normal distribution of the estimator. If
25
#'   TRUE, also alpha and beta parameters of the FOMC model are
26
#'   log-transformed, as well as k1 and k2 rate constants for the DFOP and HS
27
#'   models and the break point tb of the HS model.
28
#' @param transform_fractions Boolean specifying if formation fractions
29
#'   constants should be transformed in the model specification used in the
30
#'   fitting for better compliance with the assumption of normal distribution
31
#'   of the estimator. The default (TRUE) is to do transformations.
32
#'   The g parameter of the DFOP model is also seen as a fraction.
33
#'   If a single fraction is transformed (g parameter of DFOP or only a single
34
#'   target variable e.g. a single metabolite plus a pathway to sink), a
35
#'   logistic transformation is used [stats::qlogis()]. In other cases, i.e. if
36
#'   two or more formation fractions need to be transformed whose sum cannot
37
#'   exceed one, the [ilr] transformation is used.
38
#' @return A vector of transformed or backtransformed parameters
39
#' @importFrom stats plogis qlogis
40
#' @author Johannes Ranke
41
#' @examples
42
#'
43
#' SFO_SFO <- mkinmod(
44
#'   parent = list(type = "SFO", to = "m1", sink = TRUE),
45
#'   m1 = list(type = "SFO"), use_of_ff = "min")
46
#'
47
#' # Fit the model to the FOCUS example dataset D using defaults
48
#' FOCUS_D <- subset(FOCUS_2006_D, value != 0) # remove zero values to avoid warning
49
#' fit <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE)
50
#' fit.s <- summary(fit)
51
#' # Transformed and backtransformed parameters
52
#' print(fit.s$par, 3)
53
#' print(fit.s$bpar, 3)
54
#'
55
#' \dontrun{
56
#' # Compare to the version without transforming rate parameters (does not work
57
#' # with analytical solution, we get NA values for m1 in predictions)
58
#' fit.2 <- mkinfit(SFO_SFO, FOCUS_D, transform_rates = FALSE,
59
#'   solution_type = "deSolve", quiet = TRUE)
60
#' fit.2.s <- summary(fit.2)
61
#' print(fit.2.s$par, 3)
62
#' print(fit.2.s$bpar, 3)
63
#' }
64
#'
65
#' initials <- fit$start$value
66
#' names(initials) <- rownames(fit$start)
67
#' transformed <- fit$start_transformed$value
68
#' names(transformed) <- rownames(fit$start_transformed)
69
#' transform_odeparms(initials, SFO_SFO)
70
#' backtransform_odeparms(transformed, SFO_SFO)
71
#'
72
#' \dontrun{
73
#' # The case of formation fractions (this is now the default)
74
#' SFO_SFO.ff <- mkinmod(
75
#'   parent = list(type = "SFO", to = "m1", sink = TRUE),
76
#'   m1 = list(type = "SFO"),
77
#'   use_of_ff = "max")
78
#'
79
#' fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_D, quiet = TRUE)
80
#' fit.ff.s <- summary(fit.ff)
81
#' print(fit.ff.s$par, 3)
82
#' print(fit.ff.s$bpar, 3)
83
#' initials <- c("f_parent_to_m1" = 0.5)
84
#' transformed <- transform_odeparms(initials, SFO_SFO.ff)
85
#' backtransform_odeparms(transformed, SFO_SFO.ff)
86
#'
87
#' # And without sink
88
#' SFO_SFO.ff.2 <- mkinmod(
89
#'   parent = list(type = "SFO", to = "m1", sink = FALSE),
90
#'   m1 = list(type = "SFO"),
91
#'   use_of_ff = "max")
92
#'
93
#'
94
#' fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_D, quiet = TRUE)
95
#' fit.ff.2.s <- summary(fit.ff.2)
96
#' print(fit.ff.2.s$par, 3)
97
#' print(fit.ff.2.s$bpar, 3)
98
#' }
99
#'
100
#' @export transform_odeparms
101
transform_odeparms <- function(parms, mkinmod,
102
  transform_rates = TRUE, transform_fractions = TRUE)
103
{
104
  # We need the model specification for the names of the model
105
  # variables and the information on the sink
106 1
  spec = mkinmod$spec
107

108
  # Set up container for transformed parameters
109 1
  transparms <- numeric(0)
110

111
  # Do not transform initial values for state variables
112 1
  state.ini.optim <- parms[grep("_0$", names(parms))]
113 1
  transparms[names(state.ini.optim)] <- state.ini.optim
114

115
  # Log transformation for rate constants if requested
116 1
  k <- parms[grep("^k_", names(parms))]
117 1
  k__iore <- parms[grep("^k__iore_", names(parms))]
118 1
  k <- c(k, k__iore)
119 1
  if (length(k) > 0) {
120 1
    if(transform_rates) {
121 1
      transparms[paste0("log_", names(k))] <- log(k)
122 1
    } else transparms[names(k)] <- k
123
  }
124

125
  # Do not transform exponents in IORE models
126 1
  N <- parms[grep("^N", names(parms))]
127 1
  transparms[names(N)] <- N
128

129
  # Go through state variables and transform formation fractions if requested
130 1
  mod_vars = names(spec)
131 1
  for (box in mod_vars) {
132 1
    f <- parms[grep(paste("^f", box, sep = "_"), names(parms))]
133

134 1
    if (length(f) > 0) {
135 1
      if(transform_fractions) {
136 1
        if (spec[[box]]$sink) {
137 1
          if (length(f) == 1) {
138 1
            trans_f_name <- paste("f", box, "qlogis", sep = "_")
139 1
            transparms[trans_f_name] <- qlogis(f)
140
          } else {
141 1
            trans_f <- ilr(c(f, 1 - sum(f)))
142 1
            trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
143 1
            transparms[trans_f_names] <- trans_f
144
          }
145
        } else {
146 1
          if (length(f) > 1) {
147 1
            trans_f <- ilr(f)
148 1
            trans_f_names <- paste("f", box, "ilr", 1:length(trans_f), sep = "_")
149 1
            transparms[trans_f_names] <- trans_f
150
          }
151
        }
152
      } else {
153 1
        transparms[names(f)] <- f
154
      }
155
    }
156
  }
157

158
  # Transform also FOMC parameters alpha and beta, DFOP and HS rates k1 and k2
159
  # and HS parameter tb as well as logistic model parameters kmax, k0 and r if
160
  # transformation of rates is requested
161 1
  for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
162 1
    if (!is.na(parms[pname])) {
163 1
      if (transform_rates) {
164 1
        transparms[paste0("log_", pname)] <- log(parms[pname])
165
      } else {
166 0
        transparms[pname] <- parms[pname]
167
      }
168
    }
169
  }
170

171
  # DFOP parameter g is treated as a fraction
172 1
  if (!is.na(parms["g"])) {
173 1
    g <- parms["g"]
174 1
    if (transform_fractions) {
175 1
      transparms["g_qlogis"] <- qlogis(g)
176
    } else {
177 0
      transparms["g"] <- g
178
    }
179
  }
180

181 1
  return(transparms)
182
}
183

184
#' @rdname transform_odeparms
185
#' @export backtransform_odeparms
186
backtransform_odeparms <- function(transparms, mkinmod,
187
                                   transform_rates = TRUE,
188
                                   transform_fractions = TRUE)
189
{
190
  # We need the model specification for the names of the model
191
  # variables and the information on the sink
192 1
  spec = mkinmod$spec
193

194
  # Set up container for backtransformed parameters
195 1
  parms <- numeric(0)
196

197
  # Do not transform initial values for state variables
198 1
  state.ini.optim <- transparms[grep("_0$", names(transparms))]
199 1
  parms[names(state.ini.optim)] <- state.ini.optim
200

201
  # Exponential transformation for rate constants
202 1
  if(transform_rates) {
203 1
    trans_k <- transparms[grep("^log_k_", names(transparms))]
204 1
    trans_k__iore <- transparms[grep("^log_k__iore_", names(transparms))]
205 1
    trans_k = c(trans_k, trans_k__iore)
206 1
    if (length(trans_k) > 0) {
207 1
      k_names <- gsub("^log_k", "k", names(trans_k))
208 1
      parms[k_names] <- exp(trans_k)
209
    }
210
  } else {
211 1
    trans_k <- transparms[grep("^k_", names(transparms))]
212 1
    parms[names(trans_k)] <- trans_k
213 1
    trans_k__iore <- transparms[grep("^k__iore_", names(transparms))]
214 1
    parms[names(trans_k__iore)] <- trans_k__iore
215
  }
216

217
  # Do not transform exponents in IORE models
218 1
  N <- transparms[grep("^N", names(transparms))]
219 1
  parms[names(N)] <- N
220

221
  # Go through state variables and apply inverse transformations to formation fractions
222 1
  mod_vars = names(spec)
223 1
  for (box in mod_vars) {
224
    # Get the names as used in the model
225 1
    f_names = grep(paste("^f", box, sep = "_"), mkinmod$parms, value = TRUE)
226

227
    # Get the formation fraction parameters
228 1
    trans_f = transparms[grep(paste("^f", box, sep = "_"), names(transparms))]
229 1
    if (length(trans_f) > 0) {
230 1
      if(transform_fractions) {
231 1
        if (any(grepl("qlogis", names(trans_f)))) {
232 1
          parms[f_names] <- plogis(trans_f)
233
        } else {
234 1
          f <- invilr(trans_f)
235 1
          if (spec[[box]]$sink) {
236 1
            parms[f_names] <- f[1:length(f)-1]
237
          } else {
238 1
            parms[f_names] <- f
239
          }
240
        }
241
      } else {
242 1
        parms[names(trans_f)] <- trans_f
243
      }
244
    }
245
  }
246

247
  # Transform parameters also for FOMC, DFOP, HS and logistic models
248 1
  for (pname in c("alpha", "beta", "k1", "k2", "tb", "kmax", "k0", "r")) {
249 1
    if (transform_rates) {
250 1
      pname_trans = paste0("log_", pname)
251 1
      if (!is.na(transparms[pname_trans])) {
252 1
        parms[pname] <- exp(transparms[pname_trans])
253
      }
254
    } else {
255 1
      if (!is.na(transparms[pname])) {
256 0
        parms[pname] <- transparms[pname]
257
      }
258
    }
259
  }
260

261
  # DFOP parameter g is now transformed using qlogis
262 1
  if (!is.na(transparms["g_qlogis"])) {
263 1
    g_qlogis <- transparms["g_qlogis"]
264 1
    parms["g"] <- plogis(g_qlogis)
265
  }
266
  # In earlier times we used ilr for g, so we keep this around
267 1
  if (!is.na(transparms["g_ilr"])) {
268 0
    g_ilr <- transparms["g_ilr"]
269 0
    parms["g"] <- invilr(g_ilr)[1]
270
  }
271 1
  if (!is.na(transparms["g"])) {
272 0
    parms["g"] <- transparms["g"]
273
  }
274

275 1
  return(parms)
276
}
277
# vim: set ts=2 sw=2 expandtab:

Read our documentation on viewing source code .

Loading