diff options
| -rw-r--r-- | R/endpoints.R | 94 | ||||
| -rw-r--r-- | R/mkinerrmin.R | 2 | ||||
| -rw-r--r-- | R/mkinfit.R | 98 |
3 files changed, 100 insertions, 94 deletions
diff --git a/R/endpoints.R b/R/endpoints.R new file mode 100644 index 0000000..98e290d --- /dev/null +++ b/R/endpoints.R @@ -0,0 +1,94 @@ +endpoints <- function(fit, pseudoDT50 = FALSE) {
+ # Calculate dissipation times DT50 and DT90 and, if necessary, formation
+ # fractions and SFORB eigenvalues from optimised parameters
+ ep <- list()
+ obs_vars <- fit$obs_vars
+ parms.all <- fit$parms.all
+ ep$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
+ DT90 = rep(NA, length(obs_vars)),
+ row.names = obs_vars)
+ ep$ff <- vector()
+ ep$SFORB <- vector()
+ for (obs_var in obs_vars) {
+ type = names(fit$mkinmod$map[[obs_var]])[1]
+ if (type == "SFO") {
+ k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
+ k_tot = sum(parms.all[k_names])
+ DT50 = log(2)/k_tot
+ DT90 = log(10)/k_tot
+ for (k_name in k_names)
+ {
+ ep$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
+ }
+ }
+ if (type == "FOMC") {
+ alpha = parms.all["alpha"]
+ beta = parms.all["beta"]
+ DT50 = beta * (2^(1/alpha) - 1)
+ DT90 = beta * (10^(1/alpha) - 1)
+ }
+ if (type == "DFOP") {
+ k1 = parms.all["k1"]
+ k2 = parms.all["k2"]
+ g = parms.all["g"]
+ f <- function(t, x) {
+ fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
+ (fraction - (1 - x/100))^2
+ }
+ DTmax <- 1000
+ DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum
+ DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o)
+ DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum
+ DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o)
+ }
+ if (type == "HS") {
+ k1 = parms.all["k1"]
+ k2 = parms.all["k2"]
+ tb = parms.all["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)
+ }
+ if (type == "SFORB") {
+ # FOCUS kinetics (2006), p. 60 f
+ k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
+ k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
+ k_1output = sum(parms.all[k_out_names])
+ k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]
+ k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]
+
+ sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + 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
+
+ SFORB_fraction = function(t) {
+ ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
+ ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
+ }
+ f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2
+ max_DT <- 1000
+ DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum
+ if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o
+ f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2
+ DT90.o <- optimize(f_90, c(0.01, max_DT))$minimum
+ if (abs(DT90.o - max_DT) < 0.01) DT90 = NA else DT90 = DT90.o
+
+ for (k_out_name in k_out_names)
+ {
+ ep$ff[[sub("k_", "", k_out_name)]] = parms.all[[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
+ }
+ ep$distimes[obs_var, ] = c(DT50, DT90)
+ }
+ return(ep)
+}
diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R index 5415b41..c278647 100644 --- a/R/mkinerrmin.R +++ b/R/mkinerrmin.R @@ -76,5 +76,5 @@ mkinerrmin <- function(fit, alpha = 0.05) errmin[obs_var, c("err.min", "n.optim", "df")] <- errmin.tmp
}
- return(errmin)
+ return(errmin)
}
diff --git a/R/mkinfit.R b/R/mkinfit.R index b2641e6..d3ee5e7 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -189,97 +189,8 @@ mkinfit <- function(mkinmod, observed, value = c(state.ini.fixed, parms.fixed))
fit$fixed$type = c(rep("state", length(state.ini.fixed)), rep("deparm", length(parms.fixed)))
-
parms.all = backtransform_odeparms(c(fit$par, parms.fixed), mod_vars)
- # Calculate dissipation times DT50 and DT90 and, if necessary, formation
- # fractions and SFORB eigenvalues from optimised parameters
- fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)), DT90 = rep(NA, length(obs_vars)),
- row.names = obs_vars)
- fit$ff <- vector()
- fit$SFORB <- vector()
- for (obs_var in obs_vars) {
- type = names(mkinmod$map[[obs_var]])[1]
- if (type == "SFO") {
- k_names = grep(paste("k", obs_var, sep="_"), names(parms.all), value=TRUE)
- k_tot = sum(parms.all[k_names])
- DT50 = log(2)/k_tot
- DT90 = log(10)/k_tot
- for (k_name in k_names)
- {
- fit$ff[[sub("k_", "", k_name)]] = parms.all[[k_name]] / k_tot
- }
- }
- if (type == "FOMC") {
- alpha = parms.all["alpha"]
- beta = parms.all["beta"]
- DT50 = beta * (2^(1/alpha) - 1)
- DT90 = beta * (10^(1/alpha) - 1)
- }
- if (type == "DFOP") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- g = parms.all["g"]
- f <- function(t, x) {
- fraction <- g * exp( - k1 * t) + (1 - g) * exp( - k2 * t)
- (fraction - (1 - x/100))^2
- }
- DTmax <- 1000
- DT50.o <- optimize(f, c(0.001, DTmax), x=50)$minimum
- DT50 = ifelse(DTmax - DT50.o < 0.1, NA, DT50.o)
- DT90.o <- optimize(f, c(0.001, DTmax), x=90)$minimum
- DT90 = ifelse(DTmax - DT90.o < 0.1, NA, DT90.o)
- }
- if (type == "HS") {
- k1 = parms.all["k1"]
- k2 = parms.all["k2"]
- tb = parms.all["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)
- }
- if (type == "SFORB") {
- # FOCUS kinetics (2006), p. 60 f
- k_out_names = grep(paste("k", obs_var, "free", sep="_"), names(parms.all), value=TRUE)
- k_out_names = setdiff(k_out_names, paste("k", obs_var, "free", "bound", sep="_"))
- k_1output = sum(parms.all[k_out_names])
- k_12 = parms.all[paste("k", obs_var, "free", "bound", sep="_")]
- k_21 = parms.all[paste("k", obs_var, "bound", "free", sep="_")]
-
- sqrt_exp = sqrt(1/4 * (k_12 + k_21 + k_1output)^2 + k_12 * k_21 - (k_12 + 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
-
- SFORB_fraction = function(t) {
- ((k_12 + k_21 - b1)/(b2 - b1)) * exp(-b1 * t) +
- ((k_12 + k_21 - b2)/(b1 - b2)) * exp(-b2 * t)
- }
- f_50 <- function(t) (SFORB_fraction(t) - 0.5)^2
- max_DT <- 1000
- DT50.o <- optimize(f_50, c(0.01, max_DT))$minimum
- if (abs(DT50.o - max_DT) < 0.01) DT50 = NA else DT50 = DT50.o
- f_90 <- function(t) (SFORB_fraction(t) - 0.1)^2
- DT90.o <- optimize(f_90, c(0.01, max_DT))$minimum
- 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
- }
-
- # Return the eigenvalues for comparison with DFOP rate constants
- fit$SFORB[[paste(obs_var, "b1", sep="_")]] = b1
- fit$SFORB[[paste(obs_var, "b2", sep="_")]] = b2
- }
- fit$distimes[obs_var, ] = c(DT50, DT90)
- }
-
# Collect observed, predicted and residuals
data <- merge(fit$observed, fit$predicted, by = c("time", "name"))
names(data) <- c("time", "variable", "observed", "predicted")
@@ -342,10 +253,11 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) { ans$errmin <- mkinerrmin(object, alpha = 0.05)
ans$parms.all <- object$parms.all
- if (!is.null(object$ff))
- ans$ff <- object$ff
- if(distimes) ans$distimes <- object$distimes
- if(length(object$SFORB) != 0) ans$SFORB <- object$SFORB
+ ep <- endpoints(object)
+ if (!is.null(ep$ff))
+ ans$ff <- ep$ff
+ if(distimes) ans$distimes <- ep$distimes
+ if(length(ep$SFORB) != 0) ans$SFORB <- ep$SFORB
class(ans) <- c("summary.mkinfit", "summary.modFit")
return(ans)
}
|