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 .