- 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

1 files changed, 37 insertions, 22 deletions

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)