diff options
Diffstat (limited to 'R/mkinfit.R')
| -rw-r--r-- | R/mkinfit.R | 89 |
1 files changed, 71 insertions, 18 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index 7f9a55d9..cfd87829 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -234,10 +234,19 @@ mkinfit <- function(mkinmod, observed, }
}
+ # Define outtimes for model solution.
+ # Include time points at which observed data are available
+ # Make sure we include time 0, so initial values for state variables are for time 0
+ outtimes = sort(unique(c(observed$time, seq(min(observed$time),
+ max(observed$time),
+ length.out = n.outtimes))))
+
+
cost.old <- 1e100 # The first model cost should be smaller than this value
calls <- 0 # Counter for number of model solutions
out_predicted <- NA
- # Define the model cost function
+
+ # Define the model cost function for optimisation, including (back)transformations
cost <- function(P)
{
assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter
@@ -245,12 +254,6 @@ mkinfit <- function(mkinmod, observed, # Trace parameter values if requested
if(trace_parms) cat(P, "\n")
- # Time points at which observed data are available
- # Make sure we include time 0, so initial values for state variables are for time 0
- outtimes = sort(unique(c(observed$time, seq(min(observed$time),
- max(observed$time),
- length.out = n.outtimes))))
-
if(length(state.ini.optim) > 0) {
odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
@@ -313,6 +316,35 @@ mkinfit <- function(mkinmod, observed, return(mC)
}
+ # Define the model cost function for the t-test, without parameter transformation
+ cost_notrans <- function(P)
+ {
+ if(length(state.ini.optim) > 0) {
+ odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
+ names(odeini) <- c(state.ini.optim.boxnames, state.ini.fixed.boxnames)
+ } else {
+ odeini <- state.ini.fixed
+ names(odeini) <- state.ini.fixed.boxnames
+ }
+
+ odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
+
+ # Solve the system with current transformed parameter values
+ out <- mkinpredict(mkinmod, odeparms,
+ odeini, outtimes,
+ solution_type = solution_type,
+ use_compiled = use_compiled,
+ method.ode = method.ode,
+ atol = atol, rtol = rtol, ...)
+
+ mC <- modCost(out, observed, y = "value",
+ err = err, weight = weight, scaleVar = scaleVar)
+
+ return(mC)
+ }
+
+ # Define lower and upper bounds other than -Inf and Inf for parameters
+ # for which no internal transformation is requested in the call to mkinfit.
lower <- rep(-Inf, length(c(state.ini.optim, transparms.optim)))
upper <- rep(Inf, length(c(state.ini.optim, transparms.optim)))
names(lower) <- names(upper) <- c(names(state.ini.optim), names(transparms.optim))
@@ -431,6 +463,17 @@ mkinfit <- function(mkinmod, observed, bparms.fixed = c(state.ini.fixed, parms.fixed)
bparms.all = c(bparms.optim, parms.fixed)
+ # Attach the cost functions to the fit for post-hoc parameter uncertainty analysis
+ fit$cost <- cost
+ fit$cost_notrans <- cost_notrans
+
+ # Estimate the Hessian for the model cost without parameter transformations
+ # to make it possible to obtain the usual t-test
+ # Code ported from FME::modFit
+ Jac_notrans <- gradient(function(p, ...) cost_notrans(p)$residuals$res,
+ bparms.optim, centered = TRUE)
+ fit$hessian_notrans <- 2 * t(Jac_notrans) %*% Jac_notrans
+
# Collect initial parameter values in three dataframes
fit$start <- data.frame(value = c(state.ini.optim,
parms.optim))
@@ -489,29 +532,36 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, p <- length(param)
mod_vars <- names(object$mkinmod$diffs)
covar <- try(solve(0.5*object$hessian), silent = TRUE) # unscaled covariance
+ covar_notrans <- try(solve(0.5*object$hessian_notrans), silent = TRUE) # unscaled covariance
rdf <- object$df.residual
resvar <- object$ssr / rdf
if (!is.numeric(covar)) {
covar <- NULL
- se <- lci <- uci <- tval <- pval1 <- pval2 <- rep(NA, p)
+ se <- lci <- uci <- rep(NA, p)
} else {
rownames(covar) <- colnames(covar) <- pnames
se <- sqrt(diag(covar) * resvar)
lci <- param + qt(alpha/2, rdf) * se
uci <- param + qt(1-alpha/2, rdf) * se
- tval <- param/se
- pval1 <- 2 * pt(abs(tval), rdf, lower.tail = FALSE)
- pval2 <- pt(abs(tval), rdf, lower.tail = FALSE)
+ }
+
+ if (!is.numeric(covar_notrans)) {
+ covar_notrans <- NULL
+ se_notrans <- tval <- pval <- rep(NA, p)
+ } else {
+ rownames(covar_notrans) <- colnames(covar_notrans) <- bpnames
+ se_notrans <- sqrt(diag(covar_notrans) * resvar)
+ tval <- object$bparms.optim/se_notrans
+ pval <- pt(abs(tval), rdf, lower.tail = FALSE)
}
names(se) <- pnames
modVariance <- object$ssr / length(object$residuals)
- param <- cbind(param, se, lci, uci, tval, pval1, pval2)
- dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper",
- "t value", "Pr(>|t|)", "Pr(>t)"))
+ param <- cbind(param, se, lci, uci)
+ dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))
- bparam <- cbind(Estimate = object$bparms.optim, Lower = NA, Upper = NA)
+ bparam <- cbind(Estimate = object$bparms.optim, se_notrans, "t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA)
# Transform boundaries of CI for one parameter at a time,
# with the exception of sets of formation fractions (single fractions are OK).
@@ -617,7 +667,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . if(length(x$fixed$value) == 0) cat("None\n")
else print(x$fixed)
- cat("\nOptimised, transformed parameters:\n")
+ cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n")
print(signif(x$par, digits = digits))
if (x$calls > 0) {
@@ -634,8 +684,11 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("\nResidual standard error:",
format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
- cat("\nBacktransformed parameters:\n")
- print(signif(x$bpar, digits = digits))
+ cat("\nBacktransformed parameters:\n",
+ " Confidence intervals for internally transformed parameters are asymmetric.\n",
+ " t-test (unrealistically) based on the assumption of normal distribution\n",
+ " for estimators of untransformed parameters.\n", sep = "")
+ print(signif(x$bpar[, c(1, 3, 4, 5, 6)], digits = digits))
cat("\nChi2 error levels in percent:\n")
x$errmin$err.min <- 100 * x$errmin$err.min
|