From fff1fc581da5b4ff935ebd4d7ded02f750472fdc Mon Sep 17 00:00:00 2001 From: jranke Date: Tue, 27 Mar 2012 01:03:18 +0000 Subject: Start of the transition to fitting transformed parameters. Many things are broken now (see TODO list) git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@20 edb9625f-4e0d-4859-8d74-9fd3b1da38cb --- R/ilr.R | 44 +++++++++++++++ R/mkinfit.R | 96 ++++++++++++++++++--------------- R/mkinmod.R | 168 +++++++++++++++++++++------------------------------------- R/mkinstart.R | 23 -------- 4 files changed, 158 insertions(+), 173 deletions(-) create mode 100644 R/ilr.R delete mode 100644 R/mkinstart.R (limited to 'R') diff --git a/R/ilr.R b/R/ilr.R new file mode 100644 index 0000000..f2a6f1b --- /dev/null +++ b/R/ilr.R @@ -0,0 +1,44 @@ +# $Id: ilr.R 120 2011-09-02 14:25:35Z jranke $ + +# Copyright (C) 2012 René Lehmann, Johannes Ranke +# Contact: jranke@uni-bremen.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see + +ilr <- function(x) { + z <- vector() + for (i in 1:(length(x) - 1)) { + z[i] <- sqrt(i/(i+1)) * log((prod(x[1:i]))^(1/i) / x[i+1]) + } + return(z) +} + +invilr<-function(x) { + D <- length(x) + 1 + z <- c(x, 0) + y <- rep(0, D) + s <- sqrt(1:D*2:(D+1)) + q <- z/s + y[1] <- sum(q[1:D]) + for (i in 2:D) { + y[i] <- sum(q[i:D]) - sqrt((i-1)/i) * z[i-1] + } + z <- vector() + for (i in 1:D) { + z[i] <- exp(y[i])/sum(exp(y)) + } + return(z) +} diff --git a/R/mkinfit.R b/R/mkinfit.R index b742742..83896bc 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -1,7 +1,7 @@ # $Id: mkinfit.R 120 2011-09-02 14:25:35Z jranke $ -# Copyright (C) 2010-2011 Johannes Ranke -# Contact: mkin-devel@lists.berlios.de +# Copyright (C) 2010-2012 Johannes Ranke +# Contact: jranke@uni-bremen.de # The summary function is an adapted and extended version of summary.modFit # from the FME package, v 1.1 by Soetart and Petzoldt, which was in turn # inspired by summary.nls.lm @@ -22,9 +22,8 @@ # this program. If not, see mkinfit <- function(mkinmod, observed, - parms.ini = rep(0.1, length(mkinmod$parms)), + parms.ini = "auto", state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), - lower = 0, upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1], eigen = FALSE, @@ -34,14 +33,24 @@ mkinfit <- function(mkinmod, observed, ...) { mod_vars <- names(mkinmod$diffs) - # Subset dataframe with mapped (modelled) variables + # Subset dataframe of observed variables with mapped (modelled) variables observed <- subset(observed, name %in% names(mkinmod$map)) # Get names of observed variables obs_vars = unique(as.character(observed$name)) - # Name the parameters if they are not named yet - if(is.null(names(parms.ini))) names(parms.ini) <- mkinmod$parms - + # Define initial values for parameters where not specified by the user + if (parms.ini[[1]] == "auto") parms.ini = vector() + defaultpar.names <- setdiff(mkinmod$parms, names(parms.ini)) + for (parmname in defaultpar.names) { + # Default values for the FOMC model + if (parmname == "alpha") parms.ini[parmname] = 1 + if (parmname == "beta") parms.ini[parmname] = 10 + # Default values for log of rate constants, depending on the parameterisation + if (substr(parmname, 1, 2) == "k_") parms.ini[parmname] = -2.3 + # Default values for ilr of formation fractions + if (substr(parmname, 1, 2) == "f_") parms.ini[parmname] = 0.1 + } + # Name the inital parameter values if they are not named yet if(is.null(names(state.ini))) names(state.ini) <- mod_vars @@ -58,14 +67,6 @@ mkinfit <- function(mkinmod, observed, names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_") } - # Set upper limit for formation fractions to one if formation fractions are - # directly defined and if no user input for upper limit is given - if (all(upper==Inf) & any(grepl("f_", names(parms.ini)))==TRUE){ - upper=c( rep(Inf,length(parms.optim))) - upper[grep("f_", names(parms.optim))]=1 - upper=c(rep(Inf, length(state.ini.optim)), upper) - } - # Decide if the solution of the model can be based on a simple analytical # formula, the spectral decomposition of the matrix (fundamental system) # or a numeric ode solver from the deSolve package @@ -80,12 +81,25 @@ mkinfit <- function(mkinmod, observed, # if necessary if(solution == "deSolve") { mkindiff <- function(t, state, parms) { + transparms <- vector() + transparms <- c(transparms, exp(parms[grep("^k_", names(parms))])) + for (box in mod_vars) { + path_indices <- grep(paste("^f", box, sep = "_"), names(parms)) + n_paths <- length(path_indices) + if (n_paths > 0) { + f <- invilr(parms[grep(paste("^f", box, sep="_"), names(parms))]) + transparms <- c(transparms, f[1:n_paths]) + } + } + otherparms <- parms[setdiff(names(parms), names(transparms))] + parms <- c(transparms, otherparms) + time <- t diffs <- vector() for (box in mod_vars) { diffname <- paste("d", box, sep="_") - diffs[diffname] <- with(as.list(c(time,state, parms)), + diffs[diffname] <- with(as.list(c(time, state, parms)), eval(parse(text=mkinmod$diffs[[box]]))) } return(list(c(diffs))) @@ -251,7 +265,7 @@ mkinfit <- function(mkinmod, observed, } return(mC) } - fit <- modFit(cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, ...) + fit <- modFit(cost, c(state.ini.optim, parms.optim), ...) # We need to return some more data for summary and plotting fit$solution <- solution @@ -323,7 +337,7 @@ mkinfit <- function(mkinmod, observed, DT90 = log(10)/k_tot for (k_name in k_names) { - fit$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot + fit$ff[[sub("^k_", "", k_name)]] = parms.all[[k_name]] / k_tot } } if (type == "FOMC") { @@ -395,7 +409,7 @@ mkinfit <- function(mkinmod, observed, if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o for (k_out_name in k_out_names) { - fit$ff[[sub("k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output + fit$ff[[sub("^k_", "", k_out_name)]] = parms.all[[k_out_name]] / k_1output } # Return the eigenvalues for comparison with DFOP rate constants fit$SFORB[[paste(obs_var, "b1", sep="_")]] = b1 @@ -416,8 +430,11 @@ mkinfit <- function(mkinmod, observed, return(fit) } -summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov = FALSE,...) { +summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, ...) { param <- object$par + if (object$logk) { + names(param) <- sub("^k_", "log k_", names(param)) + } pnames <- names(param) p <- length(param) covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance @@ -441,32 +458,25 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ff = TRUE, cov warning(message) } else message <- "ok" - param <- cbind(param, se, tval, pt(tval, rdf, lower.tail = FALSE)) - dimnames(param) <- list(pnames, c("Estimate", "Std. Error", - "t value", "Pr(>t)")) - if(cov) - ans <- list(residuals = object$residuals, - residualVariance = resvar, - sigma = sqrt(resvar), - modVariance = modVariance, - df = c(p, rdf), cov.unscaled = covar, - cov.scaled = covar * resvar, - info = object$info, niter = object$iterations, - stopmess = message, - par = param) - else - ans <- list(residuals = object$residuals, - residualVariance = resvar, - sigma = sqrt(resvar), - modVariance = modVariance, - df = c(p, rdf), - info = object$info, niter = object$iterations, - stopmess = message, - par = param) + param <- cbind(param, se) + dimnames(param) <- list(pnames, c("Estimate", "Std. Error")) + + ans <- list(residuals = object$residuals, + residualVariance = resvar, + sigma = sqrt(resvar), + modVariance = modVariance, + df = c(p, rdf), cov.unscaled = covar, + cov.scaled = covar * resvar, + info = object$info, niter = object$iterations, + stopmess = message, + par = param) ans$diffs <- object$diffs if(data) ans$data <- object$data ans$start <- object$start + if (object$logk) { + names(ans$start) <- sub("^k_", "log k_", names(ans$start)) + } ans$fixed <- object$fixed ans$errmin <- object$errmin if(distimes) ans$distimes <- object$distimes diff --git a/R/mkinmod.R b/R/mkinmod.R index 31c778c..b49c813 100644 --- a/R/mkinmod.R +++ b/R/mkinmod.R @@ -1,7 +1,7 @@ # $Id: mkinmod.R 71 2010-09-12 01:13:36Z jranke $ -# Copyright (C) 2010 Johannes Ranke -# Contact: mkin-devel@lists.berlios.de +# Copyright (C) 2010-2012 Johannes Ranke +# Contact: jranke@uni-bremen.de # This file is part of the R package mkin @@ -33,21 +33,27 @@ mkinmod <- function(...) if(spec[[1]]$type %in% c("FOMC", "DFOP", "HS")) { mat = FALSE if(!is.null(spec[[1]]$to)) { - message <- paste("Only constant formation fractions over time are implemented.", - "Depending on the reason for the time dependence of degradation, this may be unrealistic", + message <- paste( + "Only constant formation fractions over time are implemented.", + "Depending on the reason for the time dependence of degradation", + "this may be unrealistic", sep="\n") warning(message) } else message <- "ok" } else mat = TRUE # Establish list of differential equations + # as well as map from observed compartments to differential equations for (varname in obs_vars) { if(is.null(spec[[varname]]$type)) stop( "Every argument to mkinmod must be a list containing a type component") if(!spec[[varname]]$type %in% c("SFO", "FOMC", "DFOP", "HS", "SFORB")) stop( "Available types are SFO, FOMC, DFOP, HS and SFORB only") - new_parms <- vector() + if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS") & match(varname, obs_vars) != 1) { + stop(paste("Types FOMC, DFOP and HS are only implemented for the first compartment,", + "which is assumed to be the source compartment")) + } # New (sub)compartments (boxes) needed for the model type new_boxes <- switch(spec[[varname]]$type, @@ -62,124 +68,75 @@ mkinmod <- function(...) # Start a new differential equation for each new box new_diffs <- paste("d_", new_boxes, " =", sep="") + names(new_diffs) <- new_boxes + diffs <- c(diffs, new_diffs) + } + + # Create content of differential equations and build parameter list + for (varname in obs_vars) + { + # Add first-order term to first (or only) box for SFO and SFORB + box_1 = map[[varname]][[1]] # This is the only box unless type is SFORB + if(spec[[varname]]$type %in% c("SFO", "SFORB")) { + k_compound <- paste("k", box_1, sep="_") + origin_term <- paste(k_compound, "*", box_1) + parms <- c(parms, k_compound) + } - # Turn on sink if not specified otherwise - if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE # Construct and add FOMC term and add FOMC parameters if needed if(spec[[varname]]$type == "FOMC") { - if(match(varname, obs_vars) != 1) { - stop("Type FOMC is only possible for the first compartment, which is assumed to be the source compartment") - } - if(spec[[varname]]$sink == FALSE) { - stop("Turning off the sink for the FOMC model is not implemented") - } # From p. 53 of the FOCUS kinetics report - nonlinear_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]]) - new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) - new_parms <- c("alpha", "beta") - ff <- vector() + origin_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", box_1) + parms <- c(parms, "alpha", "beta") } # Construct and add DFOP term and add DFOP parameters if needed if(spec[[varname]]$type == "DFOP") { - if(match(varname, obs_vars) != 1) { - stop("Type DFOP is only possible for the first compartment, which is assumed to be the source compartment") - } - if(spec[[varname]]$sink == FALSE) { - stop("Turning off the sink for the DFOP model is not implemented") - } # From p. 57 of the FOCUS kinetics report - nonlinear_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", new_boxes[[1]]) - new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) - new_parms <- c("k1", "k2", "g") - ff <- vector() + origin_term <- paste("((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) *", box_1) + parms <- c(parms, "k1", "k2", "g") } # Construct and add HS term and add HS parameters if needed if(spec[[varname]]$type == "HS") { - if(match(varname, obs_vars) != 1) { - stop("Type HS is only possible for the first compartment, which is assumed to be the source compartment") - } - if(spec[[varname]]$sink == FALSE) { - stop("Turning off the sink for the HS model is not implemented") - } # From p. 55 of the FOCUS kinetics report - nonlinear_term <- paste("ifelse(time <= tb, k1, k2)", "*", new_boxes[[1]]) - new_diffs[[1]] <- paste(new_diffs[[1]], "-", nonlinear_term) - new_parms <- c("k1", "k2", "tb") - ff <- vector() + origin_term <- paste("ifelse(time <= tb, k1, k2)", "*", box_1) + parms <- c(parms, "k1", "k2", "tb") } - # Construct terms for transfer to sink and add if appropriate - - if(spec[[varname]]$sink) { - # Add first-order sink term to first (or only) box for SFO and SFORB - if(spec[[varname]]$type %in% c("SFO", "SFORB")) { - k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_") - sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]]) - new_diffs[[1]] <- paste(new_diffs[[1]], sink_term) - new_parms <- k_compound_sink - } - } + # Add origin decline term to box 1 (usually the only box, unless type is SFORB) + diffs[[box_1]] <- paste(diffs[[box_1]], "-", origin_term) - # Add reversible binding if appropriate + # Add reversible binding for SFORB models if(spec[[varname]]$type == "SFORB") { - k_free_bound <- paste("k", varname, "free", "bound", sep="_") - k_bound_free <- paste("k", varname, "bound", "free", sep="_") - reversible_binding_terms <- c( - paste("-", k_free_bound, "*", new_boxes[[1]], "+", k_bound_free, "*", new_boxes[[2]]), - paste("+", k_free_bound, "*", new_boxes[[1]], "-", k_bound_free, "*", new_boxes[[2]])) - new_diffs <- paste(new_diffs, reversible_binding_terms) - new_parms <- c(new_parms, k_free_bound, k_bound_free) + box_2 = map[[varname]][[2]] + k_free_bound <- paste("k", varname, "free", "bound", sep="_") + k_bound_free <- paste("k", varname, "bound", "free", sep="_") + reversible_binding_term_1 <- paste("-", k_free_bound, "*", box_1, "+", + k_bound_free, "*", box_2) + diffs[[box_1]] <- paste(diffs[[box_1]], reversible_binding_term_1) + reversible_binding_term_2 <- paste("+", k_free_bound, "*", box_1, "-", + k_bound_free, "*", box_2) + diffs[[box_2]] <- paste(diffs[[box_2]], reversible_binding_term_2) + parms <- c(parms, k_free_bound, k_bound_free) } - # Add observed variable to model - parms <- c(parms, new_parms) - names(new_diffs) <- new_boxes - diffs <- c(diffs, new_diffs) - } - # Transfer between compartments - for (varname in obs_vars) { + # Transfer between compartments to <- spec[[varname]]$to if(!is.null(to)) { - origin_box <- switch(spec[[varname]]$type, - SFO = varname, - FOMC = varname, - DFOP = varname, - HS = varname, - SFORB = paste(varname, "free", sep="_")) - fraction_left <- NULL + # Name of box from which transfer takes place + origin_box <- box_1 + + # Add transfer terms to listed compartments for (target in to) { target_box <- switch(spec[[target]]$type, SFO = target, SFORB = paste(target, "free", sep="_")) - if(spec[[varname]]$type %in% c("SFO", "SFORB")) { - k_from_to <- paste("k", origin_box, target_box, sep="_") - diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", - k_from_to, "*", origin_box) - diffs[[target_box]] <- paste(diffs[[target_box]], "+", - k_from_to, "*", origin_box) - parms <- c(parms, k_from_to) - } - if(spec[[varname]]$type %in% c("FOMC", "DFOP", "HS")) { - fraction_to_target = paste("f_to", target, sep="_") - fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", - sep="") - if(is.null(fraction_left)) { - fraction_really_to_target = fraction_to_target - fraction_left = fraction_not_to_target - } else { - fraction_really_to_target = paste(fraction_left, " * ", - fraction_to_target, sep="") - fraction_left = paste(fraction_left, " * ", - fraction_not_to_target, sep="") - } - ff[target_box] = fraction_really_to_target - diffs[[target_box]] <- paste(diffs[[target_box]], "+", - ff[target_box], "*", nonlinear_term) - parms <- c(parms, fraction_to_target) - } + fraction_to_target = paste("f", origin_box, "to", target, sep="_") + diffs[[target_box]] <- paste(diffs[[target_box]], "+", + fraction_to_target, "*", origin_term) + parms <- c(parms, fraction_to_target) } } } @@ -193,26 +150,23 @@ mkinmod <- function(...) for (from in boxes) { for (to in boxes) { if (from == to) { - k.candidate = paste("k", from, c(boxes, "sink"), sep="_") - k.candidate = sub("free.*bound", "free_bound", k.candidate) - k.candidate = sub("bound.*free", "bound_free", k.candidate) - k.effective = intersect(model$parms, k.candidate) - m[from,to] = ifelse(length(k.effective) > 0, - paste("-", k.effective, collapse = " "), "0") + k.candidate = paste("k", from, sep="_") + m[from,to] = ifelse(k.candidate %in% model$parms, + paste("-", k.candidate), "0") } else { + f.candidate = paste("f", from, "to", to, sep="_") k.candidate = paste("k", from, to, sep="_") k.candidate = sub("free.*bound", "free_bound", k.candidate) k.candidate = sub("bound.*free", "bound_free", k.candidate) - k.effective = intersect(model$parms, k.candidate) - m[to, from] = ifelse(length(k.effective) > 0, - k.effective, "0") + m[to, from] = ifelse(f.candidate %in% model$parms, + paste(f.candidate, " * k_", from, sep=""), + ifelse(k.candidate %in% model$parms, k.candidate, "0")) } } } model$coefmat <- m } - if (exists("ff")) model$ff = ff class(model) <- "mkinmod" - invisible(model) + return(model) } diff --git a/R/mkinstart.R b/R/mkinstart.R deleted file mode 100644 index 46e5bed..0000000 --- a/R/mkinstart.R +++ /dev/null @@ -1,23 +0,0 @@ -mkinstart <- function(model, data, mode = "auto") -{ - if (class(model) != "mkinmod") stop("The first argument must be a model of class mkinmod") - names <- model$parms - observed <- names(model$map) - if(!all(observed %in% levels(data$name))) stop("The data must contain the observed variables used in the model") - for (obs in observed) - { - tmp <- subset(data, name == obs) - max <- tmp[which.max(tmp$value), ] - type = names(model$map[[obs]])[[1]] - kinmodel <- ifelse(type == "SFORB", "DFOP", type) - tmp.longdata <- subset(data, name == obs & time >= max$time) - tmp.widedata <- mkin_long_to_wide(tmp.longdata, outtime = "t") - names(tmp.widedata) <- c("t", "parent") - tmp.fit <- kinfit( - kindata = tmp.widedata, - kinmodels = kinmodel, - parent.0.user = max$value) - if(class(tmp.fit[[kinmodel]]) == "try-error") stop(paste("Automatic generation of starting parameters failed\nkinfit failed to find a", kinmodel, "fit for", obs)) - tmp.results <- kinresults(tmp.fit) - } -} -- cgit v1.2.1