Direct error model fitting works

- No IRLS required
- Removed optimization algorithms other than Port
- Removed the dependency on FME
- Fitting the error model 'obs' is much faster for the FOCUS_2006_D
  dataset and the FOMC_SFO model (1 second versus 3.4 seconds)
- Vignettes build slower. Compiled models needs 3 minutes instead of 1.5
- For other vignettes, the trend is less clear. Some fits are faster,
  even for error_model = "const". FOCUS_Z is faster (34.9 s versus
  44.1 s)
- Standard errors and confidence intervals are slightly smaller
- Removed code for plotting during the fit, as I hardly ever used it
- Merged the two cost functions (using transformed and untransformed
  parameters) into one log-likelihood function

1 files changed, 133 insertions, 366 deletions

diff --git a/R/mkinfit.R b/R/mkinfit.R
index 4ac54ce2..796b25c7 100644
--- a/R/mkinfit.R
+++ b/R/mkinfit.R
@@ -30,19 +30,12 @@ mkinfit <- function(mkinmod, observed,
   solution_type = c("auto", "analytical", "eigen", "deSolve"),

   method.ode = "lsoda",

   use_compiled = "auto",

-  method.modFit = c("Port", "Marq", "SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B"),

-  maxit.modFit = "auto",

-  control.modFit = list(),

+  control = list(eval.max = 300, iter.max = 200),

   transform_rates = TRUE,

   transform_fractions = TRUE,

   plot = FALSE, quiet = FALSE,

-  err = NULL,

-  weight = c("none", "manual", "std", "mean", "tc"),

-  tc = c(sigma_low = 0.5, rsd_high = 0.07),

-  scaleVar = FALSE,

   atol = 1e-8, rtol = 1e-10, n.outtimes = 100,

-  reweight.method = NULL,

-  reweight.tol = 1e-8, reweight.max.iter = 10,

+  error_model = c("const", "obs", "tc", "obs_tc"),

   trace_parms = FALSE,

   ...)

 {

@@ -61,19 +54,6 @@ mkinfit <- function(mkinmod, observed,
     }

   }

 

-  # Check optimisation method and set maximum number of iterations if specified

-  method.modFit = match.arg(method.modFit)

-  if (maxit.modFit != "auto") {

-    if (method.modFit == "Marq") control.modFit$maxiter = maxit.modFit

-    if (method.modFit == "Port") {

-      control.modFit$iter.max = maxit.modFit

-      control.modFit$eval.max = maxit.modFit

-    }

-    if (method.modFit %in% c("SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {

-        control.modFit$maxit = maxit.modFit

-    }

-  }

-

   # Get the names of the state variables in the model

   mod_vars <- names(mkinmod$diffs)

 

@@ -96,6 +76,9 @@ mkinfit <- function(mkinmod, observed,
     observed$time <- observed$time - t_of_max

   }

 

+  # Number observations used for fitting

+  n_observed <- nrow(observed)

+

   # Define starting values for parameters where not specified by the user

   if (parms.ini[[1]] == "auto") parms.ini = vector()

 

@@ -254,32 +237,42 @@ mkinfit <- function(mkinmod, observed,
     }

   }

 

+  # Get the error model and specify starting values as well as fixed components

+  err_mod <- match.arg(error_model)

+  if (err_mod == "const") {

+    errparms = c(sigma = 1)

+  }

+  if (err_mod == "obs") {

+    errparms <- rep(1, length(obs_vars))

+    names(errparms) = paste0("sigma_", obs_vars)

+  }

+  if (err_mod == "tc") {

+    errparms <- c(sigma_low = 0.5, rsd_high = 0.07)

+  }

+  if (err_mod == "obs_tc") {

+    errparms <- rep(1, length(obs_vars))

+    names(errparms) = paste0("sigma_", obs_vars)

+    errparms <- c(errparms, a = 0.1)

+  }

+

   # Define outtimes for model solution.

   # Include time points at which observed data are available

   outtimes = sort(unique(c(observed$time, seq(min(observed$time),

                                               max(observed$time),

                                               length.out = n.outtimes))))

 

-  weight.ini <- weight <- match.arg(weight)

-  if (weight.ini == "tc") {

-     observed$err = sigma_twocomp(observed$value, tc["sigma_low"], tc["rsd_high"])

-     err <- "err"

-  } else {

-    if (!is.null(err)) weight.ini = "manual"

-  }

-

