diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/mkinfit.R | 27 | ||||
| -rw-r--r-- | R/summary.mkinfit.R | 16 | ||||
| -rw-r--r-- | R/summary_DFOP_FOCUS_C.txt | 82 |
3 files changed, 111 insertions, 14 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R index d7b1b7f4..ec2d3412 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -233,7 +233,7 @@ if(getRversion() >= '2.15.1') utils::globalVariables(c("name", "time", "value")) #' fit.FOMC = mkinfit("FOMC", FOCUS_2006_D, quiet = TRUE, error_model = "tc") #' fit.FOMC_SFO <- mkinfit(FOMC_SFO, FOCUS_D, quiet = TRUE, #' parms.ini = fit.FOMC$bparms.ode, error_model = "tc") -#' +#' } #' @export mkinfit <- function(mkinmod, observed, parms.ini = "auto", @@ -258,6 +258,8 @@ mkinfit <- function(mkinmod, observed, { call <- match.call() + summary_warnings <- character() + # Derive the name used for the model if (is.character(mkinmod)) mkinmod_name <- mkinmod else mkinmod_name <- deparse(substitute(mkinmod)) @@ -289,7 +291,9 @@ mkinfit <- function(mkinmod, observed, # Also remove zero values to avoid instabilities (e.g. of the 'tc' error model) if (any(observed$value == 0)) { - warning("Observations with value of zero were removed from the data") + zero_warning <- "Observations with value of zero were removed from the data" + summary_warnings <- c(summary_warnings, zero_warning) + warning(zero_warning) observed <- subset(observed, value != 0) } @@ -848,8 +852,9 @@ mkinfit <- function(mkinmod, observed, fit$error_model_algorithm <- error_model_algorithm if (fit$convergence != 0) { - fit$warning = paste0("Optimisation did not converge:\n", fit$message) - warning(fit$warning) + convergence_warning = paste0("Optimisation did not converge:\n", fit$message) + summary_warnings <- c(warnings, convergence_warning) + warning(convergence_warning) } else { if(!quiet) message("Optimisation successfully terminated.\n") } @@ -918,14 +923,22 @@ mkinfit <- function(mkinmod, observed, fit$errparms <- errparms fit$df.residual <- n_observed - length(c(degparms, errparms)) + # Assign the class here so method dispatch works for residuals + class(fit) <- c("mkinfit") + # Check for normal distribution of residuals - fit$shapiro.p <- shapiro.test(residuals.mkinfit(fit, standardized = TRUE))$p.value - if (fit$shapiro.p < 0.05) warning("The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05") + fit$shapiro.p <- shapiro.test(residuals(fit, standardized = TRUE))$p.value + if (fit$shapiro.p < 0.05) { + shapiro_warning <- paste("Shapiro-Wilk test for standardized residuals: p = ", signif(fit$shapiro.p, 3)) + warning(shapiro_warning) + summary_warnings <- c(summary_warnings, shapiro_warning) + } + + fit$summary_warnings <- summary_warnings fit$date <- date() fit$version <- as.character(utils::packageVersion("mkin")) fit$Rversion <- paste(R.version$major, R.version$minor, sep=".") - class(fit) <- c("mkinfit") return(fit) } diff --git a/R/summary.mkinfit.R b/R/summary.mkinfit.R index 2c291ffd..f9858c32 100644 --- a/R/summary.mkinfit.R +++ b/R/summary.mkinfit.R @@ -122,7 +122,7 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, date.fit = object$date, date.summary = date(), solution_type = object$solution_type, - warning = object$warning, + warnings = object$summary_warnings, use_of_ff = object$mkinmod$use_of_ff, error_model_algorithm = object$error_model_algorithm, df = c(p, rdf), @@ -190,8 +190,6 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . cat("Date of fit: ", x$date.fit, "\n") cat("Date of summary:", x$date.summary, "\n") - if (!is.null(x$warning)) cat("\n\nWarning:", x$warning, "\n\n") - cat("\nEquations:\n") nice_diffs <- gsub("^(d.*) =", "\\1/dt =", x[["diffs"]]) writeLines(strwrap(nice_diffs, exdent = 11)) @@ -223,16 +221,20 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . if(length(x$fixed$value) == 0) cat("None\n") else print(x$fixed) + # We used to only have one warning - show this for summarising old objects + if (!is.null(x[["warning"]])) cat("\n\nWarning:", x$warning, "\n\n") + + if (length(x$warnings) > 0) { + cat("\n\nWarning(s):", "\n") + cat(x$warnings, sep = "\n") + } + if (!is.null(x$AIC)) { cat("\nResults:\n\n") print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik, row.names = " ")) } - if (!is.null(x$shapiro.p) && x$shapiro.p < 0.05) { - warning("The p-value for the Shapiro-Wilk test of normality on standardized residuals is < 0.05") - } - cat("\nOptimised, transformed parameters with symmetric confidence intervals:\n") print(signif(x$par, digits = digits)) diff --git a/R/summary_DFOP_FOCUS_C.txt b/R/summary_DFOP_FOCUS_C.txt new file mode 100644 index 00000000..ab64a588 --- /dev/null +++ b/R/summary_DFOP_FOCUS_C.txt @@ -0,0 +1,82 @@ +mkin version used for fitting: Dummy 0.0 for testing +R version used for fitting: Dummy R version for testing +Date of fit: Dummy date for testing +Date of summary: Dummy date for testing + +Equations: +d_parent/dt = - ((k1 * g * exp(-k1 * time) + k2 * (1 - g) * exp(-k2 * + time)) / (g * exp(-k1 * time) + (1 - g) * exp(-k2 * time))) + * parent + +Model predictions using solution type analytical + +Fitted using test 0 model solutions performed in test time 0 s + +Error model: Constant variance + +Error model algorithm: OLS + +Starting values for parameters to be optimised: + value type +parent_0 85.10 state +k1 0.10 deparm +k2 0.01 deparm +g 0.50 deparm + +Starting values for the transformed parameters actually optimised: + value lower upper +parent_0 85.100000 -Inf Inf +log_k1 -2.302585 -Inf Inf +log_k2 -4.605170 -Inf Inf +g_ilr 0.000000 -Inf Inf + +Fixed parameter values: +None + +Results: + + AIC BIC logLik + 29.02372 30.00984 -9.511861 + +Optimised, transformed parameters with symmetric confidence intervals: + Estimate Std. Error Lower Upper +parent_0 85.0000 0.66620 83.1500 86.8500 +log_k1 -0.7775 0.03380 -0.8713 -0.6836 +log_k2 -4.0260 0.13100 -4.3890 -3.6620 +g_ilr 1.2490 0.05811 1.0870 1.4100 +sigma 0.6962 0.16410 0.2406 1.1520 + +Parameter correlation: +[1] "Correlation matrix is platform dependent, not tested" + +Backtransformed parameters: +Confidence intervals for internally transformed parameters are asymmetric. +t-test (unrealistically) based on the assumption of normal distribution +for estimators of untransformed parameters. + Estimate t value Pr(>t) Lower Upper +parent_0 85.00000 127.600 1.131e-08 83.15000 86.85000 +k1 0.45960 29.580 3.887e-06 0.41840 0.50480 +k2 0.01785 7.636 7.901e-04 0.01241 0.02568 +g 0.85390 83.310 6.221e-08 0.82310 0.88020 +sigma 0.69620 4.243 6.618e-03 0.24060 1.15200 + +FOCUS Chi2 error levels in percent: + err.min n.optim df +All data 2.661 4 5 +parent 2.661 4 5 + +Estimated disappearance times: + DT50 DT90 DT50_k1 DT50_k2 +parent 1.887 21.25 1.508 38.83 + +Data: + time variable observed predicted residual + 0 parent 85.1 85.003 0.09726 + 1 parent 57.9 58.039 -0.13912 + 3 parent 29.9 30.054 -0.15351 + 7 parent 14.6 13.866 0.73388 + 14 parent 9.7 9.787 -0.08657 + 28 parent 6.6 7.532 -0.93205 + 63 parent 4.0 4.033 -0.03269 + 91 parent 3.9 2.447 1.45348 + 119 parent 0.6 1.484 -0.88424 |