diff options
Diffstat (limited to 'R/mkinfit.R')
| -rw-r--r-- | R/mkinfit.R | 499 |
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:
|