-  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 for optimisation, including (back)transformations

-  cost <- function(P)

+  # Define log-likelihood function for optimisation, including (back)transformations

+  nlogLik <- function(P, trans = TRUE)

   {

     assign("calls", calls+1, inherits=TRUE) # Increase the model solution counter

 

     # Trace parameter values if requested

     if(trace_parms) cat(P, "\n")

 

+    # Initial states for t0

     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)

@@ -288,13 +281,18 @@ mkinfit <- function(mkinmod, observed,
       names(odeini) <- state.ini.fixed.boxnames

     }

 

-    odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], transparms.fixed)

+    odeparms.optim <- P[(length(state.ini.optim) + 1):(length(P) - length(errparms))]

 

-    parms <- backtransform_odeparms(odeparms, mkinmod,

-                                    transform_rates = transform_rates,

-                                    transform_fractions = transform_fractions)

+    if (trans == TRUE) {

+      odeparms <- c(odeparms.optim, transparms.fixed)

+      parms <- backtransform_odeparms(odeparms, mkinmod,

+                                      transform_rates = transform_rates,

+                                      transform_fractions = transform_fractions)

+    } else {

+      parms <- c(odeparms.optim, parms.fixed)

+    }

 

-    # Solve the system with current transformed parameter values

+    # Solve the system with current parameter values

     out <- mkinpredict(mkinmod, parms,

                        odeini, outtimes,

                        solution_type = solution_type,

@@ -302,82 +300,55 @@ mkinfit <- function(mkinmod, observed,
                        method.ode = method.ode,

                        atol = atol, rtol = rtol, ...)

 

-    assign("out_predicted", out, inherits=TRUE)

-

-    mC <- modCost(out, observed, y = "value",

-      err = err, weight = weight, scaleVar = scaleVar)

-

-    # Report and/or plot if the model is improved

-    if (mC$model < cost.old) {

-      if(!quiet) cat("Model cost at call ", calls, ": ", mC$model, "\n")

-

-      # Plot the data and current model output if requested

-      if(plot) {

-        outtimes_plot = seq(min(observed$time), max(observed$time), length.out=100)

-

-        out_plot <- mkinpredict(mkinmod, parms,

-                                odeini, outtimes_plot,

-                                solution_type = solution_type,

-                                use_compiled = use_compiled,

-                                method.ode = method.ode,

-                                atol = atol, rtol = rtol, ...)

-

-        plot(0, type="n",

-          xlim = range(observed$time), ylim = c(0, max(observed$value, na.rm=TRUE)),

-          xlab = "Time", ylab = "Observed")

-        col_obs <- pch_obs <- 1:length(obs_vars)

-        lty_obs <- rep(1, length(obs_vars))

-        names(col_obs) <- names(pch_obs) <- names(lty_obs) <- obs_vars

-        for (obs_var in obs_vars) {

-          points(subset(observed, name == obs_var, c(time, value)),

-                 pch = pch_obs[obs_var], col = col_obs[obs_var])

-        }

-        matlines(out_plot$time, out_plot[-1], col = col_obs, lty = lty_obs)

-        legend("topright", inset=c(0.05, 0.05), legend=obs_vars,

-          col=col_obs, pch=pch_obs, lty=1:length(pch_obs))

-      }

+    out_long <- mkin_wide_to_long(out, time = "time")

 

-      assign("cost.old", mC$model, inherits=TRUE)

-    }

-    return(mC)

-  }

+    assign("out_predicted", out_long, inherits=TRUE)

 

-  # 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

+    if (err_mod == "const") {

+      observed$std <- P["sigma"]

+    }

+    if (err_mod == "obs") {

+      std_names <- paste0("sigma_", observed$name)

+      observed$std <- P[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(P["sigma_low"]^2 + tmp$value.y^2 * P["rsd_high"]^2)

+    }

+    if (err_mod == "obs_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), ]

+      std_names <- paste0("sigma_", observed$name)

+      observed$std <- sqrt(P[std_names] * (1 + P["a"] * tmp$value.y^2))

     }

 

-    odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)

-

-    # Solve the system with current 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)

+    data_log_lik <- merge(observed[c("name", "time", "value", "std")], out_long,

+                         by = c("name", "time"))

 

-    return(mC)

+    nlogLik <- - with(data_log_lik,

+      sum(dnorm(x = value.x, mean = value.y, sd = std, log = TRUE)))

