diff options
| author | Johannes Ranke <jranke@uni-bremen.de> | 2018-01-19 16:01:47 +0100 |
|---|---|---|
| committer | Johannes Ranke <jranke@uni-bremen.de> | 2018-01-30 10:43:58 +0100 |
| commit | a37ba8f8898e4629dfc9d2558fc19a180551de2d (patch) | |
| tree | dcae8090ab31e2e21ebdcb8f8fe11ada521523cc | |
| parent | f18213520f20aba947093e53113c44b689e8b98d (diff) | |
Reweighting with two-component error model
Static documentation except articles rebuilt by pkgdown
| -rw-r--r-- | DESCRIPTION | 2 | ||||
| -rw-r--r-- | NEWS.md | 8 | ||||
| -rw-r--r-- | R/mkinfit.R | 109 | ||||
| -rw-r--r-- | R/sigma_rl.R | 3 | ||||
| -rw-r--r-- | README.html | 4 | ||||
| -rw-r--r-- | README.md | 2 | ||||
| -rw-r--r-- | _pkgdown.yml | 1 | ||||
| -rw-r--r-- | check.log | 12 | ||||
| -rw-r--r-- | docs/index.html | 2 | ||||
| -rw-r--r-- | docs/news/index.html | 10 | ||||
| -rw-r--r-- | docs/reference/index.html | 6 | ||||
| -rw-r--r-- | docs/reference/mkinfit.html | 576 | ||||
| -rw-r--r-- | docs/reference/sigma_rl.html | 167 | ||||
| -rw-r--r-- | man/mkinfit.Rd | 28 | ||||
| -rw-r--r-- | man/sigma_rl.Rd | 25 | ||||
| -rw-r--r-- | test.log | 6 | ||||
| -rw-r--r-- | tests/testthat/test_irls.R | 46 |
17 files changed, 679 insertions, 328 deletions
diff --git a/DESCRIPTION b/DESCRIPTION index 27d81005..e3c7de4e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -2,7 +2,7 @@ Package: mkin Type: Package Title: Kinetic Evaluation of Chemical Degradation Data Version: 0.9.47.1 -Date: 2018-01-05 +Date: 2018-01-30 Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de", comment = c(ORCID = "0000-0003-4371-6538")), @@ -1,3 +1,11 @@ +# mkin 0.9.47.1 (2017-01-30) + +- 'test_data_from_UBA_2014': Added this list of datasets containing experimental data used in the expertise from 2014 + +- 'mkinfit': Added the iterative reweighting method 'tc' using the two-component error model from Rocke and Lorenzato. NA values in the data are not returned any more. + +- 'summary.mkinfit': Improved output regarding weighting method. No predictions are returned for NA values in the model (see above). + # mkin 0.9.46.3 (2017-11-16) - `README.md`, `vignettes/mkin.Rmd`: URLs were updated diff --git a/R/mkinfit.R b/R/mkinfit.R index 55734bac..adc20a9f 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -1,4 +1,4 @@ -# Copyright (C) 2010-2016 Johannes Ranke
+# Copyright (C) 2010-2018 Johannes Ranke
# Portions of this code are copyright (C) 2013 Eurofins Regulatory AG
# Contact: jranke@uni-bremen.de
# The summary function is an adapted and extended version of summary.modFit
@@ -36,7 +36,10 @@ mkinfit <- function(mkinmod, observed, transform_rates = TRUE,
transform_fractions = TRUE,
plot = FALSE, quiet = FALSE,
- err = NULL, weight = "none", scaleVar = FALSE,
+ err = NULL,
+ weight = c("none", "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,
@@ -77,8 +80,9 @@ mkinfit <- function(mkinmod, observed, # Get the names of observed variables
obs_vars <- names(mkinmod$spec)
- # Subset observed data with names of observed data in the model
+ # Subset observed data with names of observed data in the model and remove NA values
observed <- subset(observed, name %in% obs_vars)
+ observed <- subset(observed, !is.na(value))
# Obtain data for decline from maximum mean value if requested
if (from_max_mean) {
@@ -249,11 +253,17 @@ mkinfit <- function(mkinmod, observed, # Define outtimes for model solution.
# Include time points at which observed data are available
- # Make sure we include time 0, so initial values for state variables are for time 0
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_rl(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
@@ -387,33 +397,73 @@ mkinfit <- function(mkinmod, observed, # Reiterate the fit until convergence of the variance components (IRLS)
# if requested by the user
- weight.ini = weight
- if (!is.null(err)) weight.ini = "manual"
if (!is.null(reweight.method)) {
- if (reweight.method != "obs") stop("Only reweighting method 'obs' is implemented")
- if(!quiet) {
- cat("IRLS based on variance estimates for each observed variable\n")
+ if (! reweight.method %in% c("obs", "tc")) stop("Only reweighting methods 'obs' and 'tc' are implemented")
+
+ if (reweight.method == "obs") {
+ 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 (!quiet) {
- cat("Initial variance estimates are:\n")
- print(signif(fit$var_ms_unweighted, 8))
+ if (reweight.method == "tc") {
+ # We need unweighted residuals to update the weighting
+ cost_tmp <- cost(fit$par)
+
+ tc_fit <- nls(abs(res.unweighted) ~ sigma_rl(obs, sigma_low, rsd_high),
+ start = list(sigma_low = tc["sigma_low"], rsd_high = tc["rsd_high"]),
+ data = cost_tmp$residuals,
+ algorithm = "port")
+
+ tc_fitted <- coef(tc_fit)
+ if(!quiet) {
+ cat("IRLS based on variance estimates according to the Rocke and Lorenzato 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) {
n.iter <- n.iter + 1
- sigma.old <- sqrt(fit$var_ms_unweighted)
- observed[err] <- sqrt(fit$var_ms_unweighted)[as.character(observed$name)]
+ # 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
+ observed[err] <- predict(tc_fit)
+
+ }
fit <- modFit(cost, fit$par, method = method.modFit,
control = control.modFit, lower = lower, upper = upper, ...)
- reweight.diff = sum((sqrt(fit$var_ms_unweighted) - sigma.old)^2)
+
+ if (reweight.method == "obs") {
+ sr_new <- fit$var_ms_unweighted
+ }
+ if (reweight.method == "tc") {
+ cost_tmp <- cost(fit$par)
+
+ tc_fit <- nls(abs(res.unweighted) ~ sigma_rl(obs, sigma_low, rsd_high),
+ start = as.list(tc_fitted),
+ data = cost_tmp$residuals,
+ algorithm = "port")
+
+ 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(fit$var_ms_unweighted, 8))
- cat("Sum of squared differences to last variance estimates:",
+ print(signif(sr_new, 8))
+ cat("Sum of squared differences to last variance (component) estimates:",
signif(reweight.diff, 2), "\n")
}
}
@@ -530,9 +580,11 @@ mkinfit <- function(mkinmod, observed, 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
# Return different sets of backtransformed parameters for summary and plotting
fit$bparms.optim <- bparms.optim
@@ -624,6 +676,9 @@ summary.mkinfit <- function(object, data = TRUE, distimes = TRUE, alpha = 0.05, 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),
@@ -679,9 +734,23 @@ print.summary.mkinfit <- function(x, digits = max(3, getOption("digits") - 3), . "using", x$calls, "model solutions performed in", x$time[["elapsed"]], "s\n")
cat("\nWeighting:", x$weight.ini)
- if(!is.null(x$reweight.method)) cat(" then iterative reweighting method",
- x$reweight.method)
- cat("\n")
+ 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("\nStarting values for parameters to be optimised:\n")
print(x$start)
diff --git a/R/sigma_rl.R b/R/sigma_rl.R new file mode 100644 index 00000000..2b921d29 --- /dev/null +++ b/R/sigma_rl.R @@ -0,0 +1,3 @@ + sigma_rl <- function(y, sigma_low, rsd_high) { + sqrt(sigma_low^2 + y^2 * rsd_high^2) + } diff --git a/README.html b/README.html index a38856bd..f4999a38 100644 --- a/README.html +++ b/README.html @@ -121,7 +121,7 @@ $(document).ready(function () { <div id="mkin" class="section level1"> <h1>mkin</h1> -<p><a href="https://cran.r-project.org/package=mkin"><img src="data:image/svg+xml; charset=utf-8;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHdpZHRoPSIxMDAiIGhlaWdodD0iMjAiPgogIDxsaW5lYXJHcmFkaWVudCBpZD0iYiIgeDI9IjAiIHkyPSIxMDAlIj4KICAgIDxzdG9wIG9mZnNldD0iMCIgc3RvcC1jb2xvcj0iI2JiYiIgc3RvcC1vcGFjaXR5PSIuMSIvPgogICAgPHN0b3Agb2Zmc2V0PSIxIiBzdG9wLW9wYWNpdHk9Ii4xIi8+CiAgPC9saW5lYXJHcmFkaWVudD4KICA8bWFzayBpZD0iYSI+CiAgICA8cmVjdCB3aWR0aD0iMTAwIiBoZWlnaHQ9IjIwIiByeD0iMyIgZmlsbD0iI2ZmZiIvPgogIDwvbWFzaz4KICA8ZyBtYXNrPSJ1cmwoI2EpIj4KICAgIDxwYXRoIGZpbGw9IiM1NTUiIGQ9Ik0wIDBoNDN2MjBIMHoiLz4KICAgIDxwYXRoIGZpbGw9IiM0YzEiIGQ9Ik00MyAwaDc5LjV2MjBINDN6Ii8+CiAgICA8cGF0aCBmaWxsPSJ1cmwoI2IpIiBkPSJNMCAwaDEwMHYyMEgweiIvPgogIDwvZz4KICA8ZyBmaWxsPSIjZmZmIiB0ZXh0LWFuY2hvcj0ibWlkZGxlIgogICAgIGZvbnQtZmFtaWx5PSJEZWphVnUgU2FucyxWZXJkYW5hLEdlbmV2YSxzYW5zLXNlcmlmIiBmb250LXNpemU9IjExIj4KICAgIDx0ZXh0IHg9IjIxLjUiIHk9IjE1IiBmaWxsPSIjMDEwMTAxIiBmaWxsLW9wYWNpdHk9Ii4zIj4KICAgICAgQ1JBTgogICAgPC90ZXh0PgogICAgPHRleHQgeD0iMjEuNSIgeT0iMTQiPgogICAgICBDUkFOCiAgICA8L3RleHQ+CiAgICA8dGV4dCB4PSI3MC41IiB5PSIxNSIgZmlsbD0iIzAxMDEwMSIgZmlsbC1vcGFjaXR5PSIuMyI+CiAgICAgIDAuOS40Ni4yCiAgICA8L3RleHQ+CiAgICA8dGV4dCB4PSI3MC41IiB5PSIxNCI+CiAgICAgIDAuOS40Ni4yCiAgICA8L3RleHQ+CiAgPC9nPgo8L3N2Zz4=" /></a></p> +<p><a href="https://cran.r-project.org/package=mkin"><img src="data:image/svg+xml; charset=utf-8;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHdpZHRoPSIxMDAiIGhlaWdodD0iMjAiPgogIDxsaW5lYXJHcmFkaWVudCBpZD0iYiIgeDI9IjAiIHkyPSIxMDAlIj4KICAgIDxzdG9wIG9mZnNldD0iMCIgc3RvcC1jb2xvcj0iI2JiYiIgc3RvcC1vcGFjaXR5PSIuMSIvPgogICAgPHN0b3Agb2Zmc2V0PSIxIiBzdG9wLW9wYWNpdHk9Ii4xIi8+CiAgPC9saW5lYXJHcmFkaWVudD4KICA8bWFzayBpZD0iYSI+CiAgICA8cmVjdCB3aWR0aD0iMTAwIiBoZWlnaHQ9IjIwIiByeD0iMyIgZmlsbD0iI2ZmZiIvPgogIDwvbWFzaz4KICA8ZyBtYXNrPSJ1cmwoI2EpIj4KICAgIDxwYXRoIGZpbGw9IiM1NTUiIGQ9Ik0wIDBoNDN2MjBIMHoiLz4KICAgIDxwYXRoIGZpbGw9IiM0YzEiIGQ9Ik00MyAwaDc5LjV2MjBINDN6Ii8+CiAgICA8cGF0aCBmaWxsPSJ1cmwoI2IpIiBkPSJNMCAwaDEwMHYyMEgweiIvPgogIDwvZz4KICA8ZyBmaWxsPSIjZmZmIiB0ZXh0LWFuY2hvcj0ibWlkZGxlIgogICAgIGZvbnQtZmFtaWx5PSJEZWphVnUgU2FucyxWZXJkYW5hLEdlbmV2YSxzYW5zLXNlcmlmIiBmb250LXNpemU9IjExIj4KICAgIDx0ZXh0IHg9IjIxLjUiIHk9IjE1IiBmaWxsPSIjMDEwMTAxIiBmaWxsLW9wYWNpdHk9Ii4zIj4KICAgICAgQ1JBTgogICAgPC90ZXh0PgogICAgPHRleHQgeD0iMjEuNSIgeT0iMTQiPgogICAgICBDUkFOCiAgICA8L3RleHQ+CiAgICA8dGV4dCB4PSI3MC41IiB5PSIxNSIgZmlsbD0iIzAxMDEwMSIgZmlsbC1vcGFjaXR5PSIuMyI+CiAgICAgIDAuOS40Ni4zCiAgICA8L3RleHQ+CiAgICA8dGV4dCB4PSI3MC41IiB5PSIxNCI+CiAgICAgIDAuOS40Ni4zCiAgICA8L3RleHQ+CiAgPC9nPgo8L3N2Zz4=" /></a></p> <p>The R package <strong>mkin</strong> provides calculation routines for the analysis of chemical degradation data, including <b>m</b>ulticompartment <b>kin</b>etics as needed for modelling the formation and decline of transformation products, or if several compartments are involved.</p> <div id="installation" class="section level2"> <h2>Installation</h2> @@ -162,7 +162,7 @@ $(document).ready(function () { </div> <div id="news" class="section level2"> <h2>News</h2> -<p>Yes, there is a ChangeLog, for the latest <a href="https://cran.r-project.org/package=mkin/news.html">CRAN release</a> and one for the <a href="https://github.com/jranke/mkin/blob/master/NEWS.md">github master branch</a>.</p> +<p>There is a ChangeLog, for the latest <a href="https://cran.r-project.org/package=mkin/news.html">CRAN release</a> and one for the <a href="https://github.com/jranke/mkin/blob/master/NEWS.md">github master branch</a>.</p> </div> <div id="credits-and-historical-remarks" class="section level2"> <h2>Credits and historical remarks</h2> @@ -100,7 +100,7 @@ for installation instructions and a manual. ## News -Yes, there is a ChangeLog, for the latest [CRAN release](https://cran.r-project.org/package=mkin/news.html) +There is a ChangeLog, for the latest [CRAN release](https://cran.r-project.org/package=mkin/news.html) and one for the [github master branch](https://github.com/jranke/mkin/blob/master/NEWS.md). ## Credits and historical remarks diff --git a/_pkgdown.yml b/_pkgdown.yml index f7d868ef..2630bca0 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -48,6 +48,7 @@ reference: - print.mkinmod - transform_odeparms - ilr + - sigma_rl - title: Analytical solutions desc: Parent only model solutions contents: @@ -52,12 +52,10 @@ Maintainer: ‘Johannes Ranke <jranke@uni-bremen.de>’ * checking contents of ‘data’ directory ... OK * checking data for non-ASCII characters ... OK * checking data for ASCII and uncompressed saves ... OK +* checking sizes of PDF files under ‘inst/doc’ ... OK * checking installed files from ‘inst/doc’ ... OK * checking files in ‘vignettes’ ... OK -* checking examples ... NOTE -Examples with CPU or elapsed time > 5s - user system elapsed -test_data_from_UBA_2014 5.736 1.888 5.427 +* checking examples ... OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... SKIPPED * checking for unstated dependencies in vignettes ... OK @@ -66,9 +64,5 @@ test_data_from_UBA_2014 5.736 1.888 5.427 * checking PDF version of manual ... OK * DONE -Status: 1 NOTE -See - ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’ -for details. - +Status: OK diff --git a/docs/index.html b/docs/index.html index f228bbe2..1eefd772 100644 --- a/docs/index.html +++ b/docs/index.html @@ -124,7 +124,7 @@ <div id="news" class="section level2"> <h2 class="hasAnchor"> <a href="#news" class="anchor"></a>News</h2> -<p>Yes, there is a ChangeLog, for the latest <a href="https://cran.r-project.org/package=mkin/news.html">CRAN release</a> and one for the <a href="https://github.com/jranke/mkin/blob/master/NEWS.md">github master branch</a>.</p> +<p>There is a ChangeLog, for the latest <a href="https://cran.r-project.org/package=mkin/news.html">CRAN release</a> and one for the <a href="https://github.com/jranke/mkin/blob/master/NEWS.md">github master branch</a>.</p> </div> <div id="credits-and-historical-remarks" class="section level2"> <h2 class="hasAnchor"> diff --git a/docs/news/index.html b/docs/news/index.html index f442d89c..28ec949f 100644 --- a/docs/news/index.html +++ b/docs/news/index.html @@ -103,6 +103,15 @@ </div> <div class="contents"> + <div id="mkin-0-9-47-1-2017-01-30" class="section level1"> +<h1 class="hasAnchor"> +<a href="#mkin-0-9-47-1-2017-01-30" class="anchor"></a>mkin 0.9.47.1 (2017-01-30)</h1> +<ul> +<li><p>‘test_data_from_UBA_2014’: Added this list of datasets containing experimental data used in the expertise from 2014</p></li> +<li><p>‘mkinfit’: Added the iterative reweighting method ‘tc’ using the two-component error model from Rocke and Lorenzato. NA values in the data are not returned any more.</p></li> +<li><p>‘summary.mkinfit’: Improved output regarding weighting method. No predictions are returned for NA values in the model (see above).</p></li> +</ul> +</div> <div id="mkin-0-9-46-3-2017-11-16" class="section level1"> <h1 class="hasAnchor"> <a href="#mkin-0-9-46-3-2017-11-16" class="anchor"></a>mkin 0.9.46.3 (2017-11-16)</h1> @@ -578,6 +587,7 @@ <div id="tocnav"> <h2>Contents</h2> <ul class="nav nav-pills nav-stacked"> + <li><a href="#mkin-0-9-47-1-2017-01-30">0.9.47.1</a></li> <li><a href="#mkin-0-9-46-3-2017-11-16">0.9.46.3</a></li> <li><a href="#mkin-0-9-46-2-2017-10-10">0.9.46.2</a></li> <li><a href="#mkin-0-9-46-1-2017-09-14">0.9.46.1</a></li> diff --git a/docs/reference/index.html b/docs/reference/index.html index 54886c94..4675748c 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -337,6 +337,12 @@ <p><code><a href="ilr.html">ilr</a></code> <code><a href="ilr.html">invilr</a></code> </p> </td> <td><p>Function to perform isometric log-ratio transformation</p></td> + </tr><tr> + <!-- --> + <td> + <p><code><a href="sigma_rl.html">sigma_rl</a></code> </p> + </td> + <td><p>Two component error model of Rocke and Lorenzato</p></td> </tr> </tbody><tbody> <tr> diff --git a/docs/reference/mkinfit.html b/docs/reference/mkinfit.html index cdf7c7ef..0102aecb 100644 --- a/docs/reference/mkinfit.html +++ b/docs/reference/mkinfit.html @@ -70,6 +70,9 @@ <a href="../articles/FOCUS_L.html">Example evaluation of FOCUS Laboratory Data L1 to L3</a> </li> <li> + <a href="../articles/FOCUS_Z.html">Example evaluation of FOCUS Example Dataset Z</a> + </li> + <li> <a href="../articles/compiled_models.html">Performance benefit by using compiled model definitions in mkin</a> </li> <li> @@ -83,12 +86,7 @@ </ul> <ul class="nav navbar-nav navbar-right"> - <li> - <a href="http://github.com/jranke/mkin"> - <span class="fa fa-github fa-lg"></span> - - </a> -</li> + </ul> </div><!--/.nav-collapse --> </div><!--/.container --> @@ -130,7 +128,9 @@ <span class='kw'>control.modFit</span> <span class='kw'>=</span> <span class='fu'>list</span>(), <span class='kw'>transform_rates</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>transform_fractions</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, - <span class='kw'>plot</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>, <span class='kw'>err</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, <span class='kw'>weight</span> <span class='kw'>=</span> <span class='st'>"none"</span>, + <span class='kw'>plot</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>, <span class='kw'>err</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, + <span class='kw'>weight</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='st'>"none"</span>, <span class='st'>"std"</span>, <span class='st'>"mean"</span>, <span class='st'>"tc"</span>), + <span class='kw'>tc</span> <span class='kw'>=</span> <span class='fu'>c</span>(<span class='kw'>sigma_low</span> <span class='kw'>=</span> <span class='fl'>0.5</span>, <span class='kw'>rsd_high</span> <span class='kw'>=</span> <span class='fl'>0.07</span>), <span class='kw'>scaleVar</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>, <span class='kw'>atol</span> <span class='kw'>=</span> <span class='fl'>1e-8</span>, <span class='kw'>rtol</span> <span class='kw'>=</span> <span class='fl'>1e-10</span>, <span class='kw'>n.outtimes</span> <span class='kw'>=</span> <span class='fl'>100</span>, <span class='kw'>reweight.method</span> <span class='kw'>=</span> <span class='kw'>NULL</span>, @@ -293,7 +293,13 @@ <tr> <th>weight</th> <td><p>only if <code>err</code>=<code>NULL</code>: how to weight the residuals, one of "none", - "std", "mean", see details of <code>modCost</code>.</p></td> + "std", "mean", see details of <code>modCost</code>, or "tc" for the + two component error model of Rocke and Lorenzato.</p></td> + </tr> + <tr> + <th>tc</th> + <td><p>The two components of the Rocke and Lorenzato error model as used + for (initial) weighting</p></td> </tr> <tr> <th>scaleVar</th> @@ -321,12 +327,19 @@ <th>reweight.method</th> <td><p>The method used for iteratively reweighting residuals, also known as iteratively reweighted least squares (IRLS). Default is NULL, - the other method implemented is called "obs", meaning that each - observed variable is assumed to have its own variance, this is - estimated from the fit and used for weighting the residuals - in each iteration until convergence of this estimate up to + i.e. no iterative weighting. + The first reweighting method is called "obs", meaning that each + observed variable is assumed to have its own variance. This variance + is estimated from the fit (mean squared residuals) and used for weighting + the residuals in each iteration until convergence of this estimate up to <code>reweight.tol</code> or up to the maximum number of iterations - specified by <code>reweight.max.iter</code>.</p></td> + specified by <code>reweight.max.iter</code>. + The second reweighting method is called "tc" (two-component error model). + When using this method, the two components of the error model according + to Rocke and Lorenzato (1995) are estimated from the fit and the resulting + variances are used for weighting the residuals in each iteration until + convergence of these components or up to the maximum number of iterations + specified.</p></td> </tr> <tr> <th>reweight.tol</th> @@ -373,21 +386,26 @@ numerical ODE solver. In this situation it may help to switch off the internal rate transformation.</p> + <h2 class="hasAnchor" id="source"><a class="anchor" href="#source"></a>Source</h2> + + <p>Rocke, David M. und Lorenzato, Stefan (1995) A two-component model for + measurement error in analytical chemistry. Technometrics 37(2), 176-184.</p> + <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> <pre class="examples"><div class='input'><span class='co'># Use shorthand notation for parent only degradation</span> <span class='no'>fit</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='st'>"FOMC"</span>, <span class='no'>FOCUS_2006_C</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'>summary</span>(<span class='no'>fit</span>)</div><div class='output co'>#> mkin version: 0.9.46 -#> R version: 3.4.1 -#> Date of fit: Sat Jul 29 15:14:18 2017 -#> Date of summary: Sat Jul 29 15:14:18 2017 +<span class='fu'>summary</span>(<span class='no'>fit</span>)</div><div class='output co'>#> mkin version: 0.9.47.1 +#> R version: 3.4.3 +#> Date of fit: Tue Jan 30 10:05:48 2018 +#> Date of summary: Tue Jan 30 10:05:48 2018 #> #> Equations: #> d_parent/dt = - (alpha/beta) * 1/((time/beta) + 1) * parent #> #> Model predictions using solution type analytical #> -#> Fitted with method Port using 64 model solutions performed in 0.141 s +#> Fitted with method Port using 64 model solutions performed in 0.31 s #> #> Weighting: none #> @@ -456,7 +474,7 @@ <span class='kw'>m1</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>))</div><div class='output co'>#> <span class='message'>Successfully compiled differential equation model from auto-generated C code.</span></div><div class='input'><span class='co'># Fit the model to the FOCUS example dataset D using defaults</span> <span class='fu'>print</span>(<span class='fu'>system.time</span>(<span class='no'>fit</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO</span>, <span class='no'>FOCUS_2006_D</span>, <span class='kw'>solution_type</span> <span class='kw'>=</span> <span class='st'>"eigen"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)))</div><div class='output co'>#> user system elapsed -#> 0.888 0.000 0.890 </div><div class='input'><span class='fu'>coef</span>(<span class='no'>fit</span>)</div><div class='output co'>#> parent_0 log_k_parent_sink log_k_parent_m1 log_k_m1_sink +#> 1.196 0.000 1.195 </div><div class='input'><span class='fu'>coef</span>(<span class='no'>fit</span>)</div><div class='output co'>#> parent_0 log_k_parent_sink log_k_parent_m1 log_k_m1_sink #> 99.59848 -3.03822 -2.98030 -5.24750 </div><div class='input'><span class='fu'><a href='endpoints.html'>endpoints</a></span>(<span class='no'>fit</span>)</div><div class='output co'>#> $ff #> parent_sink parent_m1 m1_sink #> 0.485524 0.514476 1.000000 @@ -535,7 +553,7 @@ #> Model cost at call 153 : 371.2134 #> Optimisation by method Port successfully terminated. #> user system elapsed -#> 0.728 0.000 0.729 </div><div class='input'><span class='fu'>coef</span>(<span class='no'>fit.deSolve</span>)</div><div class='output co'>#> parent_0 log_k_parent_sink log_k_parent_m1 log_k_m1_sink +#> 1.008 0.000 1.006 </div><div class='input'><span class='fu'>coef</span>(<span class='no'>fit.deSolve</span>)</div><div class='output co'>#> parent_0 log_k_parent_sink log_k_parent_m1 log_k_m1_sink #> 99.59848 -3.03822 -2.98030 -5.24750 </div><div class='input'><span class='fu'><a href='endpoints.html'>endpoints</a></span>(<span class='no'>fit.deSolve</span>)</div><div class='output co'>#> $ff #> parent_sink parent_m1 m1_sink #> 0.485524 0.514476 1.000000 @@ -576,10 +594,10 @@ <span class='co'># Weighted fits, including IRLS</span> <span class='no'>SFO_SFO.ff</span> <span class='kw'><-</span> <span class='fu'><a href='mkinmod.html'>mkinmod</a></span>(<span class='kw'>parent</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>, <span class='st'>"m1"</span>), <span class='kw'>m1</span> <span class='kw'>=</span> <span class='fu'><a href='mkinsub.html'>mkinsub</a></span>(<span class='st'>"SFO"</span>), <span class='kw'>use_of_ff</span> <span class='kw'>=</span> <span class='st'>"max"</span>)</div><div class='output co'>#> <span class='message'>Successfully compiled differential equation model from auto-generated C code.</span></div><div class='input'><span class='no'>f.noweight</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>FOCUS_2006_D</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'>summary</span>(<span class='no'>f.noweight</span>)</div><div class='output co'>#> mkin version: 0.9.46 -#> R version: 3.4.1 -#> Date of fit: Sat Jul 29 15:14:28 2017 -#> Date of summary: Sat Jul 29 15:14:28 2017 +<span class='fu'>summary</span>(<span class='no'>f.noweight</span>)</div><div class='output co'>#> mkin version: 0.9.47.1 +#> R version: 3.4.3 +#> Date of fit: Tue Jan 30 10:06:00 2018 +#> Date of summary: Tue Jan 30 10:06:00 2018 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -587,7 +605,7 @@ #> #> Model predictions using solution type deSolve #> -#> Fitted with method Port using 185 model solutions performed in 0.746 s +#> Fitted with method Port using 185 model solutions performed in 0.739 s #> #> Weighting: none #> @@ -653,54 +671,50 @@ #> #> Data: #> time variable observed predicted residual -#> 0 parent 99.46 9.960e+01 -1.385e-01 -#> 0 parent 102.04 9.960e+01 2.442e+00 -#> 1 parent 93.50 9.024e+01 3.262e+00 -#> 1 parent 92.50 9.024e+01 2.262e+00 -#> 3 parent 63.23 7.407e+01 -1.084e+01 -#> 3 parent 68.99 7.407e+01 -5.083e+00 -#> 7 parent 52.32 4.991e+01 2.408e+00 -#> 7 parent 55.13 4.991e+01 5.218e+00 -#> 14 parent 27.27 2.501e+01 2.257e+00 -#> 14 parent 26.64 2.501e+01 1.627e+00 -#> 21 parent 11.50 1.253e+01 -1.035e+00 -#> 21 parent 11.64 1.253e+01 -8.946e-01 -#> 35 parent 2.85 3.148e+00 -2.979e-01 -#> 35 parent 2.91 3.148e+00 -2.379e-01 -#> 50 parent 0.69 7.162e-01 -2.624e-02 -#> 50 parent 0.63 7.162e-01 -8.624e-02 -#> 75 parent 0.05 6.074e-02 -1.074e-02 -#> 75 parent 0.06 6.074e-02 -7.381e-04 -#> 100 parent NA 5.151e-03 NA -#> 100 parent NA 5.151e-03 NA -#> 120 parent NA 7.155e-04 NA -#> 120 parent NA 7.155e-04 NA -#> 0 m1 0.00 0.000e+00 0.000e+00 -#> 0 m1 0.00 0.000e+00 0.000e+00 -#> 1 m1 4.84 4.803e+00 3.704e-02 -#> 1 m1 5.64 4.803e+00 8.370e-01 -#> 3 m1 12.91 1.302e+01 -1.140e-01 -#> 3 m1 12.96 1.302e+01 -6.400e-02 -#> 7 m1 22.97 2.504e+01 -2.075e+00 -#> 7 m1 24.47 2.504e+01 -5.748e-01 -#> 14 m1 41.69 3.669e+01 5.000e+00 -#> 14 m1 33.21 3.669e+01 -3.480e+00 -#> 21 m1 44.37 4.165e+01 2.717e+00 -#> 21 m1 46.44 4.165e+01 4.787e+00 -#> 35 m1 41.22 4.331e+01 -2.093e+00 -#> 35 m1 37.95 4.331e+01 -5.363e+00 -#> 50 m1 41.19 4.122e+01 -2.831e-02 -#> 50 m1 40.01 4.122e+01 -1.208e+00 -#> 75 m1 40.09 3.645e+01 3.643e+00 -#> 75 m1 33.85 3.645e+01 -2.597e+00 -#> 100 m1 31.04 3.198e+01 -9.416e-01 -#> 100 m1 33.13 3.198e+01 1.148e+00 -#> 120 m1 25.15 2.879e+01 -3.640e+00 -#> 120 m1 33.31 2.879e+01 4.520e+00</div><div class='input'><span class='no'>f.irls</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>FOCUS_2006_D</span>, <span class='kw'>reweight.method</span> <span class='kw'>=</span> <span class='st'>"obs"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'>summary</span>(<span class='no'>f.irls</span>)</div><div class='output co'>#> mkin version: 0.9.46 -#> R version: 3.4.1 -#> Date of fit: Sat Jul 29 15:14:30 2017 -#> Date of summary: Sat Jul 29 15:14:30 2017 +#> 0 parent 99.46 99.59848 -1.385e-01 +#> 0 parent 102.04 99.59848 2.442e+00 +#> 1 parent 93.50 90.23787 3.262e+00 +#> 1 parent 92.50 90.23787 2.262e+00 +#> 3 parent 63.23 74.07319 -1.084e+01 +#> 3 parent 68.99 74.07319 -5.083e+00 +#> 7 parent 52.32 49.91206 2.408e+00 +#> 7 parent 55.13 49.91206 5.218e+00 +#> 14 parent 27.27 25.01257 2.257e+00 +#> 14 parent 26.64 25.01257 1.627e+00 +#> 21 parent 11.50 12.53462 -1.035e+00 +#> 21 parent 11.64 12.53462 -8.946e-01 +#> 35 parent 2.85 3.14787 -2.979e-01 +#> 35 parent 2.91 3.14787 -2.379e-01 +#> 50 parent 0.69 0.71624 -2.624e-02 +#> 50 parent 0.63 0.71624 -8.624e-02 +#> 75 parent 0.05 0.06074 -1.074e-02 +#> 75 parent 0.06 0.06074 -7.381e-04 +#> 0 m1 0.00 0.00000 0.000e+00 +#> 0 m1 0.00 0.00000 0.000e+00 +#> 1 m1 4.84 4.80296 3.704e-02 +#> 1 m1 5.64 4.80296 8.370e-01 +#> 3 m1 12.91 13.02400 -1.140e-01 +#> 3 m1 12.96 13.02400 -6.400e-02 +#> 7 m1 22.97 25.04476 -2.075e+00 +#> 7 m1 24.47 25.04476 -5.748e-01 +#> 14 m1 41.69 36.69002 5.000e+00 +#> 14 m1 33.21 36.69002 -3.480e+00 +#> 21 m1 44.37 41.65310 2.717e+00 +#> 21 m1 46.44 41.65310 4.787e+00 +#> 35 m1 41.22 43.31312 -2.093e+00 +#> 35 m1 37.95 43.31312 -5.363e+00 +#> 50 m1 41.19 41.21831 -2.831e-02 +#> 50 m1 40.01 41.21831 -1.208e+00 +#> 75 m1 40.09 36.44703 3.643e+00 +#> 75 m1 33.85 36.44703 -2.597e+00 +#> 100 m1 31.04 31.98163 -9.416e-01 +#> 100 m1 33.13 31.98163 1.148e+00 +#> 120 m1 25.15 28.78984 -3.640e+00 +#> 120 m1 33.31 28.78984 4.520e+00</div><div class='input'><span class='no'>f.irls</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>FOCUS_2006_D</span>, <span class='kw'>reweight.method</span> <span class='kw'>=</span> <span class='st'>"obs"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) +<span class='fu'>summary</span>(<span class='no'>f.irls</span>)</div><div class='output co'>#> mkin version: 0.9.47.1 +#> R version: 3.4.3 +#> Date of fit: Tue Jan 30 10:06:02 2018 +#> Date of summary: Tue Jan 30 10:06:02 2018 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -708,9 +722,14 @@ #> #> Model predictions using solution type deSolve #> -#> Fitted with method Port using 468 model solutions performed in 1.914 s +#> Fitted with method Port using 523 model solutions performed in 2.151 s +#> +#> Weighting: none #> -#> Weighting: none then iterative reweighting method obs +#> Iterative reweighting with method obs +#> Final mean squared residuals of observed variables: +#> parent m1 +#> 11.573408 7.407845 #> #> Starting values for parameters to be optimised: #> value type @@ -739,10 +758,10 @@ #> #> Parameter correlation: #> parent_0 log_k_parent log_k_m1 f_parent_ilr_1 -#> parent_0 1.0000 0.5083 -0.1979 -0.6147 +#> parent_0 1.0000 0.5083 -0.1979 -0.6148 #> log_k_parent 0.5083 1.0000 -0.3894 -0.6062 #> log_k_m1 -0.1979 -0.3894 1.0000 0.7417 -#> f_parent_ilr_1 -0.6147 -0.6062 0.7417 1.0000 +#> f_parent_ilr_1 -0.6148 -0.6062 0.7417 1.0000 #> #> Residual standard error: 1.054 on 36 degrees of freedom #> @@ -751,9 +770,9 @@ #> t-test (unrealistically) based on the assumption of normal distribution #> for estimators of untransformed parameters. #> Estimate t value Pr(>t) Lower Upper -#> parent_0 99.67000 55.630 8.181e-37 96.040000 1.033e+02 +#> parent_0 99.67000 55.630 8.184e-37 96.040000 1.033e+02 #> k_parent 0.09906 21.930 1.016e-22 0.090310 1.087e-01 -#> k_m1 0.00524 7.996 8.487e-10 0.004066 6.753e-03 +#> k_m1 0.00524 7.996 8.486e-10 0.004066 6.753e-03 #> f_parent_to_m1 0.51340 23.000 2.038e-23 0.468100 5.584e-01 #> #> Chi2 error levels in percent: @@ -774,54 +793,50 @@ #> #> Data: #> time variable observed predicted residual err -#> 0 parent 99.46 9.967e+01 -2.122e-01 3.402 -#> 0 parent 102.04 9.967e+01 2.368e+00 3.402 -#> 1 parent 93.50 9.027e+01 3.228e+00 3.402 -#> 1 parent 92.50 9.027e+01 2.228e+00 3.402 -#> 3 parent 63.23 7.405e+01 -1.082e+01 3.402 -#> 3 parent 68.99 7.405e+01 -5.056e+00 3.402 -#> 7 parent 52.32 4.982e+01 2.499e+00 3.402 -#> 7 parent 55.13 4.982e+01 5.309e+00 3.402 -#> 14 parent 27.27 2.490e+01 2.367e+00 3.402 -#> 14 parent 26.64 2.490e+01 1.737e+00 3.402 -#> 21 parent 11.50 1.245e+01 -9.476e-01 3.402 -#> 21 parent 11.64 1.245e+01 -8.076e-01 3.402 -#> 35 parent 2.85 3.110e+00 -2.600e-01 3.402 -#> 35 parent 2.91 3.110e+00 -2.000e-01 3.402 -#> 50 parent 0.69 7.037e-01 -1.374e-02 3.402 -#> 50 parent 0.63 7.037e-01 -7.374e-02 3.402 -#> 75 parent 0.05 5.913e-02 -9.133e-03 3.402 -#> 75 parent 0.06 5.913e-02 8.666e-04 3.402 -#> 100 parent NA 4.969e-03 NA 3.402 -#> 100 parent NA 4.969e-03 NA 3.402 -#> 120 parent NA 6.852e-04 NA 3.402 -#> 120 parent NA 6.852e-04 NA 3.402 -#> 0 m1 0.00 0.000e+00 0.000e+00 2.722 -#> 0 m1 0.00 0.000e+00 0.000e+00 2.722 -#> 1 m1 4.84 4.813e+00 2.672e-02 2.722 -#> 1 m1 5.64 4.813e+00 8.267e-01 2.722 -#> 3 m1 12.91 1.305e+01 -1.378e-01 2.722 -#> 3 m1 12.96 1.305e+01 -8.779e-02 2.722 -#> 7 m1 22.97 2.508e+01 -2.106e+00 2.722 -#> 7 m1 24.47 2.508e+01 -6.062e-01 2.722 -#> 14 m1 41.69 3.671e+01 4.983e+00 2.722 -#> 14 m1 33.21 3.671e+01 -3.497e+00 2.722 -#> 21 m1 44.37 4.165e+01 2.720e+00 2.722 -#> 21 m1 46.44 4.165e+01 4.790e+00 2.722 -#> 35 m1 41.22 4.329e+01 -2.069e+00 2.722 -#> 35 m1 37.95 4.329e+01 -5.339e+00 2.722 -#> 50 m1 41.19 4.119e+01 -3.376e-03 2.722 -#> 50 m1 40.01 4.119e+01 -1.183e+00 2.722 -#> 75 m1 40.09 3.644e+01 3.652e+00 2.722 -#> 75 m1 33.85 3.644e+01 -2.588e+00 2.722 -#> 100 m1 31.04 3.199e+01 -9.497e-01 2.722 -#> 100 m1 33.13 3.199e+01 1.140e+00 2.722 -#> 120 m1 25.15 2.881e+01 -3.659e+00 2.722 -#> 120 m1 33.31 2.881e+01 4.501e+00 2.722</div><div class='input'><span class='no'>f.w.mean</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>FOCUS_2006_D</span>, <span class='kw'>weight</span> <span class='kw'>=</span> <span class='st'>"mean"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'>summary</span>(<span class='no'>f.w.mean</span>)</div><div class='output co'>#> mkin version: 0.9.46 -#> R version: 3.4.1 -#> Date of fit: Sat Jul 29 15:14:31 2017 -#> Date of summary: Sat Jul 29 15:14:31 2017 +#> 0 parent 99.46 99.67218 -2.122e-01 3.402 +#> 0 parent 102.04 99.67218 2.368e+00 3.402 +#> 1 parent 93.50 90.27153 3.228e+00 3.402 +#> 1 parent 92.50 90.27153 2.228e+00 3.402 +#> 3 parent 63.23 74.04648 -1.082e+01 3.402 +#> 3 parent 68.99 74.04648 -5.056e+00 3.402 +#> 7 parent 52.32 49.82092 2.499e+00 3.402 +#> 7 parent 55.13 49.82092 5.309e+00 3.402 +#> 14 parent 27.27 24.90287 2.367e+00 3.402 +#> 14 parent 26.64 24.90287 1.737e+00 3.402 +#> 21 parent 11.50 12.44764 -9.476e-01 3.402 +#> 21 parent 11.64 12.44764 -8.076e-01 3.402 +#> 35 parent 2.85 3.11002 -2.600e-01 3.402 +#> 35 parent 2.91 3.11002 -2.000e-01 3.402 +#> 50 parent 0.69 0.70374 -1.374e-02 3.402 +#> 50 parent 0.63 0.70374 -7.374e-02 3.402 +#> 75 parent 0.05 0.05913 -9.134e-03 3.402 +#> 75 parent 0.06 0.05913 8.662e-04 3.402 +#> 0 m1 0.00 0.00000 0.000e+00 2.722 +#> 0 m1 0.00 0.00000 0.000e+00 2.722 +#> 1 m1 4.84 4.81328 2.672e-02 2.722 +#> 1 m1 5.64 4.81328 8.267e-01 2.722 +#> 3 m1 12.91 13.04779 -1.378e-01 2.722 +#> 3 m1 12.96 13.04779 -8.779e-02 2.722 +#> 7 m1 22.97 25.07615 -2.106e+00 2.722 +#> 7 m1 24.47 25.07615 -6.062e-01 2.722 +#> 14 m1 41.69 36.70729 4.983e+00 2.722 +#> 14 m1 33.21 36.70729 -3.497e+00 2.722 +#> 21 m1 44.37 41.65050 2.720e+00 2.722 +#> 21 m1 46.44 41.65050 4.790e+00 2.722 +#> 35 m1 41.22 43.28866 -2.069e+00 2.722 +#> 35 m1 37.95 43.28866 -5.339e+00 2.722 +#> 50 m1 41.19 41.19338 -3.383e-03 2.722 +#> 50 m1 40.01 41.19338 -1.183e+00 2.722 +#> 75 m1 40.09 36.43820 3.652e+00 2.722 +#> 75 m1 33.85 36.43820 -2.588e+00 2.722 +#> 100 m1 31.04 31.98971 -9.497e-01 2.722 +#> 100 m1 33.13 31.98971 1.140e+00 2.722 +#> 120 m1 25.15 28.80897 -3.659e+00 2.722 +#> 120 m1 33.31 28.80897 4.501e+00 2.722</div><div class='input'><span class='no'>f.w.mean</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>FOCUS_2006_D</span>, <span class='kw'>weight</span> <span class='kw'>=</span> <span class='st'>"mean"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) +<span class='fu'>summary</span>(<span class='no'>f.w.mean</span>)</div><div class='output co'>#> mkin version: 0.9.47.1 +#> R version: 3.4.3 +#> Date of fit: Tue Jan 30 10:06:03 2018 +#> Date of summary: Tue Jan 30 10:06:03 2018 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -829,7 +844,7 @@ #> #> Model predictions using solution type deSolve #> -#> Fitted with method Port using 155 model solutions performed in 0.692 s +#> Fitted with method Port using 155 model solutions performed in 0.675 s #> #> Weighting: mean #> @@ -895,55 +910,51 @@ #> #> Data: #> time variable observed predicted residual -#> 0 parent 99.46 99.730570 -0.270570 -#> 0 parent 102.04 99.730570 2.309430 -#> 1 parent 93.50 90.298055 3.201945 -#> 1 parent 92.50 90.298055 2.201945 -#> 3 parent 63.23 74.025028 -10.795028 -#> 3 parent 68.99 74.025028 -5.035028 -#> 7 parent 52.32 49.748382 2.571618 -#> 7 parent 55.13 49.748382 5.381618 -#> 14 parent 27.27 24.815876 2.454124 -#> 14 parent 26.64 24.815876 1.824124 -#> 21 parent 11.50 12.378849 -0.878849 -#> 21 parent 11.64 12.378849 -0.738849 -#> 35 parent 2.85 3.080219 -0.230219 -#> 35 parent 2.91 3.080219 -0.170219 -#> 50 parent 0.69 0.693958 -0.003958 -#> 50 parent 0.63 0.693958 -0.063958 -#> 75 parent 0.05 0.057888 -0.007888 -#> 75 parent 0.06 0.057888 0.002112 -#> 100 parent NA 0.004829 NA -#> 100 parent NA 0.004829 NA -#> 120 parent NA 0.000662 NA -#> 120 parent NA 0.000662 NA -#> 0 m1 0.00 0.000000 0.000000 -#> 0 m1 0.00 0.000000 0.000000 -#> 1 m1 4.84 4.821488 0.018512 -#> 1 m1 5.64 4.821488 0.818512 -#> 3 m1 12.91 13.066692 -0.156692 -#> 3 m1 12.96 13.066692 -0.106692 -#> 7 m1 22.97 25.101058 -2.131058 -#> 7 m1 24.47 25.101058 -0.631058 -#> 14 m1 41.69 36.720923 4.969077 -#> 14 m1 33.21 36.720923 -3.510923 -#> 21 m1 44.37 41.648353 2.721647 -#> 21 m1 46.44 41.648353 4.791647 -#> 35 m1 41.22 43.269225 -2.049225 -#> 35 m1 37.95 43.269225 -5.319225 -#> 50 m1 41.19 41.173639 0.016361 -#> 50 m1 40.01 41.173639 -1.163639 -#> 75 m1 40.09 36.431224 3.658776 -#> 75 m1 33.85 36.431224 -2.581224 -#> 100 m1 31.04 31.996124 -0.956124 -#> 100 m1 33.13 31.996124 1.133876 -#> 120 m1 25.15 28.824128 -3.674128 -#> 120 m1 33.31 28.824128 4.485872</div><div class='input'><span class='no'>f.w.value</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='fu'>subset</span>(<span class='no'>FOCUS_2006_D</span>, <span class='no'>value</span> <span class='kw'>!=</span> <span class='fl'>0</span>), <span class='kw'>err</span> <span class='kw'>=</span> <span class='st'>"value"</span>, +#> 0 parent 99.46 99.73057 -0.270570 +#> 0 parent 102.04 99.73057 2.309430 +#> 1 parent 93.50 90.29805 3.201945 +#> 1 parent 92.50 90.29805 2.201945 +#> 3 parent 63.23 74.02503 -10.795028 +#> 3 parent 68.99 74.02503 -5.035028 +#> 7 parent 52.32 49.74838 2.571618 +#> 7 parent 55.13 49.74838 5.381618 +#> 14 parent 27.27 24.81588 2.454124 +#> 14 parent 26.64 24.81588 1.824124 +#> 21 parent 11.50 12.37885 -0.878849 +#> 21 parent 11.64 12.37885 -0.738849 +#> 35 parent 2.85 3.08022 -0.230219 +#> 35 parent 2.91 3.08022 -0.170219 +#> 50 parent 0.69 0.69396 -0.003958 +#> 50 parent 0.63 0.69396 -0.063958 +#> 75 parent 0.05 0.05789 -0.007888 +#> 75 parent 0.06 0.05789 0.002112 +#> 0 m1 0.00 0.00000 0.000000 +#> 0 m1 0.00 0.00000 0.000000 +#> 1 m1 4.84 4.82149 0.018512 +#> 1 m1 5.64 4.82149 0.818512 +#> 3 m1 12.91 13.06669 -0.156692 +#> 3 m1 12.96 13.06669 -0.106692 +#> 7 m1 22.97 25.10106 -2.131058 +#> 7 m1 24.47 25.10106 -0.631058 +#> 14 m1 41.69 36.72092 4.969077 +#> 14 m1 33.21 36.72092 -3.510923 +#> 21 m1 44.37 41.64835 2.721647 +#> 21 m1 46.44 41.64835 4.791647 +#> 35 m1 41.22 43.26923 -2.049225 +#> 35 m1 37.95 43.26923 -5.319225 +#> 50 m1 41.19 41.17364 0.016361 +#> 50 m1 40.01 41.17364 -1.163639 +#> 75 m1 40.09 36.43122 3.658776 +#> 75 m1 33.85 36.43122 -2.581224 +#> 100 m1 31.04 31.99612 -0.956124 +#> 100 m1 33.13 31.99612 1.133876 +#> 120 m1 25.15 28.82413 -3.674128 +#> 120 m1 33.31 28.82413 4.485872</div><div class='input'><span class='no'>f.w.value</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='fu'>subset</span>(<span class='no'>FOCUS_2006_D</span>, <span class='no'>value</span> <span class='kw'>!=</span> <span class='fl'>0</span>), <span class='kw'>err</span> <span class='kw'>=</span> <span class='st'>"value"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'>summary</span>(<span class='no'>f.w.value</span>)</div><div class='output co'>#> mkin version: 0.9.46 -#> R version: 3.4.1 -#> Date of fit: Sat Jul 29 15:14:32 2017 -#> Date of summary: Sat Jul 29 15:14:32 2017 +<span class='fu'>summary</span>(<span class='no'>f.w.value</span>)</div><div class='output co'>#> mkin version: 0.9.47.1 +#> R version: 3.4.3 +#> Date of fit: Tue Jan 30 10:06:04 2018 +#> Date of summary: Tue Jan 30 10:06:04 2018 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -951,7 +962,7 @@ #> #> Model predictions using solution type deSolve #> -#> Fitted with method Port using 174 model solutions performed in 0.701 s +#> Fitted with method Port using 174 model solutions performed in 0.68 s #> #> Weighting: manual #> @@ -1062,10 +1073,10 @@ <span class='no'>errors</span> <span class='kw'><-</span> <span class='fu'>c</span>(<span class='kw'>parent</span> <span class='kw'>=</span> <span class='fl'>2</span>, <span class='kw'>m1</span> <span class='kw'>=</span> <span class='fl'>1</span>) <span class='no'>dw</span>$<span class='no'>err.man</span> <span class='kw'><-</span> <span class='no'>errors</span>[<span class='no'>FOCUS_2006_D</span>$<span class='no'>name</span>] <span class='no'>f.w.man</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>dw</span>, <span class='kw'>err</span> <span class='kw'>=</span> <span class='st'>"err.man"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>) -<span class='fu'>summary</span>(<span class='no'>f.w.man</span>)</div><div class='output co'>#> mkin version: 0.9.46 -#> R version: 3.4.1 -#> Date of fit: Sat Jul 29 15:14:33 2017 -#> Date of summary: Sat Jul 29 15:14:33 2017 +<span class='fu'>summary</span>(<span class='no'>f.w.man</span>)</div><div class='output co'>#> mkin version: 0.9.47.1 +#> R version: 3.4.3 +#> Date of fit: Tue Jan 30 10:06:05 2018 +#> Date of summary: Tue Jan 30 10:06:05 2018 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -1073,7 +1084,7 @@ #> #> Model predictions using solution type deSolve #> -#> Fitted with method Port using 297 model solutions performed in 1.223 s +#> Fitted with method Port using 297 model solutions performed in 1.178 s #> #> Weighting: manual #> @@ -1139,55 +1150,51 @@ #> #> Data: #> time variable observed predicted residual err -#> 0 parent 99.46 99.485976 -0.025976 1 -#> 0 parent 102.04 99.485976 2.554024 1 -#> 1 parent 93.50 90.186117 3.313883 1 -#> 1 parent 92.50 90.186117 2.313883 1 -#> 3 parent 63.23 74.113162 -10.883162 1 -#> 3 parent 68.99 74.113162 -5.123162 1 -#> 7 parent 52.32 50.050295 2.269705 1 -#> 7 parent 55.13 50.050295 5.079705 1 -#> 14 parent 27.27 25.179750 2.090250 1 -#> 14 parent 26.64 25.179750 1.460250 1 -#> 21 parent 11.50 12.667654 -1.167654 1 -#> 21 parent 11.64 12.667654 -1.027654 1 -#> 35 parent 2.85 3.206164 -0.356164 1 -#> 35 parent 2.91 3.206164 -0.296164 1 -#> 50 parent 0.69 0.735619 -0.045619 1 -#> 50 parent 0.63 0.735619 -0.105619 1 -#> 75 parent 0.05 0.063256 -0.013256 1 -#> 75 parent 0.06 0.063256 -0.003256 1 -#> 100 parent NA 0.005439 NA 1 -#> 100 parent NA 0.005439 NA 1 -#> 120 parent NA 0.000764 NA 1 -#> 120 parent NA 0.000764 NA 1 -#> 0 m1 0.00 0.000000 0.000000 2 -#> 0 m1 0.00 0.000000 0.000000 2 -#> 1 m1 4.84 4.787287 0.052713 2 -#> 1 m1 5.64 4.787287 0.852713 2 -#> 3 m1 12.91 12.987848 -0.077848 2 -#> 3 m1 12.96 12.987848 -0.027848 2 -#> 7 m1 22.97 24.996945 -2.026945 2 -#> 7 m1 24.47 24.996945 -0.526945 2 -#> 14 m1 41.69 36.663527 5.026473 2 -#> 14 m1 33.21 36.663527 -3.453527 2 -#> 21 m1 44.37 41.656813 2.713187 2 -#> 21 m1 46.44 41.656813 4.783187 2 -#> 35 m1 41.22 43.350312 -2.130312 2 -#> 35 m1 37.95 43.350312 -5.400312 2 -#> 50 m1 41.19 41.256365 -0.066365 2 -#> 50 m1 40.01 41.256365 -1.246365 2 -#> 75 m1 40.09 36.460567 3.629433 2 -#> 75 m1 33.85 36.460567 -2.610567 2 -#> 100 m1 31.04 31.969288 -0.929288 2 -#> 100 m1 33.13 31.969288 1.160712 2 -#> 120 m1 25.15 28.760616 -3.610616 2 -#> 120 m1 33.31 28.760616 4.549384 2</div><div class='input'><span class='no'>f.w.man.irls</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>dw</span>, <span class='kw'>err</span> <span class='kw'>=</span> <span class='st'>"err.man"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, +#> 0 parent 99.46 99.48598 -0.025976 1 +#> 0 parent 102.04 99.48598 2.554024 1 +#> 1 parent 93.50 90.18612 3.313883 1 +#> 1 parent 92.50 90.18612 2.313883 1 +#> 3 parent 63.23 74.11316 -10.883162 1 +#> 3 parent 68.99 74.11316 -5.123162 1 +#> 7 parent 52.32 50.05029 2.269705 1 +#> 7 parent 55.13 50.05029 5.079705 1 +#> 14 parent 27.27 25.17975 2.090250 1 +#> 14 parent 26.64 25.17975 1.460250 1 +#> 21 parent 11.50 12.66765 -1.167654 1 +#> 21 parent 11.64 12.66765 -1.027654 1 +#> 35 parent 2.85 3.20616 -0.356164 1 +#> 35 parent 2.91 3.20616 -0.296164 1 +#> 50 parent 0.69 0.73562 -0.045619 1 +#> 50 parent 0.63 0.73562 -0.105619 1 +#> 75 parent 0.05 0.06326 -0.013256 1 +#> 75 parent 0.06 0.06326 -0.003256 1 +#> 0 m1 0.00 0.00000 0.000000 2 +#> 0 m1 0.00 0.00000 0.000000 2 +#> 1 m1 4.84 4.78729 0.052713 2 +#> 1 m1 5.64 4.78729 0.852713 2 +#> 3 m1 12.91 12.98785 -0.077848 2 +#> 3 m1 12.96 12.98785 -0.027848 2 +#> 7 m1 22.97 24.99695 -2.026945 2 +#> 7 m1 24.47 24.99695 -0.526945 2 +#> 14 m1 41.69 36.66353 5.026473 2 +#> 14 m1 33.21 36.66353 -3.453527 2 +#> 21 m1 44.37 41.65681 2.713187 2 +#> 21 m1 46.44 41.65681 4.783187 2 +#> 35 m1 41.22 43.35031 -2.130312 2 +#> 35 m1 37.95 43.35031 -5.400312 2 +#> 50 m1 41.19 41.25637 -0.066365 2 +#> 50 m1 40.01 41.25637 -1.246365 2 +#> 75 m1 40.09 36.46057 3.629433 2 +#> 75 m1 33.85 36.46057 -2.610567 2 +#> 100 m1 31.04 31.96929 -0.929288 2 +#> 100 m1 33.13 31.96929 1.160712 2 +#> 120 m1 25.15 28.76062 -3.610616 2 +#> 120 m1 33.31 28.76062 4.549384 2</div><div class='input'><span class='no'>f.w.man.irls</span> <span class='kw'><-</span> <span class='fu'>mkinfit</span>(<span class='no'>SFO_SFO.ff</span>, <span class='no'>dw</span>, <span class='kw'>err</span> <span class='kw'>=</span> <span class='st'>"err.man"</span>, <span class='kw'>quiet</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>, <span class='kw'>reweight.method</span> <span class='kw'>=</span> <span class='st'>"obs"</span>) -<span class='fu'>summary</span>(<span class='no'>f.w.man.irls</span>)</div><div class='output co'>#> mkin version: 0.9.46 -#> R version: 3.4.1 -#> Date of fit: Sat Jul 29 15:14:36 2017 -#> Date of summary: Sat Jul 29 15:14:36 2017 +<span class='fu'>summary</span>(<span class='no'>f.w.man.irls</span>)</div><div class='output co'>#> mkin version: 0.9.47.1 +#> R version: 3.4.3 +#> Date of fit: Tue Jan 30 10:06:08 2018 +#> Date of summary: Tue Jan 30 10:06:08 2018 #> #> Equations: #> d_parent/dt = - k_parent * parent @@ -1195,9 +1202,14 @@ #> #> Model predictions using solution type deSolve #> -#> Fitted with method Port using 628 model solutions performed in 2.578 s +#> Fitted with method Port using 692 model solutions performed in 2.733 s +#> +#> Weighting: manual #> -#> Weighting: manual then iterative reweighting method obs +#> Iterative reweighting with method obs +#> Final mean squared residuals of observed variables: +#> parent m1 +#> 11.573407 7.407845 #> #> Starting values for parameters to be optimised: #> value type @@ -1226,10 +1238,10 @@ #> #> Parameter correlation: #> parent_0 log_k_parent log_k_m1 f_parent_ilr_1 -#> parent_0 1.0000 0.5083 -0.1979 -0.6147 +#> parent_0 1.0000 0.5083 -0.1979 -0.6148 #> log_k_parent 0.5083 1.0000 -0.3894 -0.6062 #> log_k_m1 -0.1979 -0.3894 1.0000 0.7417 -#> f_parent_ilr_1 -0.6147 -0.6062 0.7417 1.0000 +#> f_parent_ilr_1 -0.6148 -0.6062 0.7417 1.0000 #> #> Residual standard error: 1.054 on 36 degrees of freedom #> @@ -1238,10 +1250,10 @@ #> t-test (unrealistically) based on the assumption of normal distribution #> for estimators of untransformed parameters. #> Estimate t value Pr(>t) Lower Upper -#> parent_0 99.67000 55.630 8.179e-37 96.040000 1.033e+02 -#> k_parent 0.09906 21.930 1.015e-22 0.090310 1.087e-01 -#> k_m1 0.00524 7.996 8.488e-10 0.004066 6.753e-03 -#> f_parent_to_m1 0.51340 23.000 2.039e-23 0.468100 5.584e-01 +#> parent_0 99.67000 55.630 8.184e-37 96.040000 1.033e+02 +#> k_parent 0.09906 21.930 1.016e-22 0.090310 1.087e-01 +#> k_m1 0.00524 7.996 8.486e-10 0.004066 6.753e-03 +#> f_parent_to_m1 0.51340 23.000 2.038e-23 0.468100 5.584e-01 #> #> Chi2 error levels in percent: #> err.min n.optim df @@ -1257,54 +1269,50 @@ #> Estimated disappearance times: #> DT50 DT90 #> parent 6.997 23.24 -#> m1 132.281 439.43 +#> m1 132.282 439.43 #> #> Data: #> time variable observed predicted residual err.ini err -#> 0 parent 99.46 9.967e+01 -2.122e-01 1 3.402 -#> 0 parent 102.04 9.967e+01 2.368e+00 1 3.402 -#> 1 parent 93.50 9.027e+01 3.228e+00 1 3.402 -#> 1 parent 92.50 9.027e+01 2.228e+00 1 3.402 -#> 3 parent 63.23 7.405e+01 -1.082e+01 1 3.402 -#> 3 parent 68.99 7.405e+01 -5.056e+00 1 3.402 -#> 7 parent 52.32 4.982e+01 2.499e+00 1 3.402 -#> 7 parent 55.13 4.982e+01 5.309e+00 1 3.402 -#> 14 parent 27.27 2.490e+01 2.367e+00 1 3.402 -#> 14 parent 26.64 2.490e+01 1.737e+00 1 3.402 -#> 21 parent 11.50 1.245e+01 -9.477e-01 1 3.402 -#> 21 parent 11.64 1.245e+01 -8.077e-01 1 3.402 -#> 35 parent 2.85 3.110e+00 -2.600e-01 1 3.402 -#> 35 parent 2.91 3.110e+00 -2.000e-01 1 3.402 -#> 50 parent 0.69 7.037e-01 -1.375e-02 1 3.402 -#> 50 parent 0.63 7.037e-01 -7.375e-02 1 3.402 -#> 75 parent 0.05 5.913e-02 -9.134e-03 1 3.402 -#> 75 parent 0.06 5.913e-02 8.659e-04 1 3.402 -#> 100 parent NA 4.969e-03 NA 1 3.402 -#> 100 parent NA 4.969e-03 NA 1 3.402 -#> 120 parent NA 6.852e-04 NA 1 3.402 -#> 120 parent NA 6.852e-04 NA 1 3.402 -#> 0 m1 0.00 0.000e+00 0.000e+00 2 2.722 -#> 0 m1 0.00 0.000e+00 0.000e+00 2 2.722 -#> 1 m1 4.84 4.813e+00 2.672e-02 2 2.722 -#> 1 m1 5.64 4.813e+00 8.267e-01 2 2.722 -#> 3 m1 12.91 1.305e+01 -1.378e-01 2 2.722 -#> 3 m1 12.96 1.305e+01 -8.778e-02 2 2.722 -#> 7 m1 22.97 2.508e+01 -2.106e+00 2 2.722 -#> 7 m1 24.47 2.508e+01 -6.061e-01 2 2.722 -#> 14 m1 41.69 3.671e+01 4.983e+00 2 2.722 -#> 14 m1 33.21 3.671e+01 -3.497e+00 2 2.722 -#> 21 m1 44.37 4.165e+01 2.719e+00 2 2.722 -#> 21 m1 46.44 4.165e+01 4.789e+00 2 2.722 -#> 35 m1 41.22 4.329e+01 -2.069e+00 2 2.722 -#> 35 m1 37.95 4.329e+01 -5.339e+00 2 2.722 -#> 50 m1 41.19 4.119e+01 -3.394e-03 2 2.722 -#> 50 m1 40.01 4.119e+01 -1.183e+00 2 2.722 -#> 75 m1 40.09 3.644e+01 3.652e+00 2 2.722 -#> 75 m1 33.85 3.644e+01 -2.588e+00 2 2.722 -#> 100 m1 31.04 3.199e+01 -9.497e-01 2 2.722 -#> 100 m1 33.13 3.199e+01 1.140e+00 2 2.722 -#> 120 m1 25.15 2.881e+01 -3.659e+00 2 2.722 -#> 120 m1 33.31 2.881e+01 4.501e+00 2 2.722</div><div class='input'> +#> 0 parent 99.46 99.67218 -2.122e-01 1 3.402 +#> 0 parent 102.04 99.67218 2.368e+00 1 3.402 +#> 1 parent 93.50 90.27153 3.228e+00 1 3.402 +#> 1 parent 92.50 90.27153 2.228e+00 1 3.402 +#> 3 parent 63.23 74.04648 -1.082e+01 1 3.402 +#> 3 parent 68.99 74.04648 -5.056e+00 1 3.402 +#> 7 parent 52.32 49.82092 2.499e+00 1 3.402 +#> 7 parent 55.13 49.82092 5.309e+00 1 3.402 +#> 14 parent 27.27 24.90288 2.367e+00 1 3.402 +#> 14 parent 26.64 24.90288 1.737e+00 1 3.402 +#> 21 parent 11.50 12.44765 -9.476e-01 1 3.402 +#> 21 parent 11.64 12.44765 -8.076e-01 1 3.402 +#> 35 parent 2.85 3.11002 -2.600e-01 1 3.402 +#> 35 parent 2.91 3.11002 -2.000e-01 1 3.402 +#> 50 parent 0.69 0.70375 -1.375e-02 1 3.402 +#> 50 parent 0.63 0.70375 -7.375e-02 1 3.402 +#> 75 parent 0.05 0.05913 -9.134e-03 1 3.402 +#> 75 parent 0.06 0.05913 8.662e-04 1 3.402 +#> 0 m1 0.00 0.00000 0.000e+00 2 2.722 +#> 0 m1 0.00 0.00000 0.000e+00 2 2.722 +#> 1 m1 4.84 4.81328 2.672e-02 2 2.722 +#> 1 m1 5.64 4.81328 8.267e-01 2 2.722 +#> 3 m1 12.91 13.04779 -1.378e-01 2 2.722 +#> 3 m1 12.96 13.04779 -8.779e-02 2 2.722 +#> 7 m1 22.97 25.07615 -2.106e+00 2 2.722 +#> 7 m1 24.47 25.07615 -6.062e-01 2 2.722 +#> 14 m1 41.69 36.70729 4.983e+00 2 2.722 +#> 14 m1 33.21 36.70729 -3.497e+00 2 2.722 +#> 21 m1 44.37 41.65050 2.720e+00 2 2.722 +#> 21 m1 46.44 41.65050 4.790e+00 2 2.722 +#> 35 m1 41.22 43.28866 -2.069e+00 2 2.722 +#> 35 m1 37.95 43.28866 -5.339e+00 2 2.722 +#> 50 m1 41.19 41.19339 -3.386e-03 2 2.722 +#> 50 m1 40.01 41.19339 -1.183e+00 2 2.722 +#> 75 m1 40.09 36.43820 3.652e+00 2 2.722 +#> 75 m1 33.85 36.43820 -2.588e+00 2 2.722 +#> 100 m1 31.04 31.98971 -9.497e-01 2 2.722 +#> 100 m1 33.13 31.98971 1.140e+00 2 2.722 +#> 120 m1 25.15 28.80897 -3.659e+00 2 2.722 +#> 120 m1 33.31 28.80897 4.501e+00 2 2.722</div><div class='input'> </div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> @@ -1319,6 +1327,8 @@ <li><a href="#note">Note</a></li> <li><a href="#note">Note</a></li> + + <li><a href="#source">Source</a></li> <li><a href="#examples">Examples</a></li> </ul> diff --git a/docs/reference/sigma_rl.html b/docs/reference/sigma_rl.html new file mode 100644 index 00000000..36f2a0c1 --- /dev/null +++ b/docs/reference/sigma_rl.html @@ -0,0 +1,167 @@ +<!-- Generated by pkgdown: do not edit by hand --> +<!DOCTYPE html> +<html> + <head> + <meta charset="utf-8"> +<meta http-equiv="X-UA-Compatible" content="IE=edge"> +<meta name="viewport" content="width=device-width, initial-scale=1.0"> + +<title>Two component error model of Rocke and Lorenzato — sigma_rl • mkin</title> + +<!-- jquery --> +<script src="https://code.jquery.com/jquery-3.1.0.min.js" integrity="sha384-nrOSfDHtoPMzJHjVTdCopGqIqeYETSXhZDFyniQ8ZHcVy08QesyHcnOUpMpqnmWq" crossorigin="anonymous"></script> +<!-- Bootstrap --> + +<link href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css" rel="stylesheet" integrity="sha384-BVYiiSIFeK1dGmJRAkycuHAHRg32OmUcww7on3RYdg4Va+PmSTsz/K68vbdEjh4u" crossorigin="anonymous"> +<script src="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/js/bootstrap.min.js" integrity="sha384-Tc5IQib027qvyjSMfHjOMaLkfuWVxZxUPnCJA7l2mCWNIpG9mGCD8wGNIcPD7Txa" crossorigin="anonymous"></script> + +<!-- Font Awesome icons --> +<link href="https://maxcdn.bootstrapcdn.com/font-awesome/4.6.3/css/font-awesome.min.css" rel="stylesheet" integrity="sha384-T8Gy5hrqNKT+hzMclPo118YTQO6cYprQmhrYwIiQ/3axmI1hQomh7Ud2hPOy8SP1" crossorigin="anonymous"> + + +<!-- pkgdown --> +<link href="../pkgdown.css" rel="stylesheet"> +<script src="../jquery.sticky-kit.min.js"></script> +<script src="../pkgdown.js"></script> + +<!-- mathjax --> +<script src='https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script> + +<!--[if lt IE 9]> +<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script> +<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script> +<![endif]--> + + + </head> + + <body> + <div class="container template-reference-topic"> + <header> + <div class="navbar navbar-default navbar-fixed-top" role="navigation"> + <div class="container"> + <div class="navbar-header"> + <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar"> + <span class="icon-bar"></span> + <span class="icon-bar"></span> + <span class="icon-bar"></span> + </button> + <a class="navbar-brand" href="../index.html">mkin</a> + </div> + <div id="navbar" class="navbar-collapse collapse"> + <ul class="nav navbar-nav"> + <li> + <a href="../reference/index.html">Functions and data</a> +</li> +<li class="dropdown"> + <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false"> + Articles + + <span class="caret"></span> + </a> + <ul class="dropdown-menu" role="menu"> + <li> + <a href="../articles/mkin.html">Introduction to mkin</a> + </li> + <li> + <a href="../articles/FOCUS_D.html">Example evaluation of FOCUS Example Dataset D</a> + </li> + <li> + <a href="../articles/FOCUS_L.html">Example evaluation of FOCUS Laboratory Data L1 to L3</a> + </li> + <li> + <a href="../articles/FOCUS_Z.html">Example evaluation of FOCUS Example Dataset Z</a> + </li> + <li> + <a href="../articles/compiled_models.html">Performance benefit by using compiled model definitions in mkin</a> + </li> + <li> + <a href="../articles/twa.html">Calculation of time weighted average concentrations with mkin</a> + </li> + </ul> +</li> +<li> + <a href="../news/index.html">News</a> +</li> + </ul> + + <ul class="nav navbar-nav navbar-right"> + + </ul> + </div><!--/.nav-collapse --> + </div><!--/.container --> +</div><!--/.navbar --> + + + </header> + + <div class="row"> + <div class="col-md-9 contents"> + <div class="page-header"> + <h1>Two component error model of Rocke and Lorenzato</h1> + </div> + + + <p>Function describing the standard deviation of the measurement error + in dependence of the measured value:</p> +<p>\(sigma = sqrt(sigma_low^2 + y^2 * rsd_high^2)\)</p> + + + <pre class="usage"><span class='fu'>sigma_rl</span>(<span class='no'>y</span>, <span class='no'>sigma_low</span>, <span class='no'>rsd_high</span>)</pre> + + <h2 class="hasAnchor" id="arguments"><a class="anchor" href="#arguments"></a> Arguments</h2> + <table class="ref-arguments"> + <colgroup><col class="name" /><col class="desc" /></colgroup> + <tr> + <th>y</th> + <td><p>The magnitude of the observed value</p></td> + </tr> + <tr> + <th>sigma_low</th> + <td><p>The asymptotic minimum of the standard deviation for low observed values</p></td> + </tr> + <tr> + <th>rsd_high</th> + <td><p>The coefficient describing the increase of the standard deviation with + the magnitude of the observed value</p></td> + </tr> + </table> + + <h2 class="hasAnchor" id="value"><a class="anchor" href="#value"></a>Value</h2> + + <p>The standard deviation of the response variable.</p> + + <h2 class="hasAnchor" id="references"><a class="anchor" href="#references"></a>References</h2> + + <p>Rocke, David M. und Lorenzato, Stefan (1995) A two-component model for + measurement error in analytical chemistry. Technometrics 37(2), 176-184.</p> + + + </div> + <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> + <h2>Contents</h2> + <ul class="nav nav-pills nav-stacked"> + <li><a href="#arguments">Arguments</a></li> + + <li><a href="#value">Value</a></li> + + <li><a href="#references">References</a></li> + </ul> + + </div> +</div> + + <footer> + <div class="copyright"> + <p>Developed by Johannes Ranke.</p> +</div> + +<div class="pkgdown"> + <p>Site built with <a href="http://hadley.github.io/pkgdown/">pkgdown</a>.</p> +</div> + + </footer> + </div> + + </body> +</html> diff --git a/man/mkinfit.Rd b/man/mkinfit.Rd index 8e2fbeb1..32edb28b 100644 --- a/man/mkinfit.Rd +++ b/man/mkinfit.Rd @@ -30,7 +30,9 @@ mkinfit(mkinmod, observed, control.modFit = list(), transform_rates = TRUE, transform_fractions = TRUE, - plot = FALSE, quiet = FALSE, err = NULL, weight = "none", + plot = FALSE, quiet = FALSE, err = NULL, + weight = c("none", "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, @@ -176,8 +178,11 @@ mkinfit(mkinmod, observed, } \item{weight}{ only if \code{err}=\code{NULL}: how to weight the residuals, one of "none", - "std", "mean", see details of \code{\link{modCost}}. + "std", "mean", see details of \code{\link{modCost}}, or "tc" for the + two component error model of Rocke and Lorenzato. } + \item{tc}{The two components of the Rocke and Lorenzato error model as used + for (initial) weighting}. \item{scaleVar}{ Will be passed to \code{\link{modCost}}. Default is not to scale Variables according to the number of observations. @@ -199,12 +204,19 @@ mkinfit(mkinmod, observed, \item{reweight.method}{ The method used for iteratively reweighting residuals, also known as iteratively reweighted least squares (IRLS). Default is NULL, - the other method implemented is called "obs", meaning that each - observed variable is assumed to have its own variance, this is - estimated from the fit and used for weighting the residuals - in each iteration until convergence of this estimate up to + i.e. no iterative weighting. + The first reweighting method is called "obs", meaning that each + observed variable is assumed to have its own variance. This variance + is estimated from the fit (mean squared residuals) and used for weighting + the residuals in each iteration until convergence of this estimate up to \code{reweight.tol} or up to the maximum number of iterations specified by \code{reweight.max.iter}. + The second reweighting method is called "tc" (two-component error model). + When using this method, the two components of the error model according + to Rocke and Lorenzato (1995) are estimated from the fit and the resulting + variances are used for weighting the residuals in each iteration until + convergence of these components or up to the maximum number of iterations + specified. } \item{reweight.tol}{ Tolerance for convergence criterion for the variance components @@ -243,6 +255,10 @@ mkinfit(mkinmod, observed, numerical ODE solver. In this situation it may help to switch off the internal rate transformation. } +\source{ + Rocke, David M. und Lorenzato, Stefan (1995) A two-component model for + measurement error in analytical chemistry. Technometrics 37(2), 176-184. +} \author{ Johannes Ranke } diff --git a/man/sigma_rl.Rd b/man/sigma_rl.Rd new file mode 100644 index 00000000..d1c22a77 --- /dev/null +++ b/man/sigma_rl.Rd @@ -0,0 +1,25 @@ +\name{sigma_rl} +\alias{sigma_rl} +\title{ Two component error model of Rocke and Lorenzato} +\description{ + Function describing the standard deviation of the measurement error + in dependence of the measured value: + + \eqn{sigma = sqrt(sigma_low^2 + y^2 * rsd_high^2)} +} +\usage{ +sigma_rl(y, sigma_low, rsd_high) +} +\arguments{ + \item{y}{ The magnitude of the observed value } + \item{sigma_low}{ The asymptotic minimum of the standard deviation for low observed values } + \item{rsd_high}{ The coefficient describing the increase of the standard deviation with + the magnitude of the observed value } +} +\value{ + The standard deviation of the response variable. +} +\references{ + Rocke, David M. und Lorenzato, Stefan (1995) A two-component model for + measurement error in analytical chemistry. Technometrics 37(2), 176-184. +} @@ -1,14 +1,10 @@ Loading mkin Loading required package: testthat -Loading required package: minpack.lm -Loading required package: rootSolve -Loading required package: inline -Loading required package: methods -Loading required package: parallel Testing mkin Calculation of FOCUS chi2 error levels: .. Results for FOCUS D established in expertise for UBA (Ranke 2014): ...... The t-test for significant difference from zero: .. +Iteratively reweighed least squares (IRLS) fitting: .. Model predictions with mkinpredict: ... Fitting of parent only models: ..................... Complex test case from Schaefer et al. (2007) Piacenza paper: .. diff --git a/tests/testthat/test_irls.R b/tests/testthat/test_irls.R new file mode 100644 index 00000000..0b005182 --- /dev/null +++ b/tests/testthat/test_irls.R @@ -0,0 +1,46 @@ +# Copyright (C) 2018 Johannes Ranke +# Contact: jranke@uni-bremen.de + +# This file is part of the R package mkin + +# mkin is free software: you can redistribute it and/or modify it under the +# terms of the GNU General Public License as published by the Free Software +# Foundation, either version 3 of the License, or (at your option) any later +# version. + +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more +# details. + +# You should have received a copy of the GNU General Public License along with +# this program. If not, see <http://www.gnu.org/licenses/> + +context("Iteratively reweighed least squares (IRLS) fitting") + + +m_synth_SFO_lin <- mkinmod(parent = mkinsub("SFO", "M1"), + M1 = mkinsub("SFO", "M2"), + M2 = mkinsub("SFO"), + use_of_ff = "max", quiet = TRUE) + +m_synth_DFOP_par <- mkinmod(parent = mkinsub("DFOP", c("M1", "M2")), + M1 = mkinsub("SFO"), + M2 = mkinsub("SFO"), + use_of_ff = "max", quiet = TRUE) + +SFO_lin_a <- synthetic_data_for_UBA_2014[[1]]$data + +DFOP_par_c <- synthetic_data_for_UBA_2014[[12]]$data + +test_that("Reweighting method 'obs' works", { + fit_irls_1 <- mkinfit(m_synth_SFO_lin, SFO_lin_a, reweight.method = "obs", quiet = TRUE) + parms_1 <- round(fit_irls_1$bparms.optim, c(1, 4, 4, 4, 4, 4)) + expect_equivalent(parms_1, c(102.1, 0.7389, 0.2982, 0.0203, 0.7677, 0.7246)) +}) + +test_that("Reweighting method 'tc' works", { + fit_irls_2 <- mkinfit(m_synth_DFOP_par, DFOP_par_c, reweight.method = "tc", quiet = TRUE) + parms_2 <- signif(fit_irls_2$bparms.optim, 3) + expect_equivalent(parms_2, c(99.3, 0.041, 0.00962, 0.597, 0.393, 0.298, 0.0203, 0.707)) +}) |