From 7624a2b8398b4ad665a3b0b622488e1893a5ee7c Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 21 Oct 2019 12:11:34 +0200
Subject: Refactor mkinfit, infrastructure work

mkinfit objects now include an ll() function to calculate the
log-likelihood. Part of the code was refactored, hopefully making it
easier to read and maintain. IRLS is currently the default algorithm for
the error model "obs", for no particular reason. This may be subject
to change when I get around to investigate.

Slow tests are now in a separate subdirectory and will probably
only be run by my own Makefile target.

Formatting of test logs is improved.

Roundtripping error model parameters works with a precision of 10% when
we use lots of replicates in the synthetic data (see slow tests). This
is not new in this commit, but as I think it is reasonable this
closes #7.
---
 R/mkinerrmin.R |   7 +-
 R/mkinfit.R    | 328 ++++++++++++++++++++++++++++++---------------------------
 2 files changed, 178 insertions(+), 157 deletions(-)

(limited to 'R')

diff --git a/R/mkinerrmin.R b/R/mkinerrmin.R
index c0c6fad7..ce4877d2 100644
--- a/R/mkinerrmin.R
+++ b/R/mkinerrmin.R
@@ -1,4 +1,4 @@
-# Copyright (C) 2010-2014 Johannes Ranke
+# Copyright (C) 2010-2019 Johannes Ranke
 # Contact: jranke@uni-bremen.de
 
 # This file is part of the R package mkin
@@ -62,9 +62,8 @@ mkinerrmin <- function(fit, alpha = 0.05)
     n.k.optim <- n.k.optim + length(grep(paste("^log_k", obs_var, sep="_"),
                                          names(parms.optim)))
     n.k__iore.optim <- length(grep(paste("^k__iore", obs_var, sep="_"), names(parms.optim)))
-    n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore", obs_var,
-							 sep = "_"),
-						   names(parms.optim)))
+    n.k__iore.optim <- n.k__iore.optim + length(grep(paste("^log_k__iore",
+          obs_var, sep = "_"), names(parms.optim)))
 
     n.N.optim <- length(grep(paste("^N", obs_var, sep="_"), names(parms.optim)))
 
diff --git a/R/mkinfit.R b/R/mkinfit.R
index b5e69e67..7e2b8cac 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -34,7 +34,7 @@ mkinfit <- function(mkinmod, observed,
   quiet = FALSE,
   atol = 1e-8, rtol = 1e-10, n.outtimes = 100,
   error_model = c("const", "obs", "tc"),
-  error_model_algorithm = c("d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", "OLS"),
+  error_model_algorithm = c("auto", "d_3", "direct", "twostep", "threestep", "fourstep", "IRLS", "OLS"),
   reweight.tol = 1e-8, reweight.max.iter = 10,
   trace_parms = FALSE,
   ...)
@@ -244,13 +244,21 @@ mkinfit <- function(mkinmod, observed,
     }
   }
 
-  # Get the error model
+  # Get the error model and the algorithm for fitting
   err_mod <- match.arg(error_model)
   error_model_algorithm = match.arg(error_model_algorithm)
+  if (error_model_algorithm == "OLS") {
+    if (err_mod != "const") stop("OLS is only appropriate for constant variance")
+  }
+  if (error_model_algorithm == "auto") {
+    error_model_algorithm = switch(err_mod,
+      const = "OLS", obs = "IRLS", tc = "d_3")
+  }
   errparm_names <- switch(err_mod,
     "const" = "sigma",
     "obs" = paste0("sigma_", obs_vars),
     "tc" = c("sigma_low", "rsd_high"))
+  errparm_names_optim <- if (error_model_algorithm == "OLS") NULL else errparm_names
 
   # Define starting values for the error model
   if (err.ini[1] != "auto") {
@@ -271,6 +279,11 @@ mkinfit <- function(mkinmod, observed,
     }
     names(errparms) <- errparm_names
   }
+  if (error_model_algorithm == "OLS") {
+    errparms_optim <- NULL
+  } else {
+    errparms_optim <- errparms
+  }
 
   # Define outtimes for model solution.
   # Include time points at which observed data are available
@@ -278,49 +291,51 @@ mkinfit <- function(mkinmod, observed,
                                               max(observed$time),
                                               length.out = n.outtimes))))
 