+    return(nlogLik)

   }

 

+  n_optim <- length(c(state.ini.optim, transparms.optim, errparms))

+  names_optim <- c(names(state.ini.optim),

+                   names(transparms.optim),

+                   names(errparms))

+

   # 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))

+  # for which no internal transformation is requested in the call to mkinfit

+  # and for error model parameters

+  lower <- rep(-Inf, n_optim)

+  upper <- rep(Inf, n_optim)

+  names(lower) <- names(upper) <- names_optim

 

   # IORE exponentes are not transformed, but need a lower bound of zero

   index_N <- grep("^N", names(lower))

   lower[index_N] <- 0

+

   if (!transform_rates) {

     index_k <- grep("^k_", names(lower))

     lower[index_k] <- 0

@@ -396,151 +367,44 @@ mkinfit <- function(mkinmod, observed,
     upper[other_fraction_parms] <- 1

   }

 

+  if (err_mod == "const") {

+    lower["sigma"] <- 0

+  }

+  if (err_mod == "obs") {

+    index_sigma <- grep("^sigma_", names(lower))

+    lower[index_sigma] <- 0

+  }

+  if (err_mod == "tc") {

+    lower["sigma_low"] <- 0

+    lower["rsd_high"] <- 0

+  }

+  if (err_mod == "obs_tc") {

+    index_sigma <- grep("^sigma_", names(lower))

+    lower[index_sigma] <- 0

+    lower["a"] <- 0

+  }

+

   # Show parameter names if tracing is requested

-  if(trace_parms) cat(names(c(state.ini.optim, transparms.optim)), "\n")

+  if(trace_parms) cat(names_optim, "\n")

 

   # Do the fit and take the time

   fit_time <- system.time({

-    fit <- modFit(cost, c(state.ini.optim, transparms.optim),

-                  method = method.modFit, control = control.modFit,

-                  lower = lower, upper = upper, ...)

-

-    # Reiterate the fit until convergence of the variance components (IRLS)

-    # if requested by the user

-

-    if (!is.null(reweight.method)) {

-      if (! reweight.method %in% c("obs", "tc")) stop("Only reweighting methods 'obs' and 'tc' are implemented")

-

-      if (reweight.method  == "obs") {

-        tc_fit <- NA

-        if(!quiet) {

-          cat("IRLS based on variance estimates for each observed variable\n")

-          cat("Initial variance estimates are:\n")

-          print(signif(fit$var_ms_unweighted, 8))

-        }

-      }

-      if (reweight.method  == "tc") {

-        tc_fit <- .fit_error_model_mad_obs(cost(fit$par)$residuals, tc, 0)

-

-        if (is.character(tc_fit)) {

-          if (!quiet) {

-            cat(tc_fit, ".\n", "No reweighting will be performed.")

-          }

-          tc_fitted <- c(sigma_low = NA, rsd_high = NA)

-        } else {

-          tc_fitted <- coef(tc_fit)

-          if(!quiet) {

-            cat("IRLS based on variance estimates according to the two component error model\n")

-            cat("Initial variance components are:\n")

-            print(signif(tc_fitted))

-          }

-        }

-      }

-      reweight.diff = 1

-      n.iter <- 0

-      if (!is.null(err)) observed$err.ini <- observed[[err]]

-      err = "err.irls"

-

-      while (reweight.diff > reweight.tol &

-             n.iter < reweight.max.iter &

-             !is.character(tc_fit)) {

-        n.iter <- n.iter + 1

-        # Store squared residual predictors used for weighting in sr_old and define new weights

-        if (reweight.method == "obs") {

-          sr_old <- fit$var_ms_unweighted

-          observed[err] <- sqrt(fit$var_ms_unweighted[as.character(observed$name)])

-        }

-        if (reweight.method == "tc") {

-          sr_old <- tc_fitted

-

-          tmp_predicted <- mkin_wide_to_long(out_predicted, time = "time")

-          tmp_data <- suppressMessages(join(observed, tmp_predicted, by = c("time", "name")))

-

-          #observed[err] <- predict(tc_fit, newdata = data.frame(mod = tmp_data[[4]]))

-          observed[err] <- predict(tc_fit, newdata = data.frame(obs = observed$value))

-

-        }

-        fit <- modFit(cost, fit$par, method = method.modFit,

-                      control = control.modFit, lower = lower, upper = upper, ...)

+    fit <- nlminb(c(state.ini.optim, transparms.optim, errparms),

+                  nlogLik, control = control,

+                  lower = lower, upper = upper, ...)})

 

-        if (reweight.method == "obs") {

-          sr_new <- fit$var_ms_unweighted

-        }

-        if (reweight.method == "tc") {

-          tc_fit <- .fit_error_model_mad_obs(cost(fit$par)$residuals, tc_fitted, n.iter)

-

-          if (is.character(tc_fit)) {

-            if (!quiet) {

-              cat(tc_fit, ".\n")

-            }

-            break

-          } else {

-            tc_fitted <- coef(tc_fit)

-            sr_new <- tc_fitted

-          }

-        }

-

-        reweight.diff = sum((sr_new - sr_old)^2)

-        if (!quiet) {

-          cat("Iteration", n.iter, "yields variance estimates:\n")

-          print(signif(sr_new, 8))

-          cat("Sum of squared differences to last variance (component) estimates:",

-              signif(reweight.diff, 2), "\n")

-        }

-      }

-    }

-  })

