From 04d3370364a0a7881c4ae815c20fd37b74f5639a Mon Sep 17 00:00:00 2001
From: jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb>
Date: Tue, 5 Mar 2013 01:00:20 +0000
Subject: - Calculate confidence intervals for parameters based on the t
 distribution   Thanks to the CAKE developers at Tessella for the call to qt()
 - Calculate asymetric confidence intervals for backtransformed parameters -
 Overhaul and recompile vignettes

git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@74 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
---
 R/mkinfit.R | 59 +++++++++++++++++++++++++++++++++++++----------------------
 1 file changed, 37 insertions(+), 22 deletions(-)

(limited to 'R')

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 1061165..12a67bc 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -213,30 +213,42 @@ mkinfit <- function(mkinmod, observed,
   return(fit)
 }
 
-summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
+summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, ...) {
   param  <- object$par
   pnames <- names(param)
   p      <- length(param)
+  mod_vars <- names(object$mkinmod$diffs)
   covar  <- try(solve(0.5*object$hessian), silent = TRUE)   # unscaled covariance
-  if (!is.numeric(covar)) {
-    message <- "Cannot estimate covariance; system is singular"
-    warning(message)
-    covar <- matrix(data = NA, nrow = p, ncol = p)
-  } else message <- "ok"
-
-  rownames(covar) <- colnames(covar) <- pnames
   rdf    <- object$df.residual
   resvar <- object$ssr / rdf
-  se     <- sqrt(diag(covar) * resvar)
+  if (!is.numeric(covar)) {
+    covar <- NULL
+    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
+
+  }
+
   names(se) <- pnames
-  tval      <- param / se
   modVariance <- object$ssr / length(object$residuals)
 
-  param <- cbind(param, se)
-  dimnames(param) <- list(pnames, c("Estimate", "Std. Error"))
-
-  bparam <- as.matrix(object$bparms.optim)
-  dimnames(bparam) <- list(pnames, c("Estimate"))
+  param <- cbind(param, se, lci, uci)
+  dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))
+
+  blci <- buci <- numeric()
+  # Only use lower end of CI for one parameter at a time
+  for (pname in pnames) {
+    par.lower <- par.upper <- object$par
+    par.lower[pname] <- param[pname, "Lower"]
+    par.upper[pname] <- param[pname, "Upper"]
+    blci[pname] <- backtransform_odeparms(par.lower, mod_vars)[pname]
+    buci[pname] <- backtransform_odeparms(par.upper, mod_vars)[pname]
+  }
+  bparam <- cbind(object$bparms.optim, blci, buci)
+  dimnames(bparam) <- list(pnames, c("Estimate", "Lower", "Upper"))
 
   ans <- list(
     version = as.character(packageVersion("mkin")),
@@ -248,9 +260,11 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, ...) {
     residualVariance = resvar,
     sigma = sqrt(resvar),
     modVariance = modVariance,
-    df = c(p, rdf), cov.unscaled = covar,
+    df = c(p, rdf), 
+    cov.unscaled = covar,
     cov.scaled = covar * resvar,
-    info = object$info, niter = object$iterations,
+    info = object$info, 
+    niter = object$iterations,
     stopmess = message,
     par = param,
     bpar = bparam)
@@ -293,10 +307,10 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   else print(x$fixed)
   
   cat("\nOptimised, transformed parameters:\n")
-  printCoefmat(x$par, digits = digits, ...)
+  print(signif(x$par, digits = digits))
 
   cat("\nBacktransformed parameters:\n")
-  printCoefmat(x$bpar, digits = digits, ...)
+  print(signif(x$bpar, digits = digits))
 
   cat("\nResidual standard error:",
       format(signif(x$sigma, digits)), "on", rdf, "degrees of freedom\n")
@@ -323,12 +337,13 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
     print(x$SFORB, digits=digits,...)
   }
 
-  printcor <- is.numeric(x$cov.unscaled)
-  if (printcor){
+  cat("\nParameter correlation:\n")
+  if (!is.null(x$cov.unscaled)){
     Corr <- cov2cor(x$cov.unscaled)
     rownames(Corr) <- colnames(Corr) <- rownames(x$par)
-    cat("\nParameter correlation:\n")
     print(Corr, digits = digits, ...)
+  } else {
+    cat("Could not estimate covariance matrix; singular system:\n")
   }
 
   printdata <- !is.null(x$data)
-- 
cgit v1.2.1