-  # Define log-likelihood function for optimisation, including (back)transformations
-  nlogLik <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...)
+  # Define the objective function for optimisation, including (back)transformations
+  cost_function <- function(P, trans = TRUE, OLS = FALSE, fixed_degparms = FALSE, fixed_errparms = FALSE, update_data = TRUE, ...)
   {
     assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter
 
     # Trace parameter values if requested and if we are actually optimising
     if(trace_parms & update_data) cat(P, "\n")
 
+    # Determine local parameter values for the cost estimation
     if (is.numeric(fixed_degparms)) {
-      degparms <- fixed_degparms
-      errparms <- P # This version of errparms is local to the function
+      cost_degparms <- fixed_degparms
+      cost_errparms <- P
       degparms_fixed = TRUE
     } else {
       degparms_fixed = FALSE
     }
 
     if (is.numeric(fixed_errparms)) {
-      degparms <- P
-      errparms <- fixed_errparms # Local to the function
+      cost_degparms <- P
+      cost_errparms <- fixed_errparms
       errparms_fixed = TRUE
     } else {
       errparms_fixed = FALSE
     }
 
     if (OLS) {
-      degparms <- P
+      cost_degparms <- P
+      cost_errparms <- numeric(0)
     }
 
     if (!OLS & !degparms_fixed & !errparms_fixed) {
-      degparms <- P[1:(length(P) - length(errparms))]
-      errparms <- P[(length(degparms) + 1):length(P)]
+      cost_degparms <- P[1:(length(P) - length(errparms))]
+      cost_errparms <- P[(length(cost_degparms) + 1):length(P)]
     }
 
     # Initial states for t0
     if(length(state.ini.optim) > 0) {
-      odeini <- c(degparms[1:length(state.ini.optim)], state.ini.fixed)
+      odeini <- c(cost_degparms[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.optim <- degparms[(length(state.ini.optim) + 1):length(degparms)]
+    odeparms.optim <- cost_degparms[(length(state.ini.optim) + 1):length(cost_degparms)]
 
     if (trans == TRUE) {
       odeparms <- c(odeparms.optim, transparms.fixed)
@@ -342,53 +357,55 @@ mkinfit <- function(mkinmod, observed,
     out_long <- mkin_wide_to_long(out, time = "time")
 
     if (err_mod == "const") {
-      observed$std <- errparms["sigma"]
+      observed$std <- if (OLS) NA else cost_errparms["sigma"]
     }
     if (err_mod == "obs") {
       std_names <- paste0("sigma_", observed$name)
-      observed$std <- errparms[std_names]
+      observed$std <- cost_errparms[std_names]
     }
     if (err_mod == "tc") {
       tmp <- merge(observed, out_long, by = c("time", "name"))
       tmp$name <- ordered(tmp$name, levels = obs_vars)
       tmp <- tmp[order(tmp$name, tmp$time), ]
-      observed$std <- sqrt(errparms["sigma_low"]^2 + tmp$value.y^2 * errparms["rsd_high"]^2)
+      observed$std <- sqrt(cost_errparms["sigma_low"]^2 + tmp$value.y^2 * cost_errparms["rsd_high"]^2)
     }
 
-    data_log_lik <- merge(observed[c("name", "time", "value", "std")], out_long,
+    cost_data <- merge(observed[c("name", "time", "value", "std")], out_long,
                          by = c("name", "time"), suffixes = c(".observed", ".predicted"))
 
     if (OLS) {
-      nlogLik <- with(data_log_lik, sum((value.observed - value.predicted)^2))
+      # Cost is the sum of squared residuals
+      cost <- with(cost_data, sum((value.observed - value.predicted)^2))
     } else {
-      nlogLik <- - with(data_log_lik,
+      # Cost is the negative log-likelihood
+      cost <- - with(cost_data,
         sum(dnorm(x = value.observed, mean = value.predicted, sd = std, log = TRUE)))
     }
 
-    # We update the current likelihood and data during the optimisation, not
+    # We update the current cost and data during the optimisation, not
     # during hessian calculations
     if (update_data) {
 
       assign("out_predicted", out_long, inherits = TRUE)
-      assign("data_errmod", data_log_lik, inherits = TRUE)
+      assign("current_data", cost_data, inherits = TRUE)
 
-      if (nlogLik < nlogLik.current) {
-        assign("nlogLik.current", nlogLik, inherits = TRUE)
+      if (cost < cost.current) {
+        assign("cost.current", cost, inherits = TRUE)
         if (!quiet) cat(ifelse(OLS, "Sum of squared residuals", "Negative log-likelihood"),
-                        " at call ", calls, ": ", nlogLik.current, "\n", sep = "")
+                        " at call ", calls, ": ", cost.current, "\n", sep = "")
       }
     }
-    return(nlogLik)
+    return(cost)
   }
 
-  n_optim <- length(c(state.ini.optim, transparms.optim, errparm_names))
   names_optim <- c(names(state.ini.optim),
                    names(transparms.optim),
-                   errparm_names)
+                   errparm_names_optim)
+  n_optim <- length(names_optim)
 
   # Define lower and upper bounds other than -Inf and Inf for parameters
   # for which no internal transformation is requested in the call to mkinfit
-  # and for error model parameters
+  # and for optimised error model parameters
   lower <- rep(-Inf, n_optim)
   upper <- rep(Inf, n_optim)
   names(lower) <- names(upper) <- names_optim
@@ -416,7 +433,9 @@ mkinfit <- function(mkinmod, observed,
   }
 
   if (err_mod == "const") {
-    lower["sigma"] <- 0
+    if (error_model_algorithm != "OLS") {
+      lower["sigma"] <- 0
+    }
   }
   if (err_mod == "obs") {
     index_sigma <- grep("^sigma_", names(lower))
@@ -427,11 +446,11 @@ mkinfit <- function(mkinmod, observed,
     lower["rsd_high"] <- 0
   }
 
-  # Counter for likelihood evaluations
+  # Counter for cost function evaluations
   calls = 0
-  nlogLik.current <- Inf
+  cost.current <- Inf
   out_predicted <- NA
-  data_errmod <- NA
+  current_data <- NA
 
   # Show parameter names if tracing is requested
   if(trace_parms) cat(names_optim, "\n")
@@ -441,132 +460,127 @@ mkinfit <- function(mkinmod, observed,
   # Do the fit and take the time until the hessians are calculated
   fit_time <- system.time({
     degparms <- c(state.ini.optim, transparms.optim)
-
-    if (err_mod == "const") {
-      error_model_algorithm = "OLS"
+    n_degparms <- length(degparms)
+    degparms_index <- seq(1, n_degparms)
+    errparms_index <- seq(n_degparms + 1, length.out = length(errparms))
+
+    if (error_model_algorithm == "d_3") {
+      if (!quiet) message("Directly optimising the complete model")
+      parms.start <- c(degparms, errparms)
+      fit_direct <- nlminb(parms.start, cost_function,
+        lower = lower[names(parms.start)],
+        upper = upper[names(parms.start)],
+        control = control, ...)
+      fit_direct$logLik <- - cost.current
+      if (error_model_algorithm == "direct") {
+        degparms <- fit_direct$par[degparms_index]
+        errparms <- fit_direct$par[errparms_index]
+      } else {
+        cost.current <- Inf # reset to avoid conflict with the OLS step
+      }
+    }
+    if (error_model_algorithm != "direct") {
       if (!quiet) message("Ordinary least squares optimisation")
-      fit <- nlminb(degparms, nlogLik, control = control,
+      fit <- nlminb(degparms, cost_function, control = control,
         lower = lower[names(degparms)],
         upper = upper[names(degparms)], OLS = TRUE, ...)
       degparms <- fit$par
 
       # Get the maximum likelihood estimate for sigma at the optimum parameter values
-      data_errmod$residual <- data_errmod$value.observed - data_errmod$value.predicted
-      sigma_mle <- sqrt(sum(data_errmod$residual^2)/nrow(data_errmod))
+      current_data$residual <- current_data$value.observed - current_data$value.predicted
+      sigma_mle <- sqrt(sum(current_data$residual^2)/nrow(current_data))
 
-      errparms <- c(sigma = sigma_mle)
-      nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
-      fit$logLik <- - nlogLik.current
-    } else {
-      if (error_model_algorithm == "d_3") {
-        if (!quiet) message("Directly optimising the complete model")
-        parms.start <- c(degparms, errparms)
-        fit_direct <- nlminb(parms.start, nlogLik,
-          lower = lower[names(parms.start)],
-          upper = upper[names(parms.start)],
-          control = control, ...)
-        fit_direct$logLik <- - nlogLik.current
-        nlogLik.current <- Inf # reset to avoid conflict with the OLS step
+      # Use that estimate for the constant variance, or as first guess if err_mod = "obs"
+      if (err_mod != "tc") {
+        errparms[names(errparms)] <- sigma_mle
       }
-      if (error_model_algorithm != "direct") {
-        if (!quiet) message("Ordinary least squares optimisation")
-        fit <- nlminb(degparms, nlogLik, control = control,
-          lower = lower[names(degparms)],
-          upper = upper[names(degparms)], OLS = TRUE, ...)
-        degparms <- fit$par
-        # Get the maximum likelihood estimate for sigma at the optimum parameter values
-        data_errmod$residual <- data_errmod$value.observed - data_errmod$value.predicted
-        sigma_mle <- sqrt(sum(data_errmod$residual^2)/nrow(data_errmod))
+      fit$par <- c(fit$par, errparms)
+
+      cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
+      fit$logLik <- - cost.current
+    }
+    if (error_model_algorithm %in% c("threestep", "fourstep", "d_3")) {
+      if (!quiet) message("Optimising the error model")
+      fit <- nlminb(errparms, cost_function, control = control,
+        lower = lower[names(errparms)],
+        upper = upper[names(errparms)],
+        fixed_degparms = degparms, ...)
+      errparms <- fit$par
+    }
+    if (error_model_algorithm == "fourstep") {
+      if (!quiet) message("Optimising the degradation model")
+      fit <- nlminb(degparms, cost_function, control = control,
+        lower = lower[names(degparms)],
+        upper = upper[names(degparms)],
+        fixed_errparms = errparms, ...)
+      degparms <- fit$par
+    }
+    if (error_model_algorithm %in%
+      c("direct", "twostep", "threestep", "fourstep", "d_3")) {
+      if (!quiet) message("Optimising the complete model")
+      parms.start <- c(degparms, errparms)
+      fit <- nlminb(parms.start, cost_function,
+        lower = lower[names(parms.start)],
+        upper = upper[names(parms.start)],
+        control = control, ...)
+      degparms <- fit$par[degparms_index]
+      errparms <- fit$par[errparms_index]
+      fit$logLik <- - cost.current
 
-        nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
-        fit$logLik <- - nlogLik.current
+      if (error_model_algorithm == "d_3") {
+        d_3_messages = c(
+           same = "Direct fitting and three-step fitting yield approximately the same likelihood",
+           threestep = "Three-step fitting yielded a higher likelihood than direct fitting",
+           direct = "Direct fitting yielded a higher likelihood than three-step fitting")
+        rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
+        if (rel_diff < 0.0001) {
+          if (!quiet) message(d_3_messages["same"])
+          fit$d_3_message <- d_3_messages["same"]
+        } else {
+          if (fit$logLik > fit_direct$logLik) {
+            if (!quiet) message(d_3_messages["threestep"])
+            fit$d_3_message <- d_3_messages["threestep"]
+          } else {
+            if (!quiet) message(d_3_messages["direct"])
+            fit <- fit_direct
+            fit$d_3_message <- d_3_messages["direct"]
+          }
+        }
       }
-      if (error_model_algorithm %in% c("threestep", "fourstep", "d_3")) {
+    }
+    if (err_mod != "const" & error_model_algorithm == "IRLS") {
+      reweight.diff <- 1
+      n.iter <- 0
+      errparms_last <- errparms
+
+      while (reweight.diff > reweight.tol &
+             n.iter < reweight.max.iter) {
+
         if (!quiet) message("Optimising the error model")
-        fit <- nlminb(errparms, nlogLik, control = control,
+        fit <- nlminb(errparms, cost_function, control = control,
           lower = lower[names(errparms)],
           upper = upper[names(errparms)],
           fixed_degparms = degparms, ...)
         errparms <- fit$par
-      }
-      if (error_model_algorithm == "fourstep") {
+
         if (!quiet) message("Optimising the degradation model")
-        fit <- nlminb(degparms, nlogLik, control = control,
+        fit <- nlminb(degparms, cost_function, control = control,
           lower = lower[names(degparms)],
           upper = upper[names(degparms)],
           fixed_errparms = errparms, ...)
         degparms <- fit$par
-      }
-      if (error_model_algorithm %in% c("direct", "twostep", "threestep",
-                                       "fourstep", "d_3")) {
-        if (!quiet) message("Optimising the complete model")
-        parms.start <- c(degparms, errparms)
-        fit <- nlminb(parms.start, nlogLik,
-          lower = lower[names(parms.start)],
-          upper = upper[names(parms.start)],
-          control = control, ...)
-        fit$logLik <- - nlogLik.current
 
-        d_3_messages = c(
-           same = "Direct fitting and three-step fitting yield approximately the same likelihood",
-           threestep = "Three-step fitting yielded a higher likelihood than direct fitting",
-           direct = "Direct fitting yielded a higher likelihood than three-step fitting")
-        if (error_model_algorithm == "d_3") {
-          rel_diff <- abs((fit_direct$logLik - fit$logLik))/-mean(c(fit_direct$logLik, fit$logLik))
-          if (rel_diff < 0.0001) {
-            if (!quiet) message(d_3_messages["same"])
-            fit$d_3_message <- d_3_messages["same"]
-          } else {
-            if (fit$logLik > fit_direct$logLik) {
-              if (!quiet) message(d_3_messages["threestep"])
-              fit$d_3_message <- d_3_messages["threestep"]
-            } else {
-              if (!quiet) message(d_3_messages["direct"])
-              fit <- fit_direct
-              fit$d_3_message <- d_3_messages["direct"]
-            }
-          }
-        }
-      }
-      if (err_mod != "const" & error_model_algorithm == "IRLS") {
-        reweight.diff <- 1
-        n.iter <- 0
+        reweight.diff <- dist(rbind(errparms, errparms_last))
         errparms_last <- errparms
 
-        while (reweight.diff > reweight.tol &
-               n.iter < reweight.max.iter) {
-
-          if (!quiet) message("Optimising the error model")
-          fit <- nlminb(errparms, nlogLik, control = control,
-            lower = lower[names(errparms)],
-            upper = upper[names(errparms)],
-            fixed_degparms = degparms, ...)
-          errparms <- fit$par
-
-          if (!quiet) message("Optimising the degradation model")
-          fit <- nlminb(degparms, nlogLik, control = control,
-            lower = lower[names(degparms)],
-            upper = upper[names(degparms)],
-            fixed_errparms = errparms, ...)
-          degparms <- fit$par
-
-          reweight.diff <- dist(rbind(errparms, errparms_last))
-          errparms_last <- errparms
-
-          fit$par <- c(fit$par, errparms)
-          nlogLik.current <- nlogLik(c(degparms, errparms), OLS = FALSE)
-          fit$logLik <- - nlogLik.current
-        }
+        fit$par <- c(fit$par, errparms)
+        cost.current <- cost_function(c(degparms, errparms), OLS = FALSE)
+        fit$logLik <- - cost.current
       }
     }
-    fit$error_model_algorithm <- error_model_algorithm
 
-    # We include the error model in the parameter uncertainty analysis, also
-    # for constant variance, to get a confidence interval for it
-    if (err_mod == "const") {
-      fit$par <- c(fit$par, sigma = sigma_mle)
-    }
-    fit$hessian <- try(numDeriv::hessian(nlogLik, fit$par, update_data = FALSE), silent = TRUE)
+    fit$hessian <- try(numDeriv::hessian(cost_function, c(degparms, errparms), OLS = FALSE,
+        update_data = FALSE), silent = TRUE)
 
     # Backtransform parameters
     bparms.optim = backtransform_odeparms(fit$par, mkinmod,
@@ -575,10 +589,12 @@ mkinfit <- function(mkinmod, observed,
     bparms.fixed = c(state.ini.fixed, parms.fixed)
     bparms.all = c(bparms.optim, parms.fixed)
 
-    fit$hessian_notrans <- try(numDeriv::hessian(nlogLik, c(bparms.optim, fit$par[names(errparms)]),
-                                       trans = FALSE, update_data = FALSE), silent = TRUE)
+    fit$hessian_notrans <- try(numDeriv::hessian(cost_function, c(bparms.all, errparms),
+        OLS = FALSE, trans = FALSE, update_data = FALSE), silent = TRUE)
   })
 
+  fit$error_model_algorithm <- error_model_algorithm
+
   if (fit$convergence != 0) {
     fit$warning = paste0("Optimisation did not converge:\n", fit$message)
     warning(fit$warning)
@@ -604,18 +620,24 @@ mkinfit <- function(mkinmod, observed,
   fit$obs_vars <- obs_vars
   fit$predicted <- out_predicted
 
-  # Attach the negative log-likelihood function for post-hoc parameter uncertainty analysis
-  fit$nlogLik <- nlogLik
+  # Residual sum of squares as a function of the fitted parameters
+  fit$rss <- function(P) cost_function(P, OLS = TRUE, update_data = FALSE)
+
+  # Log-likelihood with possibility to fix degparms or errparms
+  fit$ll <- function(P, fixed_degparms = FALSE, fixed_errparms = FALSE) {
+    - cost_function(P, fixed_degparms = fixed_degparms,
+      fixed_errparms = fixed_errparms, OLS = FALSE, update_data = FALSE)
+  }
 
   # Collect initial parameter values in three dataframes
   fit$start <- data.frame(value = c(state.ini.optim,
-                                    parms.optim, errparms))
+                                    parms.optim, errparms_optim))
   fit$start$type = c(rep("state", length(state.ini.optim)),
                      rep("deparm", length(parms.optim)),
-                     rep("error", length(errparms)))
+                     rep("error", length(errparms_optim)))
 
   fit$start_transformed = data.frame(
-      value = c(state.ini.optim, transparms.optim, errparms),
+      value = c(state.ini.optim, transparms.optim, errparms_optim),
       lower = lower,
       upper = upper)
 
@@ -624,14 +646,14 @@ mkinfit <- function(mkinmod, observed,
                      rep("deparm", length(parms.fixed)))
 
   # Sort observed, predicted and residuals
-  data_errmod$name <- ordered(data_errmod$name, levels = obs_vars)
+  current_data$name <- ordered(current_data$name, levels = obs_vars)
 
-  data <- data_errmod[order(data_errmod$name, data_errmod$time), ]
+  ordered_data <- current_data[order(current_data$name, current_data$time), ]
 
-  fit$data <- data.frame(time = data$time,
-                         variable = data$name,
-                         observed = data$value.observed,
-                         predicted = data$value.predicted)
+  fit$data <- data.frame(time = ordered_data$time,
+                         variable = ordered_data$name,
+                         observed = ordered_data$value.observed,
+                         predicted = ordered_data$value.predicted)
 
   fit$data$residual <- fit$data$observed - fit$data$predicted
 
@@ -649,8 +671,8 @@ mkinfit <- function(mkinmod, observed,
                         state.ini.fixed)
   names(fit$bparms.state) <- gsub("_0$", "", names(fit$bparms.state))
 
-  fit$errparms.optim <- fit$par[names(errparms)]
-  fit$df.residual <- n_observed - n_optim
+  fit$errparms <- errparms
+  fit$df.residual <- n_observed - length(c(degparms, errparms))
 
   fit$date <- date()
   fit$version <- as.character(utils::packageVersion("mkin"))
@@ -664,7 +686,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
   param  <- object$par
   pnames <- names(param)
   bpnames <- names(object$bparms.optim)
-  epnames <- names(object$errparms.optim)
+  epnames <- names(object$errparms)
   p      <- length(param)
   mod_vars <- names(object$mkinmod$diffs)
   covar  <- try(solve(object$hessian), silent = TRUE)
@@ -736,9 +758,9 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
     use_of_ff = object$mkinmod$use_of_ff,
     error_model_algorithm = object$error_model_algorithm,
     df = c(p, rdf),
-    cov.unscaled = covar,
+    covar = covar,
+    covar_notrans = covar_notrans,
     err_mod = object$err_mod,
-    #cov.scaled = covar * resvar,
     niter = object$iterations,
     calls = object$calls,
     time = object$time,
@@ -760,8 +782,8 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
   ans$errmin <- mkinerrmin(object, alpha = 0.05)
 
   if (object$calls > 0) {
-    if (!is.null(ans$cov.unscaled)){
-      Corr <- cov2cor(ans$cov.unscaled)
+    if (!is.null(ans$covar)){
+      Corr <- cov2cor(ans$covar)
       rownames(Corr) <- colnames(Corr) <- rownames(ans$par)
       ans$Corr <- Corr
     } else {
@@ -831,7 +853,7 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
 
   if (x$calls > 0) {
     cat("\nParameter correlation:\n")
-    if (!is.null(x$cov.unscaled)){
+    if (!is.null(x$covar)){
       print(x$Corr, digits = digits, ...)
     } else {
       cat("No covariance matrix")
-- 
cgit v1.2.1