-

-  # Check for convergence

-  if (method.modFit == "Marq") {

-    if (!fit$info %in% c(1, 2, 3)) {

-      fit$warning = paste0("Optimisation by method ", method.modFit,

-                           " did not converge.\n",

-                           "The message returned by nls.lm is:\n",

-                                    fit$message)

-      warning(fit$warning)

-    }

-    else {

-      if(!quiet) cat("Optimisation by method", method.modFit, "successfully terminated.\n")

-    }

-  }

-  if (method.modFit == "Port") {

-    if (fit$convergence != 0) {

-      fit$warning = paste0("Optimisation by method ", method.modFit,

-                           " did not converge:\n",

-                           if(is.character(fit$counts)) fit$counts # FME bug

-                           else fit$message)

-      warning(fit$warning)

-    } else {

-      if(!quiet) cat("Optimisation by method", method.modFit, "successfully terminated.\n")

-    }

-  }

-  if (method.modFit %in% c("SANN", "Nelder-Mead", "BFGS", "CG", "L-BFGS-B")) {

-    if (fit$convergence != 0) {

-      fit$warning = paste0("Optimisation by method ", method.modFit,

-                           " did not converge.\n",

-                           "Convergence code returned by optim() is", fit$convergence)

-      warning(fit$warning)

-    } else {

-      if(!quiet) cat("Optimisation by method", method.modFit, "successfully terminated.\n")

-    }

-  }

-

-  # Return number of iterations for SANN method and alert user to check if

-  # the approximation to the optimum is sufficient

-  if (method.modFit == "SANN") {

-    fit$iter = maxit.modFit

-    fit$warning <- paste0("Termination of the SANN algorithm does not imply convergence.\n",

-      "Make sure the approximation of the optimum is adequate.")

+  if (fit$convergence != 0) {

+    fit$warning = paste0("Optimisation did not converge:\n", fit$message)

     warning(fit$warning)

+  } else {

+    if(!quiet) cat("Optimisation successfully terminated.\n")

   }

 

   # We need to return some more data for summary and plotting

   fit$solution_type <- solution_type

   fit$transform_rates <- transform_rates

   fit$transform_fractions <- transform_fractions

-  fit$method.modFit <- method.modFit

-  fit$maxit.modFit <- maxit.modFit

+  fit$control <- control

   fit$calls <- calls

   fit$time <- fit_time

 

@@ -550,40 +414,30 @@ mkinfit <- function(mkinmod, observed,
   # We need data and predictions for summary and plotting

   fit$observed <- observed

   fit$obs_vars <- obs_vars

-  fit$predicted <- mkin_wide_to_long(out_predicted, time = "time")

+  fit$predicted <- out_predicted

 

   # Backtransform parameters

   bparms.optim = backtransform_odeparms(fit$par, fit$mkinmod,

-                                        transform_rates = transform_rates,

-                                        transform_fractions = transform_fractions)

+    transform_rates = transform_rates,

+    transform_fractions = transform_fractions)

   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 <- try(gradient(function(p, ...) cost_notrans(p)$residuals$res,

-                              bparms.optim, centered = TRUE), silent = TRUE)

-  if (inherits(Jac_notrans, "try-error")) {

-    warning("Calculation of the Jacobian failed for the cost function of the untransformed model.\n",

-            "No t-test results will be available")

-    fit$hessian_notrans <- NA

-  } else {

-    fit$hessian_notrans <- 2 * t(Jac_notrans) %*% Jac_notrans

-  }

+  # Attach the negative log-likelihood function for post-hoc parameter uncertainty analysis

+  fit$nlogLik <- nlogLik

+

+  fit$hessian <- hessian(nlogLik, fit$par)

+  fit$hessian_notrans <- hessian(nlogLik, c(bparms.optim, fit$par[names(errparms)]), trans = FALSE)

 

   # Collect initial parameter values in three dataframes

   fit$start <- data.frame(value = c(state.ini.optim,

-                                    parms.optim))

+                                    parms.optim, errparms))

   fit$start$type = c(rep("state", length(state.ini.optim)),

-                     rep("deparm", length(parms.optim)))

