From 6733555d7a9315c55001770bacc4c61c4d4f39d5 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Sun, 21 Jun 2015 01:46:51 +0200 Subject: Do the t-test for untransformed parameters --- R/mkinfit.R | 89 ++++++++++++++++++++++++++++++++++++++++++++++++------------- 1 file changed, 71 insertions(+), 18 deletions(-) (limited to 'R') 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 -- cgit v1.2.1