Do the t-test for untransformed parametersv0.9-36

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