+                     rep("deparm", length(parms.optim)),

+                     rep("error", length(errparms)))

 

   fit$start_transformed = data.frame(

-      value = c(state.ini.optim, transparms.optim),

+      value = c(state.ini.optim, transparms.optim, errparms),

       lower = lower,

       upper = upper)

 

@@ -596,28 +450,16 @@ mkinfit <- function(mkinmod, observed,
   data$name <- ordered(data$name, levels = obs_vars)

   data <- data[order(data$name, data$time), ]

 

-  # Add a column named "value" holding the observed values for the case

-  # that this column was selected for manual weighting, so it can be

-  # shown in the summary as "err"

-  data$value <- data$value.x

-

   fit$data <- data.frame(time = data$time,

                          variable = data$name,

                          observed = data$value.x,

                          predicted = data$value.y)

 

   fit$data$residual <- fit$data$observed - fit$data$predicted

-  if (!is.null(data$err.ini)) fit$data$err.ini <- data$err.ini

-  if (!is.null(err)) fit$data[["err"]] <- data[[err]]

 

   fit$atol <- atol

   fit$rtol <- rtol

-  fit$weight.ini <- weight.ini

-  fit$tc.ini <- tc

-  fit$reweight.method <- reweight.method

-  fit$reweight.tol <- reweight.tol

-  fit$reweight.max.iter <- reweight.max.iter

-  if (exists("tc_fitted")) fit$tc_fitted <- tc_fitted

+  fit$err_mod <- err_mod

 

   # Return different sets of backtransformed parameters for summary and plotting

   fit$bparms.optim <- bparms.optim

@@ -629,6 +471,9 @@ 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$date <- date()

   fit$version <- as.character(utils::packageVersion("mkin"))

   fit$Rversion <- paste(R.version$major, R.version$minor, sep=".")

@@ -641,39 +486,41 @@ 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)

   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

+  covar  <- try(solve(object$hessian), silent = TRUE)

+  covar_notrans  <- try(solve(object$hessian_notrans), silent = TRUE)

+  rdf <- object$df.residual

+

   if (!is.numeric(covar) | is.na(covar[1])) {

     covar <- NULL

     se <- lci <- uci <- rep(NA, p)

   } else {

     rownames(covar) <- colnames(covar) <- pnames

-    se     <- sqrt(diag(covar) * resvar)

+    se     <- sqrt(diag(covar))

     lci    <- param + qt(alpha/2, rdf) * se

     uci    <- param + qt(1-alpha/2, rdf) * se

   }

-

+  

+  beparms.optim <- c(object$bparms.optim, object$par[epnames])

   if (!is.numeric(covar_notrans) | is.na(covar_notrans[1])) {

     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)

+    rownames(covar_notrans) <- colnames(covar_notrans) <- c(bpnames, epnames)

+    se_notrans <- sqrt(diag(covar_notrans))

+    tval  <- beparms.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)

   dimnames(param) <- list(pnames, c("Estimate", "Std. Error", "Lower", "Upper"))

 

-  bparam <- cbind(Estimate = object$bparms.optim, se_notrans, "t value" = tval, "Pr(>t)" = pval, Lower = NA, Upper = NA)

+  bparam <- cbind(Estimate = beparms.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).

@@ -699,6 +546,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
       bparam[names(bpu), "Upper"] <- bpu

     }

   }

+  bparam[epnames, c("Lower", "Upper")] <- param[epnames, c("Lower", "Upper")]

 

   ans <- list(

     version = as.character(utils::packageVersion("mkin")),

@@ -706,26 +554,14 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05,
     date.fit = object$date,

     date.summary = date(),

     solution_type = object$solution_type,

-    method.modFit = object$method.modFit,

     warning = object$warning,

     use_of_ff = object$mkinmod$use_of_ff,

-    weight.ini = object$weight.ini,

-    reweight.method = object$reweight.method,

-    tc.ini = object$tc.ini,

-    var_ms_unweighted = object$var_ms_unweighted,

-    tc_fitted = object$tc_fitted,

-    residuals = object$residuals,

-    residualVariance = resvar,

-    sigma = sqrt(resvar),

-    modVariance = modVariance,

     df = c(p, rdf),

     cov.unscaled = covar,

-    cov.scaled = covar * resvar,

-    info = object$info,

+    #cov.scaled = covar * resvar,

     niter = object$iterations,

     calls = object$calls,

     time = object$time,

-    stopmess = message,

     par = param,

     bpar = bparam)

 

@@ -789,24 +625,8 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
   cat("\nFitted with method", x$method.modFit,

       "using", x$calls, "model solutions performed in", x$time[["elapsed"]],  "s\n")

 

-  cat("\nWeighting:", x$weight.ini)

-  if (x$weight.ini == "tc")  {

-    cat(" with variance components\n")

-    print(x$tc.ini)

-  } else {

-    cat ("\n")

-  }

-  if(!is.null(x$reweight.method)) {

-    cat("\nIterative reweighting with method", x$reweight.method, "\n")

-    if (x$reweight.method == "obs") {

-      cat("Final mean squared residuals of observed variables:\n")

-      print(x$var_ms_unweighted)

-    }

-    if (x$reweight.method == "tc") {

-      cat("Final components of fitted standard deviation:\n")

-      print(x$tc_fitted)

-    }

-  }

+  cat("\nError model:\n")

+  print(x$err_mod)

 

   cat("\nStarting values for parameters to be optimised:\n")

   print(x$start)

@@ -830,9 +650,6 @@ 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")

   cat("Confidence intervals for internally transformed parameters are asymmetric.\n")

   if ((x$version) < "0.9-36") {

@@ -874,54 +691,4 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), .
 

   invisible(x)

 }

-

-# Fit the median absolute deviation against the observed values,

-# using the current error model for weighting

-.fit_error_model_mad_obs <- function(tmp_res, tc, iteration) {

-  mad_agg <- aggregate(tmp_res$res.unweighted,

-                       by = list(name = tmp_res$name, time = tmp_res$x),

-                       FUN = function(x) mad(x, center = 0))

-  names(mad_agg) <- c("name", "time", "mad")

-  error_data <- suppressMessages(

-    join(data.frame(name = tmp_res$name,

-                    time = tmp_res$x,

-                    obs = tmp_res$obs),

-         mad_agg))

-  error_data_complete <- na.omit(error_data)

-

-  tc_fit <- tryCatch(

-    nls(mad ~ sigma_twocomp(obs, sigma_low, rsd_high),

-      start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),

-      weights = 1/sigma_twocomp(error_data_complete$obs,

-                                tc["sigma_low"],

-                                tc["rsd_high"])^2,

-      data = error_data_complete,

-      lower = 0,

-      algorithm = "port"),

-    error = function(e) paste("Fitting the error model failed in iteration", iteration))

-  return(tc_fit)

-}

-# Alternative way to fit the error model, fitting to modelled instead of

-# observed values

-# .fit_error_model_mad_mod <- function(tmp_res, tc) {

-#   mad_agg <- aggregate(tmp_res$res.unweighted,

-#                        by = list(name = tmp_res$name, time = tmp_res$x),

-#                        FUN = function(x) mad(x, center = 0))

-#   names(mad_agg) <- c("name", "time", "mad")

-#   mod_agg <- aggregate(tmp_res$mod,

-#                        by = list(name = tmp_res$name, time = tmp_res$x),

-#                        FUN = mean)

-#   names(mod_agg) <- c("name", "time", "mod")

-#   mod_mad <- merge(mod_agg, mad_agg)

-# 

-#   tc_fit <- tryCatch(

-#     nls(mad ~ sigma_twocomp(mod, sigma_low, rsd_high),

-#       start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),

-#       data = mod_mad,

-#       weights = 1/mod_mad$mad,

-#       lower = 0,

-#       algorithm = "port"),

-#     error = "Fitting the error model failed in iteration")

-#   return(tc_fit)

-# }

 # vim: set ts=2 sw=2 expandtab: