diff options
| author | Johannes Ranke <jranke@uni-bremen.de> | 2018-07-17 17:29:14 +0200 |
|---|---|---|
| committer | Johannes Ranke <jranke@uni-bremen.de> | 2018-07-17 17:38:29 +0200 |
| commit | 280d36230052de4f94e384648c1283031fbc9840 (patch) | |
| tree | df0ba9e07386b593cc396b8b6976210d42ee1a46 | |
| parent | e636c17f0d354a8e74546fc1469431dbe502dc76 (diff) | |
Fix inverse predictions for replicate measurements
For details, see NEWS.md
45 files changed, 1272 insertions, 296 deletions
diff --git a/.Rbuildignore b/.Rbuildignore index c7ccaca..78e037c 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,11 +1,6 @@ GNUmakefile -out$ -toc$ -bbl$ -blg$ -aux$ -log$ -vignettes/chemCal.tex -vignettes/chemCal.pdf -vignettes/chemCal.R +build.log$ chemCal_.*.tar.gz +docs +_pkgdown.yml +README.html diff --git a/DESCRIPTION b/DESCRIPTION index 7610ff0..2e1f741 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,11 @@ Package: chemCal Version: 0.2.1 -Date: 2018-07-05 +Date: 2018-07-17 Title: Calibration Functions for Analytical Chemistry Authors@R: c(person("Johannes", "Ranke", role = c("aut", "cre", "cph"), email = "jranke@uni-bremen.de", comment = c(ORCID = "0000-0003-4371-6538"))) -Suggests: MASS, knitr, testthat +Suggests: MASS, knitr, testthat, investr Description: Simple functions for plotting linear calibration functions and estimating standard errors for measurements according to the Handbook of Chemometrics and Qualimetrics: Part A diff --git a/GNUmakefile b/GNUmakefile index 3f16fbc..7992c5e 100644 --- a/GNUmakefile +++ b/GNUmakefile @@ -8,27 +8,44 @@ TGZ := $(PKGNAME)_$(PKGVERS).tar.gz # If no alternate bin folder is specified, the default is to use the folder # containing the first instance of R on the PATH. RBIN ?= $(shell dirname "`which R`") +# +# Vignettes are listed in the build target +pkgfiles = \ + GNUmakefile \ + .Rbuildignore \ + data/* \ + DESCRIPTION \ + man/* \ + NAMESPACE \ + NEWS.md \ + README.html \ + R/* \ + tests/* \ + tests/testthat* + +all: build + +$(TGZ): $(pkgfiles) vignettes + "$(RBIN)/R" CMD build . 2>&1 | tee build.log + +build: $(TGZ) README.html: README.md "$(RBIN)/Rscript" -e "rmarkdown::render('README.md', output_format = 'html_document')" -build: - "$(RBIN)/R" CMD build . - install: build "$(RBIN)/R" CMD INSTALL $(TGZ) check: build "$(RBIN)/R" CMD check --as-cran $(TGZ) -vignettes/%.html: vignettes/%.Rmd - "$(RBIN)/Rscript" -e "tools::buildVignette(file = 'vignettes/$*.Rnw', dir = 'vignettes')" +vignettes/%.html: vignettes/%.Rmd vignettes/refs.bib + "$(RBIN)/Rscript" -e "tools::buildVignette(file = 'vignettes/$*.Rmd', dir = 'vignettes')" vignettes: vignettes/chemCal.html pd: "$(RBIN)/Rscript" -e "pkgdown::build_site()" - cp vignettes/chemCal.pdf docs/articles git add -A git commit -m 'Static documentation rebuilt by pkgdown::build_site()' -e @@ -40,3 +57,9 @@ winbuilder: build test: install NOT_CRAN=true "$(RBIN)/Rscript" -e 'devtools::test()' + +clean: + $(RM) -r vignettes/*_cache + $(RM) -r vignettes/*_files + $(RM) -r vignettes/*.R + $(RM) Rplots.pdf @@ -1,12 +1,16 @@ -# chemCal 0.2.1 (2018-07-16) +# chemCal 0.2.1 (2018-07-17) -- 'lod' and 'loq': In the lists that are returned, return the list component 'y' without names, because we always have a single element in 'y' and the name '1' is meaningless +- 'inverse.predict': Do not work on the means of the calibration standards any more, as this ignores the variability of y values about the means. Thanks to Anna Burniol Figols for pointing out this issue -- Use testthat for tests to simplify further development +- Use testthat for tests to simplify development. Adapt the tests using data with replicate standard measurements to work on the means in order to show the relation to 'inverse.predict' from earlier versions. Include comparisons with investr::calibrate(method = 'Wald') for unweighted regressions. Include tests with more precision to check for changes in numerical output across versions. -- Convert vignette to html +- 'lod' and 'loq': In the lists that are returned, return the list component 'y' without names, because we always only have a single element in 'y' (previously the name '1' was returned). -- Add another example dataset, from an online course at the University of Toronto +- Convert vignette to html and explain the changes to 'inverse.predict' + +- Add two example dataset, one from an online course at the University of Toronto, one from Rocke and Lorenzato (1995) + +- Update static documentation # chemCal 0.1-33 (2014-04-24) @@ -1,5 +0,0 @@ -- lod(): At the moment, it is not possible to calculate an lod for the - case of more than one replicates (m is fixed to 1). The formula - for the prediction of y from mean(x) in Massart et al, p. 433 - could be used to generalize this. -- Write methods for nonlinear calibration functions diff --git a/_pkgdown.yml b/_pkgdown.yml index edc41e2..f5eeff8 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -2,12 +2,3 @@ title: chemCal template: bootswatch: spacelab - -navbar: - title: ~ - type: default - left: - - text: Functions and data - href: reference/index.html - - text: Short manual (pdf) - href: articles/chemCal.pdf diff --git a/build.log b/build.log new file mode 100644 index 0000000..8e55d71 --- /dev/null +++ b/build.log @@ -0,0 +1,11 @@ +* checking for file ‘./DESCRIPTION’ ... OK +* preparing ‘chemCal’: +* checking DESCRIPTION meta-information ... OK +* installing the package to build vignettes +* creating vignettes ... OK +* checking for LF line-endings in source and make files and shell scripts +* checking for empty or unneeded directories +* looking to see if a ‘data/datalist’ file should be added +* re-saving image files +* building ‘chemCal_0.2.1.tar.gz’ + diff --git a/data/rl95_toluene.rda b/data/rl95_toluene.rda Binary files differnew file mode 100644 index 0000000..4b213be --- /dev/null +++ b/data/rl95_toluene.rda diff --git a/docs/articles/chemCal.html b/docs/articles/chemCal.html new file mode 100644 index 0000000..29db7c8 --- /dev/null +++ b/docs/articles/chemCal.html @@ -0,0 +1,209 @@ +<!DOCTYPE html> +<!-- Generated by pkgdown: do not edit by hand --><html> +<head> +<meta http-equiv="Content-Type" content="text/html; charset=UTF-8"> +<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>Introduction to chemCal • chemCal</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"> +<!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/1.7.1/clipboard.min.js" integrity="sha384-cV+rhyOuRHc9Ub/91rihWcGmMmCXDeksTtCihMupQHSsi8GIIRDG0ThDc3HGQFJ3" crossorigin="anonymous"></script><!-- sticky kit --><script src="https://cdnjs.cloudflare.com/ajax/libs/sticky-kit/1.1.3/sticky-kit.min.js" integrity="sha256-c4Rlo1ZozqTPE2RLuvbusY3+SU1pQaJC0TjuhygMipw=" crossorigin="anonymous"></script><!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet"> +<script src="../pkgdown.js"></script><meta property="og:title" content="Introduction to chemCal"> +<meta property="og:description" content=""> +<meta name="twitter:card" content="summary"> +<!-- 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-article"> + <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> + <span class="navbar-brand"> + <a class="navbar-link" href="../index.html">chemCal</a> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> + </span> + </div> + + <div id="navbar" class="navbar-collapse collapse"> + <ul class="nav navbar-nav"> +<li> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> +</li> +<li> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</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 toc-ignore"> + <h1>Introduction to chemCal</h1> + <h4 class="author">Johannes Ranke</h4> + + <h4 class="date">2018-07-17</h4> + + + <div class="hidden name"><code>chemCal.Rmd</code></div> + + </div> + + + +<div id="basic-calibration-functions" class="section level1"> +<h1 class="hasAnchor"> +<a href="#basic-calibration-functions" class="anchor"></a>Basic calibration functions</h1> +<p>The <code>chemCal</code> package was first designed in the course of a lecture and lab course on “Analytics of Organic Trace Contaminants” at the University of Bremen from October to December 2004. In the fall 2005, an email exchange with Ron Wehrens led to the belief that it would be desirable to implement the inverse prediction method given in <span class="citation">Massart et al. (1997)</span> since it also covers the case of weighted regression. Studies of the IUPAC orange book and of DIN 32645 (equivalent to ISO 11843), publications by <span class="citation">Currie (1997)</span> and the Analytical Method Committee of the Royal Society of Chemistry <span class="citation">(Analytical Methods Committee 1989)</span> and a nice paper by Castells and Castillo <span class="citation">(Castells and Castillo 2000)</span> provided some further understanding of the matter.</p> +<p>At the moment, the package consists of four functions (<a href="https://pkgdown.jrwb.de/chemCal/reference/calplot.lm.html">calplot</a>, <a href="https://pkgdown.jrwb.de/chemCal/reference/lod.html">lod</a>, <a href="https://pkgdown.jrwb.de/chemCal/reference/loq.html">loq</a> and <a href="https://pkgdown.jrwb.de/chemCal/reference/inverse.predict.html">inverse.predict</a>), working on univariate linear models of class <code>lm</code> or <code>rlm</code>, plus several datasets for validation.</p> +<p>A <a href="http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877">bug report</a> and the following e-mail exchange on the r-devel mailing list about prediction intervals from weighted regression entailed some further studies on this subject. However, I did not encounter any proof or explanation of the formula cited below yet, so I can’t really confirm that Massart’s method is correct.</p> +<p>In fact, in June 2018 I was made aware of the fact that the inverse prediction method implemented in chemCal version 0.1.37 and before did not take the variance of replicate calibration standards about their means into account, nor the number of replicates when calculating the degrees of freedom. Thanks to PhD student Anna Burniol Figols for reporting this issue!</p> +<p>As a consequence, I rewrote <code>inverse.predict</code> not to automatically work with the mean responses for each calibration standard any more. The example calculations from <span class="citation">Massart et al. (1997)</span> can still be reproduced when the regression model is calculated using the means of the calibration data as shown below.</p> +</div> +<div id="usage" class="section level1"> +<h1 class="hasAnchor"> +<a href="#usage" class="anchor"></a>Usage</h1> +<p>When calibrating an analytical method, the first task is to generate a suitable model. If we want to use the <code>chemCal</code> functions, we have to restrict ourselves to univariate, possibly weighted, linear regression so far.</p> +<p>Once such a model has been created, the calibration can be graphically shown by using the <code>calplot</code> function:</p> +<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(chemCal) +m0 <-<span class="st"> </span><span class="kw">lm</span>(y <span class="op">~</span><span class="st"> </span>x, <span class="dt">data =</span> massart97ex3) +<span class="kw"><a href="../reference/calplot.lm.html">calplot</a></span>(m0)</code></pre></div> +<p><img src="chemCal_files/figure-html/unnamed-chunk-1-1.png" width="700"></p> +<p>As we can see, the scatter increases with increasing x. This is also illustrated by one of the diagnostic plots for linear models provided by R:</p> +<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">plot</span>(m0, <span class="dt">which=</span><span class="dv">3</span>)</code></pre></div> +<p><img src="chemCal_files/figure-html/unnamed-chunk-2-1.png" width="700"></p> +<p>Therefore, in Example 8 in <span class="citation">Massart et al. (1997)</span>, weighted regression is proposed which can be reproduced by the following code. Note that we are building the model on the mean values for each standard in order to be able to reproduce the results given in the book with the current version of chemCal.</p> +<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">weights <-<span class="st"> </span><span class="kw">with</span>(massart97ex3, { + yx <-<span class="st"> </span><span class="kw">split</span>(y, x) + ybar <-<span class="st"> </span><span class="kw">sapply</span>(yx, mean) + s <-<span class="st"> </span><span class="kw">round</span>(<span class="kw">sapply</span>(yx, sd), <span class="dt">digits =</span> <span class="dv">2</span>) + w <-<span class="st"> </span><span class="kw">round</span>(<span class="dv">1</span> <span class="op">/</span><span class="st"> </span>(s<span class="op">^</span><span class="dv">2</span>), <span class="dt">digits =</span> <span class="dv">3</span>) +}) +massart97ex3.means <-<span class="st"> </span><span class="kw">aggregate</span>(y <span class="op">~</span><span class="st"> </span>x, massart97ex3, mean) + +m <-<span class="st"> </span><span class="kw">lm</span>(y <span class="op">~</span><span class="st"> </span>x, <span class="dt">w =</span> weights, <span class="dt">data =</span> massart97ex3.means)</code></pre></div> +<p>If we now want to predict a new x value from measured y values, we use the <code>inverse.predict</code> function:</p> +<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="../reference/inverse.predict.html">inverse.predict</a></span>(m, <span class="dv">15</span>, <span class="dt">ws=</span><span class="fl">1.67</span>)</code></pre></div> +<pre><code>## $Prediction +## [1] 5.865367 +## +## $`Standard Error` +## [1] 0.8926109 +## +## $Confidence +## [1] 2.478285 +## +## $`Confidence Limits` +## [1] 3.387082 8.343652</code></pre> +<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="../reference/inverse.predict.html">inverse.predict</a></span>(m, <span class="dv">90</span>, <span class="dt">ws =</span> <span class="fl">0.145</span>)</code></pre></div> +<pre><code>## $Prediction +## [1] 44.06025 +## +## $`Standard Error` +## [1] 2.829162 +## +## $Confidence +## [1] 7.855012 +## +## $`Confidence Limits` +## [1] 36.20523 51.91526</code></pre> +<p>The weight <code>ws</code> assigned to the measured y value has to be given by the user in the case of weighted regression, or alternatively, the approximate variance <code>var.s</code> at this location.</p> +</div> +<div id="background-for-inverse-predict" class="section level1"> +<h1 class="hasAnchor"> +<a href="#background-for-inverse-predict" class="anchor"></a>Background for <code>inverse.predict</code> +</h1> +<p>Equation 8.28 in <span class="citation">Massart et al. (1997)</span> gives a general equation for predicting the standard error <span class="math inline">\(s_{\hat{x_s}}\)</span> for an <span class="math inline">\(x\)</span> value predicted from measurements of <span class="math inline">\(y\)</span> according to the linear calibration function <span class="math inline">\(y = b_0 + b_1 \cdot x\)</span>:</p> +<span class="math display">\[\begin{equation} +s_{\hat{x_s}} = \frac{s_e}{b_1} \sqrt{\frac{1}{w_s m} + \frac{1}{\sum{w_i}} + + \frac{(\bar{y_s} - \bar{y_w})^2 \sum{w_i}} + {{b_1}^2 \left( \sum{w_i} \sum{w_i {x_i}^2} - + {\left( \sum{ w_i x_i } \right)}^2 \right) }} +\end{equation}\]</span> +<p>with</p> +<span class="math display">\[\begin{equation} +s_e = \sqrt{ \frac{\sum w_i (y_i - \hat{y_i})^2}{n - 2}} +\end{equation}\]</span> +<p>In chemCal version before 0.2, I interpreted <span class="math inline">\(w_i\)</span> to be the weight for calibration standard <span class="math inline">\(i\)</span>, <span class="math inline">\(y_i\)</span> to be the mean value observed for standard <span class="math inline">\(i\)</span>, and <span class="math inline">\(n\)</span> to be the number of calibration standards. With this implementation I was able to reproduce the examples given in the book. However, as noted above, I was made aware of the fact that this way of calculation does not take the variation of the y values about the means into account. Furthermore, I noticed that for the case of unweighted linear calibration with replicate standards, <code>inverse.predict</code> produced different results than <code>calibrate</code> from the <code>investr</code> package when using the Wald method.</p> +<p>Both issues are now addressed in chemCal starting from version 0.2.1. Here, <span class="math inline">\(y_i\)</span> is calibration measurement <span class="math inline">\(i\)</span>, <span class="math inline">\(\hat{y_i}\)</span> is the estimated value for calibration measurement <span class="math inline">\(i\)</span> and <span class="math inline">\(n\)</span> is the total number of calibration measurements.</p> +<p><span class="math inline">\(w_s\)</span> is the weight attributed to the sample <span class="math inline">\(s\)</span>, <span class="math inline">\(m\)</span> is the number of replicate measurements of sample <span class="math inline">\(s\)</span>, <span class="math inline">\(\bar{y_s}\)</span> is the mean response for the sample, <span class="math inline">\(\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}\)</span> is the weighted mean of responses <span class="math inline">\(y_i\)</span>, and <span class="math inline">\(x_i\)</span> is the given <span class="math inline">\(x\)</span> value for standard <span class="math inline">\(i\)</span>.</p> +<p>The weight <span class="math inline">\(w_s\)</span> for the sample should be estimated or calculated in accordance to the weights used in the linear regression.</p> +<p>I had also adjusted the above equation in order to be able to take a different precisions in standards and samples into account. In analogy to Equation 8.26 from I am using</p> +<span class="math display">\[\begin{equation} +s_{\hat{x_s}} = \frac{1}{b_1} \sqrt{\frac{{s_s}^2}{w_s m} + + {s_e}^2 \left( \frac{1}{\sum{w_i}} + + \frac{(\bar{y_s} - \bar{y_w})^2 \sum{w_i}} + {{b_1}^2 \left( \sum{w_i} \sum{w_i {x_i}^2} - {\left( \sum{ w_i x_i } \right)}^2 \right) } \right) } +\end{equation}\]</span> +<p>where I interpret <span class="math inline">\(\frac{{s_s}^2}{w_s}\)</span> as an estimator of the variance at location <span class="math inline">\(\hat{x_s}\)</span>, which can be replaced by a user-specified value using the argument <code>var.s</code> of the function <code>inverse.predict</code>.</p> +<div id="refs" class="references"> +<div id="ref-amc89"> +<p>Analytical Methods Committee. 1989. “Robust Statistics — How Not to Reject Outliers. Part 1. Basic Concepts.” <em>The Analyst</em> 114: 1693–7.</p> +</div> +<div id="ref-castells00"> +<p>Castells, Reynaldo César, and Marcela Alejandra Castillo. 2000. “Systematic Errors: Detection and Correction by Means of Standard Calibration, Youden Calibration and Standard Additions Method in Conjunction with a Method Response Model.” <em>Analytica Chimica Acta</em> 423: 179–85.</p> +</div> +<div id="ref-currie97"> +<p>Currie, L. A. 1997. “Nomenclature in Evaluation of Analytical Methods Including Detection and Quantification Capabilities (IUPAC Recommendations 1995).” <em>Analytica Chimica Acta</em> 391: 105–26.</p> +</div> +<div id="ref-massart97"> +<p>Massart, D. L, B. G. M. Vandeginste, L. M. C. Buydens, S. De Jong, P. J. Lewi, and J Smeyers-Verbeke. 1997. <em>Handbook of Chemometrics and Qualimetrics: Part A</em>. Amsterdam: Elsevier.</p> +</div> +</div> +</div> + </div> + + <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> + <div id="tocnav"> + <h2 class="hasAnchor"> +<a href="#tocnav" class="anchor"></a>Contents</h2> + <ul class="nav nav-pills nav-stacked"> +<li><a href="#basic-calibration-functions">Basic calibration functions</a></li> + <li><a href="#usage">Usage</a></li> + <li><a href="#background-for-inverse-predict">Background for <code>inverse.predict</code></a></li> + </ul> +</div> + </div> + +</div> + + + <footer><div class="copyright"> + <p>Developed by Johannes Ranke.</p> +</div> + +<div class="pkgdown"> + <p>Site built with <a href="http://pkgdown.r-lib.org/">pkgdown</a>.</p> +</div> + + </footer> +</div> + + + + </body> +</html> diff --git a/docs/articles/chemCal_files/figure-html/unnamed-chunk-1-1.png b/docs/articles/chemCal_files/figure-html/unnamed-chunk-1-1.png Binary files differnew file mode 100644 index 0000000..62c95ca --- /dev/null +++ b/docs/articles/chemCal_files/figure-html/unnamed-chunk-1-1.png diff --git a/docs/articles/chemCal_files/figure-html/unnamed-chunk-2-1.png b/docs/articles/chemCal_files/figure-html/unnamed-chunk-2-1.png Binary files differnew file mode 100644 index 0000000..464048c --- /dev/null +++ b/docs/articles/chemCal_files/figure-html/unnamed-chunk-2-1.png diff --git a/docs/articles/chemCal_files/figure-html/unnamed-chunk-3-1.png b/docs/articles/chemCal_files/figure-html/unnamed-chunk-3-1.png Binary files differnew file mode 100644 index 0000000..464048c --- /dev/null +++ b/docs/articles/chemCal_files/figure-html/unnamed-chunk-3-1.png diff --git a/docs/articles/index.html b/docs/articles/index.html index 2e993ac..a29b7c6 100644 --- a/docs/articles/index.html +++ b/docs/articles/index.html @@ -18,23 +18,35 @@ <!-- 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"> +<!-- clipboard.js --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/1.7.1/clipboard.min.js" integrity="sha384-cV+rhyOuRHc9Ub/91rihWcGmMmCXDeksTtCihMupQHSsi8GIIRDG0ThDc3HGQFJ3" crossorigin="anonymous"></script> + +<!-- sticky kit --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/sticky-kit/1.1.3/sticky-kit.min.js" integrity="sha256-c4Rlo1ZozqTPE2RLuvbusY3+SU1pQaJC0TjuhygMipw=" crossorigin="anonymous"></script> <!-- pkgdown --> <link href="../pkgdown.css" rel="stylesheet"> -<script src="../jquery.sticky-kit.min.js"></script> <script src="../pkgdown.js"></script> + + +<meta property="og:title" content="Articles" /> + + + <!-- mathjax --> -<script src='https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML'></script> +<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-vignette-index"> + <div class="container template-article-index"> <header> <div class="navbar navbar-default navbar-fixed-top" role="navigation"> <div class="container"> @@ -44,21 +56,35 @@ <span class="icon-bar"></span> <span class="icon-bar"></span> </button> - <a class="navbar-brand" href="../index.html">chemCal</a> + <span class="navbar-brand"> + <a class="navbar-link" href="../index.html">chemCal</a> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> + </span> </div> + <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> +</li> +<li> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../news/index.html">Changelog</a> </li> </ul> <ul class="nav navbar-nav navbar-right"> </ul> + </div><!--/.nav-collapse --> </div><!--/.container --> </div><!--/.navbar --> @@ -66,17 +92,18 @@ </header> - <div class="page-header"> - <h1>Vignette reference <small>version 0.1-37.9001</small></h1> -</div> - <div class="row"> - <div class="col-md-9"> + <div class="col-md-9 contents"> + <div class="page-header"> + <h1>Articles</h1> + </div> + <div class="section "> <h3>All vignettes</h3> <p class="section-desc"></p> <ul> + <li><a href="chemCal.html">Introduction to chemCal</a></li> </ul> </div> </div> @@ -88,11 +115,14 @@ </div> <div class="pkgdown"> - <p>Site built with <a href="http://hadley.github.io/pkgdown/">pkgdown</a>.</p> + <p>Site built with <a href="http://pkgdown.r-lib.org/">pkgdown</a>.</p> </div> </footer> </div> + + </body> </html> + diff --git a/docs/authors.html b/docs/authors.html index c9890cd..ea55846 100644 --- a/docs/authors.html +++ b/docs/authors.html @@ -58,17 +58,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="reference/index.html">Functions and data</a> + <a href="index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="articles/chemCal.pdf">Short manual (pdf)</a> + <a href="articles/chemCal.html">Get started</a> +</li> +<li> + <a href="reference/index.html">Reference</a> +</li> +<li> + <a href="news/index.html">Changelog</a> </li> </ul> diff --git a/docs/index.html b/docs/index.html index 6205aa1..b2233bb 100644 --- a/docs/index.html +++ b/docs/index.html @@ -35,17 +35,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="reference/index.html">Functions and data</a> + <a href="index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="articles/chemCal.pdf">Short manual (pdf)</a> + <a href="articles/chemCal.html">Get started</a> +</li> +<li> + <a href="reference/index.html">Reference</a> +</li> +<li> + <a href="news/index.html">Changelog</a> </li> </ul> <ul class="nav navbar-nav navbar-right"></ul> @@ -73,7 +82,7 @@ <ul class="list-unstyled"> <li>Download from CRAN at <br><a href="https://cloud.r-project.org/package=chemCal">https://cloud.r-project.org/package=chemCal</a> </li> -<li>Report a bug at <br><a href="NA">NA</a> +<li>Report a bug at <br><a href="http://github.com/jranke/chemCal/issues">http://github.com/jranke/chemCal/issues</a> </li> </ul> </div> diff --git a/docs/news/index.html b/docs/news/index.html new file mode 100644 index 0000000..62ab315 --- /dev/null +++ b/docs/news/index.html @@ -0,0 +1,153 @@ +<!-- 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>Changelog • chemCal</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"> + +<!-- clipboard.js --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/1.7.1/clipboard.min.js" integrity="sha384-cV+rhyOuRHc9Ub/91rihWcGmMmCXDeksTtCihMupQHSsi8GIIRDG0ThDc3HGQFJ3" crossorigin="anonymous"></script> + +<!-- sticky kit --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/sticky-kit/1.1.3/sticky-kit.min.js" integrity="sha256-c4Rlo1ZozqTPE2RLuvbusY3+SU1pQaJC0TjuhygMipw=" crossorigin="anonymous"></script> + +<!-- pkgdown --> +<link href="../pkgdown.css" rel="stylesheet"> +<script src="../pkgdown.js"></script> + + + +<meta property="og:title" content="Changelog" /> + + + +<!-- 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-news"> + <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> + <span class="navbar-brand"> + <a class="navbar-link" href="../index.html">chemCal</a> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> + </span> + </div> + + <div id="navbar" class="navbar-collapse collapse"> + <ul class="nav navbar-nav"> + <li> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> +</li> +<li> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</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>Changelog <small></small></h1> + + </div> + + <div id="chemcal-0-2-1-2018-07-17" class="section level1"> +<h1 class="page-header"> +<a href="#chemcal-0-2-1-2018-07-17" class="anchor"></a>chemCal 0.2.1 (2018-07-17)<small> Unreleased </small> +</h1> +<ul> +<li><p>‘inverse.predict’: Do not work on the means of the calibration standards any more, as this ignores the variability of y values about the means</p></li> +<li><p>Use testthat for tests to simplify further development. Adapt the tests using data with replicate standard measurements to work on the means in order to show the relation to ‘inverse.predict’ from earlier versions. Include comparisons with investr::calibrate(method = ‘Wald’) for unweighted regressions. Include tests with more precision to check for changes in numerical output across versions.</p></li> +<li><p>‘lod’ and ‘loq’: In the lists that are returned, return the list component ‘y’ without names, because we always only have a single element in ‘y’ (previously the name ‘1’ was returned).</p></li> +<li><p>Convert vignette to html and explain the changes to ‘inverse.predict’</p></li> +<li><p>Add two example dataset, one from an online course at the University of Toronto, one from Rocke and Lorenzato (1995)</p></li> +</ul> +</div> + <div id="chemcal-0-1-33-2014-04-24" class="section level1"> +<h1 class="page-header"> +<a href="#chemcal-0-1-33-2014-04-24" class="anchor"></a>chemCal 0.1-33 (2014-04-24)<small> Unreleased </small> +</h1> +<ul> +<li>Bugfix in lod() and loq(): In the case of small absolute x values (e.g. on the order of 1e-4 and smaller), the lod or loq calculated using the default method could produce inaccurate results as the default tolerance that was used in the internal call to optimize is inappropriate in such cases. Now a reasonable default is used which can be overriden by the user. Thanks to Jérôme Ambroise for reporting the bug.</li> +</ul> +<p>For a detailed list of changes to the chemCal source please consult the commit history on <a href="https://cgit.jrwb.de/chemCal" class="uri">https://cgit.jrwb.de/chemCal</a></p> +</div> + </div> + + <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> + <div id="tocnav"> + <h2>Contents</h2> + <ul class="nav nav-pills nav-stacked"> + <li><a href="#chemcal-0-2-1-2018-07-17">0.2.1</a></li> + <li><a href="#chemcal-0-1-33-2014-04-24">0.1-33</a></li> + </ul> + </div> + </div> + +</div> + + <footer> + <div class="copyright"> + <p>Developed by Johannes Ranke.</p> +</div> + +<div class="pkgdown"> + <p>Site built with <a href="http://pkgdown.r-lib.org/">pkgdown</a>.</p> +</div> + + </footer> + </div> + + + + </body> +</html> + diff --git a/docs/pkgdown.yml b/docs/pkgdown.yml index 6e67788..3b5e371 100644 --- a/docs/pkgdown.yml +++ b/docs/pkgdown.yml @@ -1,5 +1,6 @@ pandoc: 1.19.2.4 pkgdown: 1.1.0 pkgdown_sha: ~ -articles: [] +articles: + chemCal: chemCal.html diff --git a/docs/reference/calplot.lm.html b/docs/reference/calplot.lm.html index 7eda066..a5b3bbb 100644 --- a/docs/reference/calplot.lm.html +++ b/docs/reference/calplot.lm.html @@ -62,17 +62,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> @@ -134,7 +143,7 @@ <th>alpha</th> <td><p>The error tolerance level for the confidence and prediction bands. Note that this includes both tails of the Gaussian distribution, unlike the alpha and beta parameters - used in <code>lod</code> (see note below).</p></td> + used in <code><a href='lod.html'>lod</a></code> (see note below).</p></td> </tr> <tr> <th>varfunc</th> @@ -157,7 +166,7 @@ the internally used function <code>predict.lm</code>, therefore, <code>calplot</code> does not draw prediction bands for such models.</p> <p>It is possible to compare the <code>calplot</code> prediction bands with the - <code>lod</code> values if the <code>lod()</code> alpha and beta parameters are + <code><a href='lod.html'>lod</a></code> values if the <code><a href='lod.html'>lod()</a></code> alpha and beta parameters are half the value of the <code>calplot()</code> alpha parameter.</p> diff --git a/docs/reference/chemCal-package.html b/docs/reference/chemCal-package.html index 8e6f39b..cadd357 100644 --- a/docs/reference/chemCal-package.html +++ b/docs/reference/chemCal-package.html @@ -61,17 +61,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> diff --git a/docs/reference/din32645.html b/docs/reference/din32645.html index 03170e0..dbc9504 100644 --- a/docs/reference/din32645.html +++ b/docs/reference/din32645.html @@ -61,17 +61,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> @@ -118,33 +127,17 @@ <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> - <pre class="examples"><div class='input'><span class='fu'>data</span>(<span class='no'>din32645</span>) -<span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>din32645</span>) -<span class='fu'>calplot</span>(<span class='no'>m</span>)</div><div class='img'><img src='din32645-1.png' alt='' width='700' height='433' /></div><div class='input'> + <pre class="examples"><div class='input'><span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>din32645</span>) +<span class='fu'><a href='calplot.lm.html'>calplot</a></span>(<span class='no'>m</span>)</div><div class='img'><img src='din32645-1.png' alt='' width='700' height='433' /></div><div class='input'> <span class='co'>## Prediction of x with confidence interval</span> -(<span class='no'>prediction</span> <span class='kw'><-</span> <span class='fu'>inverse.predict</span>(<span class='no'>m</span>, <span class='fl'>3500</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>))</div><div class='output co'>#> $Prediction -#> [1] 0.1054792 -#> -#> $`Standard Error` -#> [1] 0.02215619 -#> -#> $Confidence -#> [1] 0.07434261 -#> -#> $`Confidence Limits` -#> [1] 0.03113656 0.17982178 -#> </div><div class='input'> +<span class='no'>prediction</span> <span class='kw'><-</span> <span class='fu'><a href='inverse.predict.html'>inverse.predict</a></span>(<span class='no'>m</span>, <span class='fl'>3500</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) + <span class='co'># This should give 0.07434 according to test data from Dintest, which </span> <span class='co'># was collected from Procontrol 3.1 (isomehr GmbH) in this case</span> -<span class='fu'>round</span>(<span class='no'>prediction</span>$<span class='no'>Confidence</span>,<span class='fl'>5</span>)</div><div class='output co'>#> [1] 0.07434</div><div class='input'> +<span class='fu'>round</span>(<span class='no'>prediction</span>$<span class='no'>Confidence</span>, <span class='fl'>5</span>)</div><div class='output co'>#> [1] 0.07434</div><div class='input'> <span class='co'>## Critical value:</span> -(<span class='no'>crit</span> <span class='kw'><-</span> <span class='fu'>lod</span>(<span class='no'>m</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>beta</span> <span class='kw'>=</span> <span class='fl'>0.5</span>))</div><div class='output co'>#> $x -#> [1] 0.0698127 -#> -#> $y -#> 1 -#> 3155.393 -#> </div><div class='input'> +<span class='no'>crit</span> <span class='kw'><-</span> <span class='fu'><a href='lod.html'>lod</a></span>(<span class='no'>m</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>beta</span> <span class='kw'>=</span> <span class='fl'>0.5</span>) + <span class='co'># According to DIN 32645, we should get 0.07 for the critical value</span> <span class='co'># (decision limit, "Nachweisgrenze")</span> <span class='fu'>round</span>(<span class='no'>crit</span>$<span class='no'>x</span>, <span class='fl'>2</span>)</div><div class='output co'>#> [1] 0.07</div><div class='input'><span class='co'># and according to Dintest test data, we should get 0.0698 from</span> @@ -153,21 +146,16 @@ <span class='co'># In German, the smallest detectable value is the "Erfassungsgrenze", and we</span> <span class='co'># should get 0.14 according to DIN, which we achieve by using the method </span> <span class='co'># described in it:</span> -<span class='no'>lod.din</span> <span class='kw'><-</span> <span class='fu'>lod</span>(<span class='no'>m</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>beta</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>method</span> <span class='kw'>=</span> <span class='st'>"din"</span>) +<span class='no'>lod.din</span> <span class='kw'><-</span> <span class='fu'><a href='lod.html'>lod</a></span>(<span class='no'>m</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>beta</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>method</span> <span class='kw'>=</span> <span class='st'>"din"</span>) <span class='fu'>round</span>(<span class='no'>lod.din</span>$<span class='no'>x</span>, <span class='fl'>2</span>)</div><div class='output co'>#> [1] 0.14</div><div class='input'> <span class='co'>## Limit of quantification</span> <span class='co'># This accords to the test data coming with the test data from Dintest again, </span> <span class='co'># except for the last digits of the value cited for Procontrol 3.1 (0.2121)</span> -(<span class='no'>loq</span> <span class='kw'><-</span> <span class='fu'>loq</span>(<span class='no'>m</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>))</div><div class='output co'>#> $x -#> [1] 0.2119575 -#> -#> $y -#> 1 -#> 4528.787 -#> </div><div class='input'><span class='fu'>round</span>(<span class='no'>loq</span>$<span class='no'>x</span>,<span class='fl'>4</span>)</div><div class='output co'>#> [1] 0.212</div><div class='input'> +<span class='no'>loq</span> <span class='kw'><-</span> <span class='fu'><a href='loq.html'>loq</a></span>(<span class='no'>m</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>) +<span class='fu'>round</span>(<span class='no'>loq</span>$<span class='no'>x</span>, <span class='fl'>4</span>)</div><div class='output co'>#> [1] 0.212</div><div class='input'> <span class='co'># A similar value is obtained using the approximation </span> <span class='co'># LQ = 3.04 * LC (Currie 1999, p. 120)</span> -<span class='fl'>3.04</span> * <span class='fu'>lod</span>(<span class='no'>m</span>,<span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>beta</span> <span class='kw'>=</span> <span class='fl'>0.5</span>)$<span class='no'>x</span></div><div class='output co'>#> [1] 0.2122306</div></pre> +<span class='fl'>3.04</span> * <span class='fu'><a href='lod.html'>lod</a></span>(<span class='no'>m</span>, <span class='kw'>alpha</span> <span class='kw'>=</span> <span class='fl'>0.01</span>, <span class='kw'>beta</span> <span class='kw'>=</span> <span class='fl'>0.5</span>)$<span class='no'>x</span></div><div class='output co'>#> [1] 0.2122306</div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> <h2>Contents</h2> diff --git a/docs/reference/index.html b/docs/reference/index.html index ce4f346..d314ec2 100644 --- a/docs/reference/index.html +++ b/docs/reference/index.html @@ -58,17 +58,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> @@ -152,6 +161,18 @@ <p><code><a href="massart97ex3.html">massart97ex3</a></code> </p> </td> <td><p>Calibration data from Massart et al. (1997), example 3</p></td> + </tr><tr> + + <td> + <p><code><a href="rl95_toluene.html">rl95_toluene</a></code> </p> + </td> + <td><p>Toluene amounts by GC/MS as reported by Rocke and Lorenzato (1995)</p></td> + </tr><tr> + + <td> + <p><code><a href="utstats14.html">utstats14</a></code> </p> + </td> + <td><p>Example data for calibration with replicates from University of Toronto</p></td> </tr> </tbody> </table> diff --git a/docs/reference/inverse.predict.html b/docs/reference/inverse.predict.html index 672e322..b107154 100644 --- a/docs/reference/inverse.predict.html +++ b/docs/reference/inverse.predict.html @@ -70,17 +70,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> @@ -164,7 +173,11 @@ <h2 class="hasAnchor" id="note"><a class="anchor" href="#note"></a>Note</h2> - <p>The function was validated with examples 7 and 8 from Massart et al. (1997).</p> + <p>The function was validated with examples 7 and 8 from Massart et al. (1997). + Note that the behaviour of inverse.predict changed with chemCal version + 0.2.1. Confidence intervals for x values obtained from calibrations with + replicate measurements did not take the variation about the means into account. + Please refer to the vignette for details.</p> <h2 class="hasAnchor" id="references"><a class="anchor" href="#references"></a>References</h2> @@ -175,7 +188,6 @@ <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> <pre class="examples"><div class='input'><span class='co'># This is example 7 from Chapter 8 in Massart et al. (1997)</span> -<span class='fu'>data</span>(<span class='no'>massart97ex1</span>) <span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>massart97ex1</span>) <span class='fu'>inverse.predict</span>(<span class='no'>m</span>, <span class='fl'>15</span>) <span class='co'># 6.1 +- 4.9</span></div><div class='output co'>#> $Prediction #> [1] 6.09381 @@ -210,7 +222,44 @@ #> #> $`Confidence Limits` #> [1] 40.77134 47.10832 -#> </div></pre> +#> </div><div class='input'> +<span class='co'># For reproducing the results for replicate standard measurements in example 8,</span> +<span class='co'># we need to do the calibration on the means when using chemCal > 0.2</span> +<span class='no'>weights</span> <span class='kw'><-</span> <span class='fu'>with</span>(<span class='no'>massart97ex3</span>, { + <span class='no'>yx</span> <span class='kw'><-</span> <span class='fu'>split</span>(<span class='no'>y</span>, <span class='no'>x</span>) + <span class='no'>ybar</span> <span class='kw'><-</span> <span class='fu'>sapply</span>(<span class='no'>yx</span>, <span class='no'>mean</span>) + <span class='no'>s</span> <span class='kw'><-</span> <span class='fu'>round</span>(<span class='fu'>sapply</span>(<span class='no'>yx</span>, <span class='no'>sd</span>), <span class='kw'>digits</span> <span class='kw'>=</span> <span class='fl'>2</span>) + <span class='no'>w</span> <span class='kw'><-</span> <span class='fu'>round</span>(<span class='fl'>1</span> / (<span class='no'>s</span>^<span class='fl'>2</span>), <span class='kw'>digits</span> <span class='kw'>=</span> <span class='fl'>3</span>) +}) + +<span class='no'>massart97ex3.means</span> <span class='kw'><-</span> <span class='fu'>aggregate</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='no'>massart97ex3</span>, <span class='no'>mean</span>) + +<span class='no'>m3.means</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>w</span> <span class='kw'>=</span> <span class='no'>weights</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>massart97ex3.means</span>) + +<span class='fu'>inverse.predict</span>(<span class='no'>m3.means</span>, <span class='fl'>15</span>, <span class='kw'>ws</span> <span class='kw'>=</span> <span class='fl'>1.67</span>) <span class='co'># 5.9 +- 2.5</span></div><div class='output co'>#> $Prediction +#> [1] 5.865367 +#> +#> $`Standard Error` +#> [1] 0.8926109 +#> +#> $Confidence +#> [1] 2.478285 +#> +#> $`Confidence Limits` +#> [1] 3.387082 8.343652 +#> </div><div class='input'><span class='fu'>inverse.predict</span>(<span class='no'>m3.means</span>, <span class='fl'>90</span>, <span class='kw'>ws</span> <span class='kw'>=</span> <span class='fl'>0.145</span>) <span class='co'># 44.1 +- 7.9</span></div><div class='output co'>#> $Prediction +#> [1] 44.06025 +#> +#> $`Standard Error` +#> [1] 2.829162 +#> +#> $Confidence +#> [1] 7.855012 +#> +#> $`Confidence Limits` +#> [1] 36.20523 51.91526 +#> </div><div class='input'> +</div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> <h2>Contents</h2> diff --git a/docs/reference/lod.html b/docs/reference/lod.html index 7519f19..b3d9574 100644 --- a/docs/reference/lod.html +++ b/docs/reference/lod.html @@ -68,17 +68,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> @@ -184,18 +193,16 @@ <h2 class="hasAnchor" id="see-also"><a class="anchor" href="#see-also"></a>See also</h2> - <div class='dont-index'><p>Examples for <code>din32645</code></p></div> + <div class='dont-index'><p>Examples for <code><a href='din32645.html'>din32645</a></code></p></div> <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> - <pre class="examples"><div class='input'><span class='fu'>data</span>(<span class='no'>din32645</span>) -<span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>din32645</span>) + <pre class="examples"><div class='input'><span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>din32645</span>) <span class='fu'>lod</span>(<span class='no'>m</span>)</div><div class='output co'>#> $x #> [1] 0.08655484 #> #> $y -#> 1 -#> 3317.154 +#> [1] 3317.154 #> </div><div class='input'> <span class='co'># The critical value (decision limit, German Nachweisgrenze) can be obtained</span> <span class='co'># by using beta = 0.5:</span> @@ -203,8 +210,7 @@ #> [1] 0.0698127 #> #> $y -#> 1 -#> 3155.393 +#> [1] 3155.393 #> </div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> diff --git a/docs/reference/loq.html b/docs/reference/loq.html index 09567db..eee0925 100644 --- a/docs/reference/loq.html +++ b/docs/reference/loq.html @@ -67,17 +67,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> @@ -108,7 +117,7 @@ $$L = k c(L)$$ where c(L) is half of the length of the confidence interval at the limit L (DIN 32645, equivalent to ISO 11843). c(L) is internally estimated by - <code>inverse.predict</code>, and L is obtained by iteration.</p> + <code><a href='inverse.predict.html'>inverse.predict</a></code>, and L is obtained by iteration.</p> </div> @@ -150,7 +159,7 @@ <th>w.loq</th> <td><p>The weight that should be attributed to the LOQ. Defaults to one for unweighted regression, and to the mean of the weights - for weighted regression. See <code>massart97ex3</code> for + for weighted regression. See <code><a href='massart97ex3.html'>massart97ex3</a></code> for an example how to take advantage of knowledge about the variance function.</p></td> </tr> @@ -180,27 +189,23 @@ <h2 class="hasAnchor" id="see-also"><a class="anchor" href="#see-also"></a>See also</h2> - <div class='dont-index'><p>Examples for <code>din32645</code></p></div> + <div class='dont-index'><p>Examples for <code><a href='din32645.html'>din32645</a></code></p></div> <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> - <pre class="examples"><div class='input'><span class='fu'>data</span>(<span class='no'>massart97ex3</span>) -<span class='fu'>attach</span>(<span class='no'>massart97ex3</span>) -<span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>) + <pre class="examples"><div class='input'><span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>massart97ex1</span>) <span class='fu'>loq</span>(<span class='no'>m</span>)</div><div class='output co'>#> $x #> [1] 13.97764 #> #> $y -#> 1 -#> 30.6235 +#> [1] 30.6235 #> </div><div class='input'> <span class='co'># We can get better by using replicate measurements</span> <span class='fu'>loq</span>(<span class='no'>m</span>, <span class='kw'>n</span> <span class='kw'>=</span> <span class='fl'>3</span>)</div><div class='output co'>#> $x #> [1] 9.971963 #> #> $y -#> 1 -#> 22.68539 +#> [1] 22.68539 #> </div></pre> </div> <div class="col-md-3 hidden-xs hidden-sm" id="sidebar"> diff --git a/docs/reference/massart97ex1.html b/docs/reference/massart97ex1.html index ddd1ce1..844023f 100644 --- a/docs/reference/massart97ex1.html +++ b/docs/reference/massart97ex1.html @@ -61,17 +61,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> diff --git a/docs/reference/massart97ex3.html b/docs/reference/massart97ex3.html index 434fc7b..e9fa06a 100644 --- a/docs/reference/massart97ex3.html +++ b/docs/reference/massart97ex3.html @@ -61,17 +61,26 @@ </button> <span class="navbar-brand"> <a class="navbar-link" href="../index.html">chemCal</a> - <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.1.37.9001</span> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> </span> </div> <div id="navbar" class="navbar-collapse collapse"> <ul class="nav navbar-nav"> <li> - <a href="../reference/index.html">Functions and data</a> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> </li> <li> - <a href="../articles/chemCal.pdf">Short manual (pdf)</a> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</a> </li> </ul> @@ -100,7 +109,7 @@ </div> - <pre class="usage"><span class='fu'>data</span>(<span class='no'>massart97ex3</span>)</pre> + <pre class="usage"><span class='no'>massart97ex3</span></pre> <h2 class="hasAnchor" id="format"><a class="anchor" href="#format"></a>Format</h2> @@ -115,18 +124,21 @@ <h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2> - <pre class="examples"><div class='input'><span class='fu'>data</span>(<span class='no'>massart97ex3</span>) -<span class='fu'>attach</span>(<span class='no'>massart97ex3</span>)</div><div class='output co'>#> <span class='message'>The following objects are masked from massart97ex3 (pos = 3):</span> -#> <span class='message'></span> -#> <span class='message'> x, y</span></div><div class='input'><span class='no'>yx</span> <span class='kw'><-</span> <span class='fu'>split</span>(<span class='no'>y</span>, <span class='no'>x</span>) -<span class='no'>ybar</span> <span class='kw'><-</span> <span class='fu'>sapply</span>(<span class='no'>yx</span>, <span class='no'>mean</span>) -<span class='no'>s</span> <span class='kw'><-</span> <span class='fu'>round</span>(<span class='fu'>sapply</span>(<span class='no'>yx</span>, <span class='no'>sd</span>), <span class='kw'>digits</span> <span class='kw'>=</span> <span class='fl'>2</span>) -<span class='no'>w</span> <span class='kw'><-</span> <span class='fu'>round</span>(<span class='fl'>1</span> / (<span class='no'>s</span>^<span class='fl'>2</span>), <span class='kw'>digits</span> <span class='kw'>=</span> <span class='fl'>3</span>) -<span class='no'>weights</span> <span class='kw'><-</span> <span class='no'>w</span>[<span class='fu'>factor</span>(<span class='no'>x</span>)] -<span class='no'>m</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>w</span> <span class='kw'>=</span> <span class='no'>weights</span>) -<span class='fu'>calplot</span>(<span class='no'>m</span>)</div><div class='output co'>#> <span class='warning'>Warning: Assuming constant prediction variance even though model fit is weighted</span></div><div class='img'><img src='massart97ex3-1.png' alt='' width='700' height='433' /></div><div class='input'> + <pre class="examples"><div class='input'><span class='co'># For reproducing the results for replicate standard measurements in example 8,</span> +<span class='co'># we need to do the calibration on the means when using chemCal > 0.2</span> +<span class='no'>weights</span> <span class='kw'><-</span> <span class='fu'>with</span>(<span class='no'>massart97ex3</span>, { + <span class='no'>yx</span> <span class='kw'><-</span> <span class='fu'>split</span>(<span class='no'>y</span>, <span class='no'>x</span>) + <span class='no'>ybar</span> <span class='kw'><-</span> <span class='fu'>sapply</span>(<span class='no'>yx</span>, <span class='no'>mean</span>) + <span class='no'>s</span> <span class='kw'><-</span> <span class='fu'>round</span>(<span class='fu'>sapply</span>(<span class='no'>yx</span>, <span class='no'>sd</span>), <span class='kw'>digits</span> <span class='kw'>=</span> <span class='fl'>2</span>) + <span class='no'>w</span> <span class='kw'><-</span> <span class='fu'>round</span>(<span class='fl'>1</span> / (<span class='no'>s</span>^<span class='fl'>2</span>), <span class='kw'>digits</span> <span class='kw'>=</span> <span class='fl'>3</span>) +}) + +<span class='no'>massart97ex3.means</span> <span class='kw'><-</span> <span class='fu'>aggregate</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='no'>massart97ex3</span>, <span class='no'>mean</span>) + +<span class='no'>m3.means</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>w</span> <span class='kw'>=</span> <span class='no'>weights</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>massart97ex3.means</span>) + <span class='co'># The following concords with the book p. 200</span> -<span class='fu'>inverse.predict</span>(<span class='no'>m</span>, <span class='fl'>15</span>, <span class='kw'>ws</span> <span class='kw'>=</span> <span class='fl'>1.67</span>) <span class='co'># 5.9 +- 2.5</span></div><div class='output co'>#> $Prediction +<span class='fu'><a href='inverse.predict.html'>inverse.predict</a></span>(<span class='no'>m3.means</span>, <span class='fl'>15</span>, <span class='kw'>ws</span> <span class='kw'>=</span> <span class='fl'>1.67</span>) <span class='co'># 5.9 +- 2.5</span></div><div class='output co'>#> $Prediction #> [1] 5.865367 #> #> $`Standard Error` @@ -137,7 +149,7 @@ #> #> $`Confidence Limits` #> [1] 3.387082 8.343652 -#> </div><div class='input'><span class='fu'>inverse.predict</span>(<span class='no'>m</span>, <span class='fl'>90</span>, <span class='kw'>ws</span> <span class='kw'>=</span> <span class='fl'>0.145</span>) <span class='co'># 44.1 +- 7.9</span></div><div class='output co'>#> $Prediction +#> </div><div class='input'><span class='fu'><a href='inverse.predict.html'>inverse.predict</a></span>(<span class='no'>m3.means</span>, <span class='fl'>90</span>, <span class='kw'>ws</span> <span class='kw'>=</span> <span class='fl'>0.145</span>) <span class='co'># 44.1 +- 7.9</span></div><div class='output co'>#> $Prediction #> [1] 44.06025 #> #> $`Standard Error` @@ -151,33 +163,30 @@ #> </div><div class='input'> <span class='co'># The LOD is only calculated for models from unweighted regression</span> <span class='co'># with this version of chemCal</span> -<span class='no'>m0</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>) -<span class='fu'>lod</span>(<span class='no'>m0</span>)</div><div class='output co'>#> $x +<span class='no'>m0</span> <span class='kw'><-</span> <span class='fu'>lm</span>(<span class='no'>y</span> ~ <span class='no'>x</span>, <span class='kw'>data</span> <span class='kw'>=</span> <span class='no'>massart97ex3</span>) +<span class='fu'><a href='lod.html'>lod</a></span>(<span class='no'>m0</span>)</div><div class='output co'>#> $x #> [1] 5.407085 #> #> $y -#> 1 -#> 13.63911 +#> [1] 13.63911 #> </div><div class='input'> <span class='co'># Limit of quantification from unweighted regression</span> -<span class='fu'>loq</span>(<span class='no'>m0</span>)</div><div class='output co'>#> $x -#> [1] 13.97764 +<span class='fu'><a href='loq.html'>loq</a></span>(<span class='no'>m0</span>)</div><div class='output co'>#> $x +#> [1] 9.627349 #> #> $y -#> 1 -#> 30.6235 +#> [1] 22.00246 #> </div><div class='input'> <span class='co'># For calculating the limit of quantification from a model from weighted</span> <span class='co'># regression, we need to supply weights, internally used for inverse.predict</span> <span class='co'># If we are not using a variance function, we can use the weight from</span> <span class='co'># the above example as a first approximation (x = 15 is close to our</span> <span class='co'># loq approx 14 from above).</span> -<span class='fu'>loq</span>(<span class='no'>m</span>, <span class='kw'>w.loq</span> <span class='kw'>=</span> <span class='fl'>1.67</span>)</div><div class='output co'>#> $x +<span class='fu'><a href='loq.html'>loq</a></span>(<span class='no'>m3.means</span>, <span class='kw'>w.loq</span> <span class='kw'>=</span> <span class='fl'>1.67</span>)</div><div class='output co'>#> $x #> [1] 7.346195 #> #> $y -#> 1 -#> 17.90777 +#> [1] 17.90777 #> </div><div class='input'># The weight for the loq should therefore be derived at x = 7.3 instead # of 15, but the graphical procedure of Massart (p. 201) to derive the # variances on which the weights are based is quite inaccurate anyway. diff --git a/docs/reference/rl95_toluene.html b/docs/reference/rl95_toluene.html new file mode 100644 index 0000000..b252760 --- /dev/null +++ b/docs/reference/rl95_toluene.html @@ -0,0 +1,153 @@ +<!-- 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>Toluene amounts by GC/MS as reported by Rocke and Lorenzato (1995) — rl95_toluene • chemCal</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"> + +<!-- clipboard.js --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/1.7.1/clipboard.min.js" integrity="sha384-cV+rhyOuRHc9Ub/91rihWcGmMmCXDeksTtCihMupQHSsi8GIIRDG0ThDc3HGQFJ3" crossorigin="anonymous"></script> + +<!-- sticky kit --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/sticky-kit/1.1.3/sticky-kit.min.js" integrity="sha256-c4Rlo1ZozqTPE2RLuvbusY3+SU1pQaJC0TjuhygMipw=" crossorigin="anonymous"></script> + +<!-- pkgdown --> +<link href="../pkgdown.css" rel="stylesheet"> +<script src="../pkgdown.js"></script> + + + +<meta property="og:title" content="Toluene amounts by GC/MS as reported by Rocke and Lorenzato (1995) — rl95_toluene" /> + +<meta property="og:description" content="Dataset reproduced from Table 4 in Rocke and Lorenzato (1995)." /> +<meta name="twitter:card" content="summary" /> + + + +<!-- 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> + <span class="navbar-brand"> + <a class="navbar-link" href="../index.html">chemCal</a> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> + </span> + </div> + + <div id="navbar" class="navbar-collapse collapse"> + <ul class="nav navbar-nav"> + <li> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> +</li> +<li> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</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>Toluene amounts by GC/MS as reported by Rocke and Lorenzato (1995)</h1> + + <div class="hidden name"><code>rl95_toluene.Rd</code></div> + </div> + + <div class="ref-description"> + + <p>Dataset reproduced from Table 4 in Rocke and Lorenzato (1995).</p> + + </div> + + + <h2 class="hasAnchor" id="format"><a class="anchor" href="#format"></a>Format</h2> + + <p>A dataframe containing four replicate observations for each + of the six calibration standards.</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> + + + </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="#format">Format</a></li> + + <li><a href="#source">Source</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://pkgdown.r-lib.org/">pkgdown</a>.</p> +</div> + + </footer> + </div> + + + + </body> +</html> + diff --git a/docs/reference/utstats14.html b/docs/reference/utstats14.html new file mode 100644 index 0000000..c16955f --- /dev/null +++ b/docs/reference/utstats14.html @@ -0,0 +1,154 @@ +<!-- 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>Example data for calibration with replicates from University of Toronto — utstats14 • chemCal</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"> + +<!-- clipboard.js --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/1.7.1/clipboard.min.js" integrity="sha384-cV+rhyOuRHc9Ub/91rihWcGmMmCXDeksTtCihMupQHSsi8GIIRDG0ThDc3HGQFJ3" crossorigin="anonymous"></script> + +<!-- sticky kit --> +<script src="https://cdnjs.cloudflare.com/ajax/libs/sticky-kit/1.1.3/sticky-kit.min.js" integrity="sha256-c4Rlo1ZozqTPE2RLuvbusY3+SU1pQaJC0TjuhygMipw=" crossorigin="anonymous"></script> + +<!-- pkgdown --> +<link href="../pkgdown.css" rel="stylesheet"> +<script src="../pkgdown.js"></script> + + + +<meta property="og:title" content="Example data for calibration with replicates from University of Toronto — utstats14" /> + +<meta property="og:description" content="Dataset read into R from http://www.chem.utoronto.ca/coursenotes/analsci/stats/files/example14.xls." /> +<meta name="twitter:card" content="summary" /> + + + +<!-- 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> + <span class="navbar-brand"> + <a class="navbar-link" href="../index.html">chemCal</a> + <span class="label label-default" data-toggle="tooltip" data-placement="bottom" title="Released package">0.2.1</span> + </span> + </div> + + <div id="navbar" class="navbar-collapse collapse"> + <ul class="nav navbar-nav"> + <li> + <a href="../index.html"> + <span class="fa fa-home fa-lg"></span> + + </a> +</li> +<li> + <a href="../articles/chemCal.html">Get started</a> +</li> +<li> + <a href="../reference/index.html">Reference</a> +</li> +<li> + <a href="../news/index.html">Changelog</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>Example data for calibration with replicates from University of Toronto</h1> + + <div class="hidden name"><code>utstats14.Rd</code></div> + </div> + + <div class="ref-description"> + + <p>Dataset read into R from <a href='http://www.chem.utoronto.ca/coursenotes/analsci/stats/files/example14.xls'>http://www.chem.utoronto.ca/coursenotes/analsci/stats/files/example14.xls</a>.</p> + + </div> + + + <h2 class="hasAnchor" id="format"><a class="anchor" href="#format"></a>Format</h2> + + <p>A tibble containing three replicate observations of the response for five + calibration concentrations.</p> + + <h2 class="hasAnchor" id="source"><a class="anchor" href="#source"></a>Source</h2> + + <p>David Stone and Jon Ellis (2011) Statistics in Analytical Chemistry. Tutorial website + maintained by the Departments of Chemistry, University of Toronto. + <a href='http://www.chem.utoronto.ca/coursenotes/analsci/stats/index.html'>http://www.chem.utoronto.ca/coursenotes/analsci/stats/index.html</a></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="#format">Format</a></li> + + <li><a href="#source">Source</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://pkgdown.r-lib.org/">pkgdown</a>.</p> +</div> + + </footer> + </div> + + + + </body> +</html> + diff --git a/man/din32645.Rd b/man/din32645.Rd index 12c641a..ffcbaed 100644 --- a/man/din32645.Rd +++ b/man/din32645.Rd @@ -10,19 +10,18 @@ A dataframe containing 10 rows of x and y values. } \examples{ -data(din32645) m <- lm(y ~ x, data = din32645) calplot(m) ## Prediction of x with confidence interval -(prediction <- inverse.predict(m, 3500, alpha = 0.01)) +prediction <- inverse.predict(m, 3500, alpha = 0.01) # This should give 0.07434 according to test data from Dintest, which # was collected from Procontrol 3.1 (isomehr GmbH) in this case -round(prediction$Confidence,5) +round(prediction$Confidence, 5) ## Critical value: -(crit <- lod(m, alpha = 0.01, beta = 0.5)) +crit <- lod(m, alpha = 0.01, beta = 0.5) # According to DIN 32645, we should get 0.07 for the critical value # (decision limit, "Nachweisgrenze") @@ -40,12 +39,12 @@ round(lod.din$x, 2) ## Limit of quantification # This accords to the test data coming with the test data from Dintest again, # except for the last digits of the value cited for Procontrol 3.1 (0.2121) -(loq <- loq(m, alpha = 0.01)) -round(loq$x,4) +loq <- loq(m, alpha = 0.01) +round(loq$x, 4) # A similar value is obtained using the approximation # LQ = 3.04 * LC (Currie 1999, p. 120) -3.04 * lod(m,alpha = 0.01, beta = 0.5)$x +3.04 * lod(m, alpha = 0.01, beta = 0.5)$x } \references{ DIN 32645 (equivalent to ISO 11843), Beuth Verlag, Berlin, 1994 diff --git a/man/inverse.predict.Rd b/man/inverse.predict.Rd index 26ba6b8..373623e 100644 --- a/man/inverse.predict.Rd +++ b/man/inverse.predict.Rd @@ -52,7 +52,10 @@ } \note{ The function was validated with examples 7 and 8 from Massart et al. (1997). -} + Note that the behaviour of inverse.predict changed with chemCal version + 0.2.1. Confidence intervals for x values obtained from calibrations with + replicate measurements did not take the variation about the means into account. + Please refer to the vignette for details.} \references{ Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J., Smeyers-Verbeke, J. (1997) Handbook of Chemometrics and Qualimetrics: Part A, @@ -60,10 +63,26 @@ } \examples{ # This is example 7 from Chapter 8 in Massart et al. (1997) -data(massart97ex1) m <- lm(y ~ x, data = massart97ex1) inverse.predict(m, 15) # 6.1 +- 4.9 inverse.predict(m, 90) # 43.9 +- 4.9 inverse.predict(m, rep(90,5)) # 43.9 +- 3.2 + +# For reproducing the results for replicate standard measurements in example 8, +# we need to do the calibration on the means when using chemCal > 0.2 +weights <- with(massart97ex3, { + yx <- split(y, x) + ybar <- sapply(yx, mean) + s <- round(sapply(yx, sd), digits = 2) + w <- round(1 / (s^2), digits = 3) +}) + +massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean) + +m3.means <- lm(y ~ x, w = weights, data = massart97ex3.means) + +inverse.predict(m3.means, 15, ws = 1.67) # 5.9 +- 2.5 +inverse.predict(m3.means, 90, ws = 0.145) # 44.1 +- 7.9 + } \keyword{manip} @@ -74,7 +74,6 @@ Analytica Chimica Acta 391, 105 - 126. } \examples{ -data(din32645) m <- lm(y ~ x, data = din32645) lod(m) @@ -68,9 +68,7 @@ and therefore not tested. Feedback is welcome. } \examples{ -data(massart97ex3) -attach(massart97ex3) -m <- lm(y ~ x) +m <- lm(y ~ x, data = massart97ex1) loq(m) # We can get better by using replicate measurements diff --git a/man/massart97ex3.Rd b/man/massart97ex3.Rd index efdcf02..d7f8d00 100644 --- a/man/massart97ex3.Rd +++ b/man/massart97ex3.Rd @@ -5,29 +5,32 @@ \description{ Sample dataset from p. 188 to test the package. } -\usage{data(massart97ex3)} +\usage{massart97ex3} \format{ A dataframe containing 6 levels of x values with 5 observations of y for each level. } \examples{ -data(massart97ex3) -attach(massart97ex3) -yx <- split(y, x) -ybar <- sapply(yx, mean) -s <- round(sapply(yx, sd), digits = 2) -w <- round(1 / (s^2), digits = 3) -weights <- w[factor(x)] -m <- lm(y ~ x, w = weights) -calplot(m) +# For reproducing the results for replicate standard measurements in example 8, +# we need to do the calibration on the means when using chemCal > 0.2 +weights <- with(massart97ex3, { + yx <- split(y, x) + ybar <- sapply(yx, mean) + s <- round(sapply(yx, sd), digits = 2) + w <- round(1 / (s^2), digits = 3) +}) + +massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean) + +m3.means <- lm(y ~ x, w = weights, data = massart97ex3.means) # The following concords with the book p. 200 -inverse.predict(m, 15, ws = 1.67) # 5.9 +- 2.5 -inverse.predict(m, 90, ws = 0.145) # 44.1 +- 7.9 +inverse.predict(m3.means, 15, ws = 1.67) # 5.9 +- 2.5 +inverse.predict(m3.means, 90, ws = 0.145) # 44.1 +- 7.9 # The LOD is only calculated for models from unweighted regression # with this version of chemCal -m0 <- lm(y ~ x) +m0 <- lm(y ~ x, data = massart97ex3) lod(m0) # Limit of quantification from unweighted regression @@ -38,7 +41,7 @@ loq(m0) # If we are not using a variance function, we can use the weight from # the above example as a first approximation (x = 15 is close to our # loq approx 14 from above). -loq(m, w.loq = 1.67) +loq(m3.means, w.loq = 1.67) # The weight for the loq should therefore be derived at x = 7.3 instead # of 15, but the graphical procedure of Massart (p. 201) to derive the # variances on which the weights are based is quite inaccurate anyway. diff --git a/man/rl95_toluene.Rd b/man/rl95_toluene.Rd new file mode 100644 index 0000000..1faafdc --- /dev/null +++ b/man/rl95_toluene.Rd @@ -0,0 +1,16 @@ +\name{rl95_toluene} +\docType{data} +\alias{rl95_toluene} +\title{Toluene amounts by GC/MS as reported by Rocke and Lorenzato (1995)} +\description{ + Dataset reproduced from Table 4 in Rocke and Lorenzato (1995). +} +\format{ + A dataframe containing four replicate observations for each + of the six calibration standards. +} +\source{ + Rocke, David M. und Lorenzato, Stefan (1995) A two-component model for + measurement error in analytical chemistry. Technometrics 37(2), 176-184. +} +\keyword{datasets} diff --git a/tests/testthat/test_compare_investr.R b/tests/testthat/test_compare_investr.R new file mode 100644 index 0000000..6191c89 --- /dev/null +++ b/tests/testthat/test_compare_investr.R @@ -0,0 +1,30 @@ +context("Compare with investr::calibrate") + +library(chemCal) +library(investr) + +test_that("Unweighted regressions give same results as investr::calibrate using the Wald method", { + compare_investr <- function(object, y_sample) { + pred_chemCal <- inverse.predict(object, y_sample) + pred_investr <- calibrate(object, y_sample, interval = "Wald") + expect_equivalent(pred_chemCal[["Prediction"]], + pred_investr$estimate) + expect_equivalent(pred_chemCal[["Standard Error"]], + pred_investr$se) + expect_equivalent(pred_chemCal[["Confidence Limits"]][1], + pred_investr$lower) + expect_equivalent(pred_chemCal[["Confidence Limits"]][2], + pred_investr$upper) + } + m_tol <- lm(peak_area ~ amount, data = rl95_toluene) + compare_investr(m_tol, 1000) + + m_din <- lm(y ~ x, din32645) + compare_investr(m_din, 5000) + + m_m1 <- lm(y ~ x, massart97ex1) + compare_investr(m_m1, 15) + + m_m3 <- lm(y ~ x, massart97ex3) + compare_investr(m_m3, 15) +}) diff --git a/tests/testthat/test_inverse.predict.R b/tests/testthat/test_inverse.predict.R index 61484dc..a4e1bfd 100644 --- a/tests/testthat/test_inverse.predict.R +++ b/tests/testthat/test_inverse.predict.R @@ -23,20 +23,23 @@ test_that("Inverse predictions for unweighted regressions are stable", { }) test_that("Inverse predictions for weighted regressions are stable", { - attach(massart97ex3) - yx <- split(y, x) - ybar <- sapply(yx, mean) - s <- round(sapply(yx, sd), digits = 2) - w <- round(1 / (s^2), digits = 3) - weights <- w[factor(x)] - m3 <- lm(y ~ x, w = weights) - - p3.1 <- inverse.predict(m3, 15, ws = 1.67) + weights <- with(massart97ex3, { + yx <- split(y, x) + ybar <- sapply(yx, mean) + s <- round(sapply(yx, sd), digits = 2) + w <- round(1 / (s^2), digits = 3) + }) + + massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean) + + m3.means <- lm(y ~ x, w = weights, data = massart97ex3.means) + + p3.1 <- inverse.predict(m3.means, 15, ws = 1.67) expect_equal(signif(p3.1$Prediction, 7), 5.865367) expect_equal(signif(p3.1$`Standard Error`, 7), 0.8926109) expect_equal(signif(p3.1$Confidence, 7), 2.478285) - p3.2 <- inverse.predict(m3, 90, ws = 0.145) + p3.2 <- inverse.predict(m3.means, 90, ws = 0.145) expect_equal(signif(p3.2$Prediction, 7), 44.06025) expect_equal(signif(p3.2$`Standard Error`, 7), 2.829162) expect_equal(signif(p3.2$Confidence, 7), 7.855012) diff --git a/tests/testthat/test_lod_loq.R b/tests/testthat/test_lod_loq.R index 6ba0ad0..da0f45e 100644 --- a/tests/testthat/test_lod_loq.R +++ b/tests/testthat/test_lod_loq.R @@ -15,7 +15,11 @@ test_that("lod is stable across chemCal versions", { }) test_that("loq is stable across chemCal versions", { - m2 <- lm(y ~ x, data = massart97ex3) + # Actually it was not stable between chemCal <0.2 and + # chemCal > 0.2, so we needed to adapt the test + # to work on a model on the means + massart97ex3_means <- aggregate(y ~ x, massart97ex3, mean) + m2 <- lm(y ~ x, data = massart97ex3_means) loq_1 <- loq(m2) expect_equal(signif(loq_1$x, 7), 13.97764) expect_equal(signif(loq_1$y, 7), 30.6235) diff --git a/tests/testthat/test_massart.R b/tests/testthat/test_massart.R index 791c2e7..6e4ce75 100644 --- a/tests/testthat/test_massart.R +++ b/tests/testthat/test_massart.R @@ -19,22 +19,24 @@ test_that("Inverse predictions for example 1 are correct",{ expect_equal(round(p1.3$Confidence, 1), 3.2) }) +test_that("Inverse predictions for example data 3 are correct when regressing on means",{ + weights <- with(massart97ex3, { + yx <- split(y, x) + ybar <- sapply(yx, mean) + s <- round(sapply(yx, sd), digits = 2) + w <- round(1 / (s^2), digits = 3) + }) -test_that("Inverse predictions for example 3 are correct",{ - attach(massart97ex3) - yx <- split(y, x) - ybar <- sapply(yx, mean) - s <- round(sapply(yx, sd), digits = 2) - w <- round(1 / (s^2), digits = 3) - weights <- w[factor(x)] - m3 <- lm(y ~ x, w = weights) + massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean) + + m3.means <- lm(y ~ x, w = weights, data = massart97ex3.means) # Known values are from the book - p3.1 <- inverse.predict(m3, 15, ws = 1.67) + p3.1 <- inverse.predict(m3.means, 15, ws = 1.67) expect_equal(round(p3.1$Prediction, 1), 5.9) expect_equal(round(p3.1$Confidence, 1), 2.5) - p3.2 <- inverse.predict(m3, 90, ws = 0.145) + p3.2 <- inverse.predict(m3.means, 90, ws = 0.145) expect_equal(round(p3.2$Prediction, 1), 44.1) expect_equal(round(p3.2$Confidence, 1), 7.9) }) diff --git a/vignettes/chemCal.R b/vignettes/chemCal.R index 701db7b..d9015e9 100644 --- a/vignettes/chemCal.R +++ b/vignettes/chemCal.R @@ -1,36 +1,23 @@ -### R code from vignette source '/home/jranke/git/chemCal/vignettes/chemCal.Rnw' - -################################################### -### code chunk number 1: chemCal.Rnw:38-42 -################################################### +## ------------------------------------------------------------------------ library(chemCal) -data(massart97ex3) m0 <- lm(y ~ x, data = massart97ex3) calplot(m0) +## ------------------------------------------------------------------------ +plot(m0, which=3) -################################################### -### code chunk number 2: chemCal.Rnw:49-50 -################################################### -plot(m0,which=3) - +## ---- message = FALSE, echo = TRUE--------------------------------------- +weights <- with(massart97ex3, { + yx <- split(y, x) + ybar <- sapply(yx, mean) + s <- round(sapply(yx, sd), digits = 2) + w <- round(1 / (s^2), digits = 3) +}) +massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean) -################################################### -### code chunk number 3: chemCal.Rnw:56-63 -################################################### -attach(massart97ex3) -yx <- split(y, x) -ybar <- sapply(yx, mean) -s <- round(sapply(yx, sd), digits = 2) -w <- round(1 / (s^2), digits = 3) -weights <- w[factor(x)] -m <- lm(y ~ x, w = weights) +m <- lm(y ~ x, w = weights, data = massart97ex3.means) - -################################################### -### code chunk number 4: chemCal.Rnw:69-71 -################################################### +## ------------------------------------------------------------------------ inverse.predict(m, 15, ws=1.67) inverse.predict(m, 90, ws = 0.145) - diff --git a/vignettes/chemCal.Rmd b/vignettes/chemCal.Rmd index 2515abb..ccbbdc7 100644 --- a/vignettes/chemCal.Rmd +++ b/vignettes/chemCal.Rmd @@ -8,19 +8,14 @@ output: toc_float: true code_folding: show fig_retina: null -bibliography: refs.bib +bibliography: references.bib vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{Introduction to chemCal} + %\VignetteEncoding{UTF-8} --- -[Wissenschaftlicher Berater, Kronacher Str. 12, 79639 Grenzach-Wyhlen, Germany](http://www.jrwb.de)<br /> - -```{r, include = FALSE} -require(knitr) -opts_chunk$set(engine='R', tidy=FALSE) -``` -# Basic calibration functions for analytical chemistry +# Basic calibration functions The `chemCal` package was first designed in the course of a lecture and lab course on "Analytics of Organic Trace Contaminants" at the University of Bremen @@ -28,13 +23,17 @@ from October to December 2004. In the fall 2005, an email exchange with Ron Wehrens led to the belief that it would be desirable to implement the inverse prediction method given in @massart97 since it also covers the case of weighted regression. Studies of the IUPAC orange book and of DIN 32645 -as well as publications by Currie and the Analytical Method Committee of the -Royal Society of Chemistry and a nice paper by Castillo and Castells provided -further understanding of the matter. - -At the moment, the package consists of four functions, working on univariate -linear models of class `lm` or `rlm`, plus two datasets for -validation. +(equivalent to ISO 11843), publications by @currie97 and the Analytical +Method Committee of the Royal Society of Chemistry [@amc89] and a nice paper by +Castells and Castillo [@castells00] provided some further understanding of the matter. + +At the moment, the package consists of four functions +([calplot](https://pkgdown.jrwb.de/chemCal/reference/calplot.lm.html), +[lod](https://pkgdown.jrwb.de/chemCal/reference/lod.html), +[loq](https://pkgdown.jrwb.de/chemCal/reference/loq.html) and +[inverse.predict](https://pkgdown.jrwb.de/chemCal/reference/inverse.predict.html)), +working on univariate linear models of class `lm` or `rlm`, plus several +datasets for validation. A [bug report](http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877) and the following e-mail exchange on the r-devel mailing list about @@ -42,8 +41,21 @@ prediction intervals from weighted regression entailed some further studies on this subject. However, I did not encounter any proof or explanation of the formula cited below yet, so I can't really confirm that Massart's method is correct. +In fact, in June 2018 I was made aware of the fact that the inverse prediction +method implemented in chemCal version 0.1.37 and before did not take the +variance of replicate calibration standards about their means into account, nor +the number of replicates when calculating the degrees of freedom. Thanks to +PhD student Anna Burniol Figols for reporting this issue! + +As a consequence, I rewrote `inverse.predict` not to automatically work with +the mean responses for each calibration standard any more. The example +calculations from @massart97 can still be reproduced when the regression model +is calculated using the means of the calibration data as shown below. + +# Usage + When calibrating an analytical method, the first task is to generate a suitable -model. If we want to use the `chemCal` functions, we will have to restrict +model. If we want to use the `chemCal` functions, we have to restrict ourselves to univariate, possibly weighted, linear regression so far. Once such a model has been created, the calibration can be graphically @@ -64,16 +76,21 @@ plot(m0, which=3) ``` Therefore, in Example 8 in @massart97, weighted regression -is proposed which can be reproduced by +is proposed which can be reproduced by the following code. +Note that we are building the model on the mean values for +each standard in order to be able to reproduce the results +given in the book with the current version of chemCal. ```{r, message = FALSE, echo = TRUE} -attach(massart97ex3) -yx <- split(y, x) -ybar <- sapply(yx, mean) -s <- round(sapply(yx, sd), digits = 2) -w <- round(1 / (s^2), digits = 3) -weights <- w[factor(x)] -m <- lm(y ~ x, w = weights) +weights <- with(massart97ex3, { + yx <- split(y, x) + ybar <- sapply(yx, mean) + s <- round(sapply(yx, sd), digits = 2) + w <- round(1 / (s^2), digits = 3) +}) +massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean) + +m <- lm(y ~ x, w = weights, data = massart97ex3.means) ``` If we now want to predict a new x value from measured y values, @@ -88,7 +105,7 @@ The weight `ws` assigned to the measured y value has to be given by the user in the case of weighted regression, or alternatively, the approximate variance `var.s` at this location. -# Some theory for `inverse.predict` +# Background for `inverse.predict` Equation 8.28 in @massart97 gives a general equation for predicting the standard error $s_{\hat{x_s}}$ for an $x$ value predicted from measurements of @@ -107,20 +124,32 @@ with s_e = \sqrt{ \frac{\sum w_i (y_i - \hat{y_i})^2}{n - 2}} \end{equation} -where $w_i$ is the weight for calibration standard $i$, $y_i$ is the $y$ -value observed for standard $i$, $\hat{y_i}$ is the estimated value for -standard $i$, $n$ is the number of calibration samples, $w_s$ is the weight -attributed to the sample $s$, $m$ is the number of replicate measurements of -sample $s$, $\bar{y_s}$ is the mean response for the sample, -$\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}$ is the weighted mean of responses -$y_i$, and $x_i$ is the given $x$ value for standard $i$. +In chemCal version before 0.2, I interpreted $w_i$ to be the weight for +calibration standard $i$, $y_i$ to be the mean value observed for standard $i$, +and $n$ to be the number of calibration standards. With this implementation +I was able to reproduce the examples given in the book. However, as noted above, +I was made aware of the fact that this way of calculation does not take the +variation of the y values about the means into account. Furthermore, I noticed +that for the case of unweighted linear calibration with replicate standards, +`inverse.predict` produced different results than `calibrate` from the +`investr` package when using the Wald method. + +Both issues are now addressed in chemCal starting from version 0.2.1. Here, +$y_i$ is calibration measurement $i$, $\hat{y_i}$ is the estimated value for +calibration measurement $i$ and $n$ is the total number of calibration +measurements. + +$w_s$ is the weight attributed to the sample $s$, $m$ is the number of +replicate measurements of sample $s$, $\bar{y_s}$ is the mean response for the +sample, $\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}$ is the weighted mean of +responses $y_i$, and $x_i$ is the given $x$ value for standard $i$. The weight $w_s$ for the sample should be estimated or calculated in accordance to the weights used in the linear regression. -I adjusted the above equation in order to be able to take a different +I had also adjusted the above equation in order to be able to take a different precisions in standards and samples into account. In analogy to Equation 8.26 -from \cite{massart97} we get +from \cite{massart97} I am using \begin{equation} s_{\hat{x_s}} = \frac{1}{b_1} \sqrt{\frac{{s_s}^2}{w_s m} + diff --git a/vignettes/chemCal.html b/vignettes/chemCal.html index 222ca8f..5779746 100644 --- a/vignettes/chemCal.html +++ b/vignettes/chemCal.html @@ -11,7 +11,7 @@ <meta name="author" content="Johannes Ranke" /> -<meta name="date" content="2018-07-05" /> +<meta name="date" content="2018-07-17" /> <title>Introduction to chemCal</title> @@ -236,18 +236,22 @@ div.tocify { <h1 class="title toc-ignore">Introduction to chemCal</h1> <h4 class="author"><em>Johannes Ranke</em></h4> -<h4 class="date"><em>2018-07-05</em></h4> +<h4 class="date"><em>2018-07-17</em></h4> </div> -<p><a href="http://www.jrwb.de">Wissenschaftlicher Berater, Kronacher Str. 12, 79639 Grenzach-Wyhlen, Germany</a><br /></p> -<div id="basic-calibration-functions-for-analytical-chemistry" class="section level1"> -<h1>Basic calibration functions for analytical chemistry</h1> -<p>The <code>chemCal</code> package was first designed in the course of a lecture and lab course on “Analytics of Organic Trace Contaminants” at the University of Bremen from October to December 2004. In the fall 2005, an email exchange with Ron Wehrens led to the belief that it would be desirable to implement the inverse prediction method given in <span class="citation">Massart et al. (1997)</span> since it also covers the case of weighted regression. Studies of the IUPAC orange book and of DIN 32645 as well as publications by Currie and the Analytical Method Committee of the Royal Society of Chemistry and a nice paper by Castillo and Castells provided further understanding of the matter.</p> -<p>At the moment, the package consists of four functions, working on univariate linear models of class <code>lm</code> or <code>rlm</code>, plus two datasets for validation.</p> +<div id="basic-calibration-functions" class="section level1"> +<h1>Basic calibration functions</h1> +<p>The <code>chemCal</code> package was first designed in the course of a lecture and lab course on “Analytics of Organic Trace Contaminants” at the University of Bremen from October to December 2004. In the fall 2005, an email exchange with Ron Wehrens led to the belief that it would be desirable to implement the inverse prediction method given in <span class="citation">Massart et al. (1997)</span> since it also covers the case of weighted regression. Studies of the IUPAC orange book and of DIN 32645 (equivalent to ISO 11843), publications by <span class="citation">Currie (1997)</span> and the Analytical Method Committee of the Royal Society of Chemistry <span class="citation">(Analytical Methods Committee 1989)</span> and a nice paper by Castells and Castillo <span class="citation">(Castells and Castillo 2000)</span> provided some further understanding of the matter.</p> +<p>At the moment, the package consists of four functions (<a href="https://pkgdown.jrwb.de/chemCal/reference/calplot.lm.html">calplot</a>, <a href="https://pkgdown.jrwb.de/chemCal/reference/lod.html">lod</a>, <a href="https://pkgdown.jrwb.de/chemCal/reference/loq.html">loq</a> and <a href="https://pkgdown.jrwb.de/chemCal/reference/inverse.predict.html">inverse.predict</a>), working on univariate linear models of class <code>lm</code> or <code>rlm</code>, plus several datasets for validation.</p> <p>A <a href="http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877">bug report</a> and the following e-mail exchange on the r-devel mailing list about prediction intervals from weighted regression entailed some further studies on this subject. However, I did not encounter any proof or explanation of the formula cited below yet, so I can’t really confirm that Massart’s method is correct.</p> -<p>When calibrating an analytical method, the first task is to generate a suitable model. If we want to use the <code>chemCal</code> functions, we will have to restrict ourselves to univariate, possibly weighted, linear regression so far.</p> +<p>In fact, in June 2018 I was made aware of the fact that the inverse prediction method implemented in chemCal version 0.1.37 and before did not take the variance of replicate calibration standards about their means into account, nor the number of replicates when calculating the degrees of freedom. Thanks to PhD student Anna Burniol Figols for reporting this issue!</p> +<p>As a consequence, I rewrote <code>inverse.predict</code> not to automatically work with the mean responses for each calibration standard any more. The example calculations from <span class="citation">Massart et al. (1997)</span> can still be reproduced when the regression model is calculated using the means of the calibration data as shown below.</p> +</div> +<div id="usage" class="section level1"> +<h1>Usage</h1> +<p>When calibrating an analytical method, the first task is to generate a suitable model. If we want to use the <code>chemCal</code> functions, we have to restrict ourselves to univariate, possibly weighted, linear regression so far.</p> <p>Once such a model has been created, the calibration can be graphically shown by using the <code>calplot</code> function:</p> <pre class="r"><code>library(chemCal) m0 <- lm(y ~ x, data = massart97ex3) @@ -256,14 +260,16 @@ calplot(m0)</code></pre> <p>As we can see, the scatter increases with increasing x. This is also illustrated by one of the diagnostic plots for linear models provided by R:</p> <pre class="r"><code>plot(m0, which=3)</code></pre> <p><img src="data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAqAAAAHgCAIAAAD17khjAAAACXBIWXMAAA7DAAAOwwHHb6hkAAAgAElEQVR4nOzdd1iTZ9sG8CtMUYSAiMheIgKCyBQQXhQRFNSqVG0RcUKtuK0LcW+qWK2rrbMt7gWOKs5aoILKUIZQVJYLGaISWfn+iF9KkeAieZJw/g6O9wh3xnOifT151n2zuFwuAQAAgHSRYToAAAAAtDwUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDyAGElNTf3iiy/09fXbtWtnZWW1cOHC8vLyz/lADQ2N7du3f+y7vLy83NzcPme7AMA4FDyAuLh165ajo+Pz589XrVp18ODBUaNG7dy5s3///jU1NUxHE5bAwMAhQ4YQ0cuXL1ksVnBwMNOJAKSHHNMBAOCtZcuWmZqaxsXFKSgoEJGfn1///v0dHBxOnz7Na0EpJisrO2zYMHt7e6aDAEgP7MEDiIuMjAwbGxteu/PY29t/8803cnJvfxGvra1dtGhR165d27dv7+Hh8ffff/PGKyoqJk2apK2traioaGxsvHz5ci6X2+Qmdu/ebWdnxzv+v3v37k8IWVtbu2zZsm7durVv397JyenUqVMNn/qoeM7Ozr/99tvJkydZLNabN28SExNra2vfuxVdXd3o6OgVK1ZYWVmpqqoOHz78+fPnn/CDAEg/LgCIh4EDB7Zt23bbtm2vXr1q8gVjx45VVVXdsmVLbGyst7c3m80uLi7mcrlBQUEdO3Zcu3btyZMnZ86cSUTR0dG8t3To0GHbtm28x1FRUXJycgsXLoyNjf3222+JaOvWrU1uqG/fvq6urk0+FRISoqiouHr16piYmLFjxxLR8ePHPy1eSUnJ0KFD+/fv/+jRo7q6Oh0dnY0bN753Kzo6OnZ2dlFRURwOJzs7u1OnTqGhoR//hw0g/VDwAOIiMzPTxMSEiNq3bz9o0KDNmzfn5OTwn83KymKxWIcPH+Z9W1lZ2bZt259++onL5Q4ZMmTv3r38V1pbW8+ZM4f3mF/wlZWVampqS5cu5b9s4sSJurq6TSYRVPAPHz6UlZX94Ycf+CMDBw60sbH55Hhff/314MGDeY/5Bd/MVngv69u3L/+pkJCQXr16NflTALRyOAcPIC7Mzc2zsrIuX758/vz5uLi4sLAwIhoyZMiePXtUVVUTExNlZGQGDx7Me7GysvLDhw8VFRWJ6Pjx40TE5XILCwsvXbqUmZnZr1+/Rh9+9+7dsrIyLy+vkpIS3kifPn1++uknDofTpk2bD0yYkpJSV1f35Zdf8ke+/PLL4OBg3gH2z4n3gVvhfaCrqyv/KWVl5fr6+g/MD9Cq4Bw8gFioqanhVW+/fv3Wr19/+/bt4uLiLVu2XLhwYdasWURUUFCgoaEhLy/Pf4uGhkb79u2J6Pbt2/3799fU1Ozevfvvv/+upqb27uc/ePCAiFxdXTv+v1GjRvE+Ni4ujvX/5s2b10zI4uJiFovVsWNH/kjnzp25XO7jx48/M94HboX37Xs/AQAIV9EDiImCggITE5M//vjD29ubN9K5c+dvv/02OTn58uXLRKSlpVVaWlpfXy8j8/b38vz8fHl5eSUlJRcXl8GDB586dcrJyUlGRsbR0fHdz+f1ZUFBga6ubqOntLW1MzMzeY/V1dWbCckr2pKSEk1NTd7IkydPiEhTU/Mz433gVpp/IwA0hD14ALFgaGjIP1/OH6ytrU1NTTU2NiYiOzu7mpqaM2fO8J6qrq52dHTct29fUlISh8OJjIzs1auXjIxMZWVldnb2u59vZWUlLy9/+vRp/khUVNSXX37J5XLbtWtn/v+aL1EbGxsZGZkjR47wRw4fPmxhYaGkpPSZ8T5wK82/EQAawh48gFiQkZHZtm3bl19+6ezsPGLECG1t7adPn0ZHR9+5c+fq1atEZGtr+8UXXwQHB69evdrIyGjnzp2vX78eMWIEb6d51apVgYGBpaWlK1eurKurS0hIyMrKMjc353++pqbm9OnTZ8+eXVlZaW1tHR8fv2LFitWrV7NYrCbzPH/+/MSJEw1HunTpYmlpOX78+NmzZ1dVVVlaWh47duzUqVNHjx795Hjy8vI5OTnJyck9evTgb8jQ0FDQVgDgIzB3fR8ANHbt2jVfX19dXV1FRUUTE5NRo0alpaXxn62qqpoxY4aRkZGysrK7u3tCQgJvfN++fSYmJm3btnVycjp27NilS5fs7Ox4V6Q3vE2urq5u/fr1vF1hc3PzH3/8sb6+vskYffv2ffffirlz53K53Orq6oiICDMzs3bt2jk4OJw4ceJz4l26dMnQ0FBZWbm8vLzhbXLNbKXhy7hc7qJFi7755pvP/XMHkEYsroAJMQAAAEBy4Rw8AACAFELBAwAASCEUPAAAgBRCwQMAAEghFDwAAIAUQsEDAABIIRQ8AACAFGpFM9nNmjWLNyMYAACAaCgpKR05cqRTp06i33QrKvj4+PgpU6Z0796d6SAAANBajBo16smTJyh4oTM3N7ezs2M6BQAAtBYMLpKEc/AAAABSCAUPAAAghVDwAAAAUggFDwDQ2rm7u7P+n5+fH9NxhIvL5To6OmZlZfFHjhw50qVLFzab7enpmZmZyWC2loWCBwBoRZ48eVJVVdVoMC8vLzc3t7KysrKy8vDhw4wEE42rV69OmjQpKSmJP/Lo0aNx48bt2bPn+fPnPj4+w4cPl5pV1FHwAADSr76+/scff9TS0rKxsenYsaOnp2dKSgrvKQ6HU15ebmxsrKysrKyszOBV3yKQmJgoKyvbpk2bhiP29vaurq6ysrIzZszIzMwsLy/nP3vlypX58+d/++23O3bsePcXIzGHggcAkH4LFy6Mjo6+cuXK48ePKyoqgoODfXx8MjIyiOj+/ftycnL29vZsNrtv3745OTlMhxWiuXPnbt++XVVVlT/Sr18//kGL69evGxoastlsIqqtrf3666+/+eYbFRUVS0vLS5cudevWLT09nZncn6R13QcPANAKlZaWbt++PTc3t0OHDkQkKys7ZsyYsrKyVatW/frrrxUVFY6Ojj/++KOBgUF4ePioUaOSk5OZjiw6vOMWXC735MmTISEhW7duZbFYRLRly5Znz56lpaXJy8sT0eTJk3/77bevvvoqLS2N9wLxhz14AAApl5KSYmtry2t3Pl9fX96paGdn5/Pnz3fp0kVBQWHp0qW3bt16/vw5Q0mZUVpaGhAQMG/evKNHjw4bNow3eODAgfDwcF6783z99dc1NTV37txhKOZHQ8EDAEg5GRmZ2traRoO1tbUyMjJE9Ndff126dIn/SllZWUVFRVFHZE51dbWPj4+mpmZaWpqbmxt//MmTJ/r6+o1ebGBg8OjRI9EG/HTiW/Dl5eU5OTl1dXVMBwEAkGw9e/a8c+dOYWFhw8Fjx4717t2biDgczogRIzIyMurq6tasWePp6amsrMxQUgacPHmyrq5uw4YN9fX1HA6Hw+HwrqLX09O7d+9ew1dyudycnJx3W19siUvBT5gw4cqVK7zHz549GzRokJqampmZmbq6+tatW6XmpgUAANFTUVFZsGCBj4/P5cuXa2pqSktL16xZs23btoULFxJR3759v/vuu759+2ppad2+fXvPnj1M5xWpmzdv3rp1S6mBiooKIhozZkxERMSLFy/4r9ywYYOWlpa5uTlzYT+OuFxk98svv9jb2//vf/8jokmTJv31118bN240Nzf/66+/pk+f3rFjx4CAAKYzAgBIqtmzZxsZGc2cOTM7O7tNmzY+Pj4JCQkGBga8Z+fMmTNnzhxmE4rS48eP+Y/XrFmzZs2ad18zbty4zMxMCwuLYcOGsdnsK1euvHjx4vjx4yKM+bnEpeD5SktLT5w4cerUKX9/fyLy8fGpq6uLjIxEwQMAfI5hw4YNGzastrZWTk7s/uUXQywWKzIycuzYsXFxcRUVFbNmzRo4cKCsrCzTuT6C2P01379/n4g8PT35I25ubj/88ANziQAApAfa/aNYWlpaWloyneITics5eD5DQ0M5ObmGF4Okp6c3nJQAAAAA3kuMfpVbvHjxsWPHzMzMjI2Np0+ffu7cOSKKi4vbvHmzj48P0+kAAAAkibgU/JEjR3Jzc3Nzc+/evfv69eurV6/yxv38/KysrFatWsVsPAAAAMkiLgXPnzyIh8Ph8B7ExsZ6enpK1nUNAAAAjBO7c/A8/KV+vLy8ampqXr58yWweAAAAySIue/DNCA4OPnjw4IfMdVNWVhYXFyfo2UePHkncYn8AAACfRgIK/sPnTXz48CF/1b93FRUVpaSkNLwBDwAAQFpJQMGHhISEhIR8yCt79Ohx6NAhQc+2b9++ffv2LZcLAABAfInpOXgAAAD4HCh4AAAAKYSCBwAAkELicg5+9uzZzb8gMjJSNEkAAACkgLgUfG1t7c6dO6uqqnR1dRUVFd99AQoeAADgw4lLwUdFRfXr18/Pz+/s2bNWVlZMxwEAAJBsYnQO3tfXt23btkynAAAAkAZiVPAyMjL79+/X0dFhOggAAIDEE5dD9DxDhw5lOgIAAIA0EKM9eAAAAGgpKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAgNbuyJEjXbp0YbPZnp6emZmZvMHo6OguXbpoamoGBQVVVVUxmxA+AQoeAKBVe/To0bhx4/bs2fP8+XMfH5/hw4dzudz09PSwsLCjR48+fPiwrq5uyZIlTMeEj4aCBwBoFZKTk4cNG2Zqampvb79w4cLKykreeGJior29vaurq6ys7IwZMzIzM8vLy8+ePTtkyBBra2slJaWlS5dGR0czGx4+AQoeAED6/frrr0OGDBk4cOAff/yxY8eOkpKSnj17lpSUEFG/fv0OHz7Me9n169cNDQ3ZbHZNTQ3/vbKysoWFhbW1tcxEZ1plZWVxcTHTKT4FCh4AQMpxOJyZM2f+8ccf48aNMzExsbOz27Fjx8CBA1evXk1EysrKHTp04HK5J06cGDVq1Pr161ksVr9+/Y4fP56cnFxaWrpgwQIul/vixQumfw5RS05OdnV11dHRcXBw6NixY1RUVF1dHdOhPgIKHgBAyt28edPExMTS0rLhYHBw8IULF3iPS0tLAwIC5s2bd/To0WHDhhGRo6Pj+vXrg4KCevXq1bNnTzk5OTabzUB05qSlpfn5+YWFhVVUVBQVFcXHxx87dmzBggVM5/oIKHgAACn3+vVrZWXlRoPt27d//fo1EVVXV/v4+Ghqaqalpbm5ufGeffny5YABAzIyMrKzsz08PLp27Soj07r6YtWqVeHh4SNHjmSxWETUpUuX48eP79y5s7S0lOloH6p1/YUBALRClpaWt2/ffvXqVcPBa9euWVtbE9HJkyfr6uo2bNhQX1/P4XA4HA6Xy33w4EHXrl2zs7NfvXq1dOnSiRMnMpSdMTdv3uzfv3/DkQ4dOtja2qampjIV6WOh4AEApJy2tvagQYPGjBnz/Plz3si1a9fCw8PnzJlDRDdv3rx165ZSAxUVFVZWVitXrnR3d7ewsLC2tg4LC2P0J2CAnJxcwysNeaqrq2VlZRnJ8wlQ8AAA0m/r1q2GhoZdunRxcXHp1q3buHHjfv755169ehHRmjVruP/FO90+ZcqUJ0+ePHz4cPXq1a3t+DwRubm5HTlypOHIw4cPMzMz7ezsmIr0seSYDgAAAELXpk2byMjIRYsWZWRkqKmpmZqaysnh3//mhIeH9+rVS1ZWduLEiaqqqleuXJk2bdrixYvbtWvHdLQPhb9gAIDWQlVVlbfXDu9lYGBw48aNBQsWWFlZVVZW2tjYrF+/3t/fn+lcHwEFDwAA0ARdXd19+/YxneLTtbrTKgAAAK0BCh4AAEAKoeABAACkkAwRubm5/fPPP0wnAQAAgBYjR0SVlZW8GY5KSkoWLVpUX1/f/Ht69uwZEhIiinQAAADwSf5zFX27du3s7OzeW/BmZmbCjAQAAACf6z8Fr6SkNGHChCZfl5OTc+/ePUdHx44dO4okGAAAAHw6gRfZ5efne3t78+YfPnv2bLdu3fz8/MzMzJKSkkQYDwAAAD6FwIKfOnVqenq6u7s7Ea1cudLZ2Tk3N9fFxSUiIkKE8QAAAOBTCCz4a9euLVmyJCAgoLy8PD4+/ttvvzUxMRk9enRycrIo8wEAAMAnEFjw9fX1KioqRHT+/Hkul9unTx8ikpWVraqqEl06AAAA+CQC56J3dHTcsWOHgYHBunXrXFxcOnXq9OLFi/379+MSegAAAPEnsODXr1/fv39/V1dXJSWls2fPEpGrq2tWVtahQ4dEGA8AAAA+hcCCt7GxycvLy8jI0NfX19TUJKLFixdbWlp269ZNhPEAAKBllJWV7dy5My0tTUVFxcfHZ/DgwUwnAuFqbi76tm3b2tvb89qdiIYPH452BwCQRImJiZaWlvfv3x84cKC1tfXy5cv9/f2rq6uZzgVC1HgPfvv27e99T2hoqHDCAABAy6uvrw8MDNy5c6efnx9vJCQkZPDgwZs2bZozZw6z2UB4Ghd8eHj4e9+DggcAkCApKSlKSkr8diciGRmZuXPnzpo1CwUvxRoXfElJCSM5AABASJ49e6atrd1oUEdH59mzZ4zkAdH4uPXgMzMzo6KihBQFAACEwcjIKCsrq9FCYnfv3jUyMmIqEoiAwKvouVzugQMHbt++3fC/ifj4+JqamunTp4skGwAAtAAzMzNDQ8Ply5dHRESwWCwievLkyYIFCxYuXMh0NBAigQW/fPnyJUuWdO/ePTMzU01NzcDAgPcLYGxsrCjzAQDA54uOjh45cuSJEyfc3d3Ly8vPnTs3derUESNGMJ0LhEjgIfrdu3fPnDkzNTV1586dtra2N27cKCoqsrKyKi8vF2U+AAD4fNra2levXo2KijIxMenbt++NGzew+y71BO7BFxcX85aS69u374wZM+rr69u3b79gwYLVq1cPGTJEhAkBAKAFsFgsDw8PDw8PpoOAiAjcg1dXV8/OziYiXV1dWVlZ3jLwGhoa6enpoksHAAAAn0Rgwfv5+UVGRu7bt4/FYtnZ2UVFRRUWFu7fv//dey0AAACkVk0NVVQwHeJTCCz4yMhIHx+f48ePE9GaNWtOnDihp6f3008/RUREiDAeAAAAE8rK6PffaeRI6tSJ1q1jOs2nEHgOXlVVde/evbzHtra2jx49SkpKMjExMTY2FlU2AAAA0bp/ny5coJgY+vNPcnAgPz/asIEk89C1wIJvhM1m9+vXT6hRAAAAGFBfT7dvU0wMxcZSfj75+FBQEEVHk7Iy08k+i8CCbzhrcUMWFhbrJPNgBQAAwL+qqigujmJjKSaG1NTI35+iosjVlVgsppO1DIEFr6GhwX/M5XKfPn2akJDQqVOnoUOHiiQYAACAEDx7RmfPUmwsnT9PVlbk709XrpCZGdOxWp7Agt+zZ0+jkYqKii+++KK2tla4iQAAAFrc3btvd9azsqhPH/Lzo59+IlVVpmMJ0UcsNqOqqrpo0aItW7YILw0AAECLqa2l69dp2jTS16dBg6i4mJYsoUeP6NAhCgqS7nanD7/Ijuf169d5eXlCigIAANACSkvp4sW3F80ZG5OfHx0/TnZ2TMcSNYEF/+uvvzYaef78+ffff29vby/kSAAAAB8vL+9tqSclUe/e5O9P69aRlhbTsRgjsOAnTJjQaERGRsbGxmbbtm1CjgTQGuXn5/MmhHZ0dNTT02M6DoCEqKujlBSKiaHDh6mkhPr3p6lTydubFBWZTsY8gQXP4XBEmQOg1aqtrf3uu+9+//13Nzc3Ipo8eXJgYOC6detkZWWZjgYgrl6/posXKTaWTp0idXXy96cdO6TpDrcW8XHn4AGgxS1btiwzMzMrK4vNZhNRWVnZyJEjly9fvmTJEqajAYiZ/Hw6d+4/08wtWkS6ukzHElONC77h7e9NcnJyOn36tNDyALQu9fX127ZtS05O5rU7EampqW3fvt3Z2Xnx4sUs7I4AUIM73LKzydOTAgLot99IRYXpWOKuccGvWLGC94DD4YSHh2toaAwfPlxHR+fx48dHjx598+bN/PnzRR4SQGqVlZVxuVwDA4OGg0ZGRnV1dWVlZerq6kwFA2AYh0PXr1NMDB07RgoK5OdHS5bQ//5Hcjjw/KEa/0mFhobyHzg6Op47d05BQYE3snLlSj8/v/379/POFALA51NWVq6qquJwOG3atOEPVlVVVVVVKUv4PNgAn6KkhM6codhYunCBLC3J358uXCBzc6ZjSSSBE93ExMRMnjyZ3+5EJCcnFxoaGhMTI5JgAK2CoqKip6fn9u3bGw5u27atb9++Df/fB9AisrOzDx8+fOHChfLycqaz/FdeHm3aRP36kakpHT5MXl6UnU3Xr9PcuWj3TybwWAeXyy0qKmo0mJ+fL4fDIwAt6ocffvD09Lxz584XX3xBRMeOHbt48eLly5eZzgVSpaysLCQkJDEx0dnZuby8PC0tbcWKFe/eDi1SdXWUkECxsXTyJHE45O1NU6dS//6EX21biMC29vX1jYiIMDU1HTBgAIvF4nK5Z86ciYiIGDFihLAz1dbWlpaWqqmpycvLC3tbAIwzNja+e/fuDz/8wJtkwtXVddOmTTg+Dy0rKCjI0NAwJydHUVGRiO7du+fr66utrT1gwABRR3n1ii5dosOH6fRpMjIiPz/69Vfq2RN3uLU4gQW/cePGnJwcPz8/NTU1bW3t4uLisrIyDw+PyMhIIUVJSUnZsmVLbGzs06dPuVwui8XS1NQcMGDAtGnTbGxshLRRAHGgrKy8YMECplOA1Lp//35ycvKxY8f4e01mZmZr1qzZsmWL6Ar+wQM6f/4/d7itWUPa2iLaeqsksOBVVFSuXr168eLF5OTkx48f6+joODg4eHh4COm+nT///NPb29vAwCA4ONjQ0FBVVbWioqKwsDAmJsbZ2fnixYsuLi7C2C4AgNTLycmxtrZudEzU3t5+3rx5wt1wfT3dvv12+tj8fPLxoaAgio4myTlAdeXKlT/++OPFixfW1tZBQUFKSkpMJ/oIzZ1QZ7FYXl5eXl5eIsixYMECX1/fQ4cONTrHv2zZspCQkPnz51+9elUEMQAApI+qqmpJSUmjwWfPnvFnX2hhVVX0118UE0NHjhCbTf7+FBVFLi4k8xHrlzKutrZ2zJgxt27dCgoK0tPTu3Tp0urVq2NiYrp37850tA/VRMEHBwfb29s3s+779OnTWzxHamrqnj173r2CT0ZGJjg4mIGzRAAA0qJnz55Pnz69fv16w5ucN2/e7O/v35KbefaMzp6l2Fg6f56srMjfny5fJjOzltyECG3ZsuXZs2dpaWm8Ix+TJ0/+7bffvvrqq7S0NEmZgaqJgt+7d+/Lly8TExMFvUcYBW9iYpKQkDB06NB3n4qPjzc1NW3xLQIAtBLy8vI///zzsGHDvvnmG09Pz7Kysp07d5aWlja6P/MT8aeZy8qiPn3Iz49++kkKllo/cODAunXrGp7X+Prrr5cvX37nzh1J2YlvouArKytlZWVFfKYhLCxs/PjxRUVFw4YNMzQ0VFFRqaysLCgoiImJ2bVr1+7du0UZBgBAyvTv3z8pKWnDhg2LFy9WV1cfPHjw+PHjP/22Z/4dbsePU00N+fvTkiXk4UFSdOvTkydP9PX1Gw0aGBg8evRIggu+yftzcnJy7t275+jo2LFjR2HkGDdunLy8/Lp166KjoxuOW1lZ7du3LzAwUBgbBQBoPfT19aOioj7rI0pL6eLFtxfNGRuTnx/9/jvZ2bVQQPGip6d37949Q0ND/giXy7137967rS+2BP76lp+fP2HChK5du27evPns2bP+/v51dXVsNvv8+fMODg7CiDJ69OjAwMCCgoLi4uLS0lJ1dXVtbW09PT1JOdsBACCd8vLelnpSEvXuTf7+tHYtde7MdCzhCg4OXrRokZOTk+r/n27YsGFD586dzSVnZj2BBT916tT09PSJEycS0cqVK52dnffu3Tt16tSIiIizZ88KKQ2LxdLX12/0+xGHw6mtrcW8HwAAotPwDreCAurfn6ZOJW9vUlRkOpmIjB07NjMz08LCYtiwYWw2+8qVKy9fvjx27BjTuT6CwIK/du3a6tWrAwICysvL4+Pjf/vtNxMTk9GjR4eFhYkyHxEFBwcfPHiQy+W+95WJiYlhYWGCXllVVSV20y8DAIiV16/p4kWKjaVTp0hd/e0dbq6urXCaORaLtX79+uDg4EuXLpWXl8+ZM2fgwIEyEnWnn8CCr6+vV1FRIaLz589zudw+ffoQkaysbFVVlejSERGRp6fnB+6+29jYNHNRqLu7u6rkX9gJANDynj6lc+fo8OF/p5kLDyc9PaZjMc/S0tLS0pLpFJ9IYME7Ojru2LHDwMBg3bp1Li4unTp1evHixf79+81EflNjSEhISEjIh7xSSUnJTvDlHjIyMjidLxrFxcX79++/d++erq5uQECAlZUV04fqH58AACAASURBVIkAoCn8O9yys8nTkwIC6LffSEWF6VjQMgQebVi/fn1WVparq2tGRsaqVauIyNXV9ezZs4sWLRJhPJA8R44csbW1LS4udnV15XK53t7eK1asYDoUAPw/DofOnKHQUNLVpYAAKiujtWvp8WM6dIiCgtDu0kTgHryNjU1eXl5GRoa+vr6mpiYRLV682NLSslu3biKMBxLmyZMnkydPvnTpEn+vferUqY6Ojh4eHr1792Y2G0DrVV1NN27Q5ct0+TIlJZG9Pfn705UrhDnEpFpzsxy0bdvW3t6ew+Hk5uaampoOGzYMh7ihebGxsf379294TF5DQ2PKlCmHDh1CwQOIVF0dpaRQXBxdv05//UWdO5ObG33zDR07RkKagh7ETHMFv23btlWrVhUWFhIRl8sdMWKEo6PjrFmzhFHzs2fPbv4FwlumFlrQkydP9N65MMfAwCA+Pp6RPACtS10dZWXRX39RXBxduEDq6uTlRUFBtHcvqaszHQ5ETWDB79q1a/LkyRMnTuzTp8+oUaOIqFevXrNnz27Xrt0333zT4jlqa2t37txZVVWlq6ur2NR9lih4iaCvr3/ixIlGg1lZWQYGBozkAWgV8vIoLo7i4ujiRWKzycuLAgJo2zbq0IHpZMAoLpdrbW2dmprK/a9u3bpNmTKFy+U+e/aM9zIul7to0SJLS0uucMTGxhJRenq6kD5fWVn5p59+EtKHA09ZWZmWltaFCxf4I/fu3evcuXNKSgqDqQCk0D//cHfs4AYEcDt04BobcydN4u7dyy0oYDoWNNZkw4qGwD34Bw8e+Pr6Nhp0dXUV3p60r69v27ZthfThIBpsNvvQoUOBgYG2trY9evR48ODB2bNn161bZ2Njw3Q0AMnH21O/fp0uXyYFBXJ1JS8viowkyZkdHURJYMF37do1MTGx0ULsqampJiYmQooiIyOzf/9+HR0dIX0+iEbv3r0zMjJOnTp179693r17r127tlOnTkyHApBYeXlvr5I7d45qasjNjby8aNkyarAICkCTBBb8lClTQkND5eTk+vXrR0QFBQVnzpxZtmzZ0qVLhZemyfXgQeK0a9eOd90GAHyKR4/o+nWKi6Pz54nDod69ycuLpk4liZ1SDRghsODHjRtXWVm5bNmyxYsXE5G+vr6CgsL06dNnzJghwngAAK3D48f0559vj8CXlJCHB7m60qRJ1LNnK5wHHlpE0wX/5s2bR48ehYSETJgw4e7duw8ePOjYsaOVlRVvxhsAaFl5eXmLFy9OSEggIhcXlyVLlhgbGzMdCoTv6VO6evXtEfj8fHJyIjc3lDq0lKYLPi8vz8LCYt++faNHj3ZycnJychJxLIDW4++//x40aNB33323ZMkSIjp27FivXr1iY2MdHByYjgZCUFJCCQlvb1XPySFHR/Lyoh07yNaWJGqlMhB/TRe8ubl5nz59oqOjAwMDMXsdgFDNmDFjy5YtAQEBvG/nzJmjp6c3Y8aM69evMxtMzP3888/bt29/+vSpvr7+qlWr3N3dmU4kWGUl/f3321vV790jJyfy8qKoKHJ2JrnmZhsD+BxN/7fFYrHGjx+/bNkyBweH/v37d+jQoeEiuNOnTxdVPAAp9/r169TU1C+++KLh4PDhw8ePH19VVaWkpMRUMDHn5OR069YtHx+f3r1737hxw9PTc+LEic2sFs2AV68oIeHtOfU7d8jBgVxdac0acncnBQWmw0GrIPCXx++++46IXr58uXfv3kZPoeABWsqbN28UFRXl/rsbJycnp6CgwOFwUPBNWrt2bUpKSn5+fufOnXkj58+f9/X1DQ0N7dGjB5PJXr+m+Pi359STksjKitzcaMkS6t2bmpqgE0CoBBY8bwp6ABAqNputrKycmpracC6glJQUNputpqbGYDBxtm/fvkGDBvHbnYi8vb2NjIyioqL27Nkj6jRVVXTz5ttz6gkJZG5OXl40dy65uVGbNqIOA9AATv8AMInFYi1cuDAoKOjAgQO8tZjv3r0bFBS0cOFCpqOJr4qKCsN3pnnR1NR89OiRiBLU1lJq6ttz6g1L3dWVcNAFxAYKHoBhISEhMjIynp6eHTp0IKLnz5+vWLFiwoQJTOcSX507d75582ajwby8vJEjRwpxqw1LPTGRunZ9e5/6kSOkqirE7QJ8KhQ8APMmTpwYHBx87949FotlZmYmhyurmxUeHj506NAjR44MHz6cN/Ldd9+VlJQsWLCghbfEW1Kdd049Lu7tkuqTJtHhw1hSHcQf/h35UNXV1Qq49hWERl5e3hITkX6YwYMHT5s2bcSIEVpaWlpaWvfv33/16tWePXtaZiauJpdUDwig7duxpDpIFsyr8B55eXkBAQG8K56sra0PHTrEdCIAoA0bNuTm5o4YMUJPT2/q1KnPnz8PDAz8rE/My6OdO+nLL0lTkwYNops3KSCAcnLon39oxw4KCEC7g8RpvAevoaHR/BucnJxOnz4ttDzi5Z9//nFzc5szZ87PP//crl27+Pj4b7/99uHDh3PmzGE6GkBrZ2RktGHDhs/6CN7qq3FxdPkyqaiQlxf5+dGGDaSr20IZAZjUuOBXrFjBe8DhcMLDwzU0NIYPH66jo/P48eOjR4++efNm/vz5Ig/JmJUrV06dOnXmzJm8b93d3c+ePWttbR0aGtq+fXtmswHAp2i4pHpd3dvVV7GkOkijxgUfGhrKf+Do6Hju3Dn+ieeVK1f6+fnt37/fzc1NpBmZEx8f32hnXVdX18rK6vbt22I9LyYANFRc/PacesMl1ZcuJSMjppOBWKuqqvr+++8vXLhQXl5uY2Mzb948CwsLpkN9BIHn4GNiYiZPntzwsjI5ObnQ0NCYmBiRBBML9fX1srKyjQZlZGTq6uoYyQMAH+rRIzp8mEJCyNiYevakw4fJzo7OnKHiYjp0iCZNQrtD80pLS+3s7O7du7d8+fL9+/fb29v36dPn2LFjTOf6CAKvoudyuUVFRY0G8/PzW9UNPI6OjmfOnDEzM+OPPH36NC0tzdbWlsFUANC0J0/o2jUsqQ4tYvXq1X369NmyZQvvW2trazc3t4EDB/r5+UnKHVUC29rX1zciIsLU1HTAgAEsFovL5Z45cyYiImLEiBGizMesBQsW9OnTR1VVdfTo0XJycunp6RMnTgwLC2PjFlgAMfHkCV2+TFeu0OXLVFpKHh7k6UnTppFEHUoFMfTHH380WoqlZ8+eurq6t27dcnZ2ZirVRxFY8Bs3bszJyfHz81NTU9PW1i4uLi4rK/Pw8IiMjBRlPmZZWFicP39+5syZYWFhSkpKbdq0CQ8PnzRpEtO5AFqr168pI4PS0+nuXUpPpzt3iMOh3r3J05MmT6bu3bGnDi3l5cuXqu/MUchms1++fMlInk8gsOBVVFSuXr168eLF5OTkx48f6+joODg4eHh4tLbl4a2trePi4t68eVNZWfneewgBoCXV1lJ+Pt29SxkZdPcu3bxJeXlkbEyWlmRhQRMmkIUFWVig1EEYLCwsEhISjI2N+SNVVVUpKSkSdJ1dcyfUWSyWl5eXm5tbYWGhqakpl8ttbe3Op6ioqIjVHgGErbj4bZfz/vfWLWKzyc6OLC3Jy4umTSMrK6y7CqIxa9asMWPGmJub29nZEdGLFy+++eYbHx8fbW1tpqN9qOYKftu2batWreKtG8vlckeMGOHo6Dhr1qxWW/MA0JLKyv7t8owMSkmhujqysCA7O7Kzo9GjydaW2rVjOiW0Up6enhs2bBgyZEinTp3U1NRSUlKGDx/+/fffM53rIwgs+F27dk2ePHnixIl9+vQZNWoUEfXq1Wv27Nnt2rX75ptvRJgQAKTCmzeUm0s3b/7b6KWlZGLyttH9/cnamlpkMnmAFjJ8+HB/f//09PTy8nJra+uWWexAhAQWfGRk5JQpUzZv3lxSUsIbmTFjRllZ2Y8//oiCB4D3qKmhgoK3J855jd7w9PmkSWRpSUZGOH0OYk5RUdHe3p7pFJ9IYME/ePDA19e30aCrq2uruooeAD4U//Q5r9Gzs0lTkywsyNKS/Pxo8WIyN6d3po0CAOERWPBdu3ZNTEwcMGBAw8HU1FQTExPhpwLJxuFwtmzZkpKSYmJiEhoa2rlzZ6YTQUtrePr85k1KTaX27d/unfNWTLezIyUlplMCtGoCC37KlCmhoaFycnL9+vUjooKCgjNnzixbtmzp0qUijAeS5/Dhw0FBQfLy8np6en/88cfKlSvnz5+/fPlypnPBZ3jxgnJy/m305GR68+bt3rmFBfn7U48ehJtIAcSMwIIfN25cZWXlsmXLFi9eTET6+voKCgrTp0+fMWOGCOOBhHn69OmoUaOmTJkSFRXFG9m7d++4ceN69+7t7e3NbDb4UDU1dO/ef25X450+592uNmkSbdpEDW4OBgDxJLDgWSzW9OnTJ0yYcPfu3QcPHnTs2NHKykririEEEVu1apWGhga/3YlozJgxv/zyy+LFi1HwYqqujh4+/M/tallZ1KnT24vbAwJo8WLq1o1kBC5MBQDiSWDBl5SUqKqqKisrOzk5OTk58QYrKysvXrw4ZMgQUcUDCZOZmWn8zr6do6PjwYMHGckDTWg0mczt26Sq+vZgO28yGUtLatOG6ZQA8LkEFnzHjh1tbGwOHTrUcC21oqKiL774gsvliiQbSJ6OHTvev3+/0eDDhw9VVFQYyQNUXk537vzb6GlpVFPz7+nz0aOpRw9SVmY6JQC0vOZmsquurrazs9uxY8dXX30lskAg0cLCwlxcXK5du+bu7s4befLkSUxMDK7NFJHqasrJ+c9kMs+fk6np20b38iIHB9LSYjolAIhCcwUfHR199OjRwMDAy5cvb9q0qW3btiKLBRLKyckpMDDQ09PT09PT1dU1MzPzxIkTlpaWc+fOZTqaNOKvxdLkZDJBQViLBaA1a67gZWVlly1b5ubm9vXXXycmJh46dAiz0MN77d27d9SoUeHh4Tt27NDQ0Pj+++/DwsKYDiUtGk0mk5Hx71osfn40dy7WYgEAvuYKnsfb2zslJWXkyJH29va4Rw4+hI+Pj4+PD9MpJF+jyWTS0khZ+T+TyfTsSTiuBgACvL/giUhHR+fSpUvh4eErV64UdiCA1qWmhsrL//168IDu3KE7dyg9nerqqHt3srKi7t3pyy/J0pLYbKbjAoDEEFjw6enppqam/G/l5eXXrl3r6el5/fp1kQQDkEwNC7uigsrK/tPfDb8qKqi8nN68ITb73y99/bfH262sCFP8AsBnEFjwVlZW7w7i0Cu0RlVVVFb2ni8O5+3LHj8mRUVSU2viy8zs7QMlJWrT5u1jLS3MIQMAwtBEwQcHB9vb29fW1gp6z/Tp04UZCUDIPqqwnzwhBYWmC9vYuAUL+9mzZ7t27SKicePGdezYsaV/ZgBodZoo+L179758+TIxMVHQe1DwIF7eW9j8thZhYX8UOzu7W7duycjIENG8efPs7OySk5OFvVEAkG5NFHxlZaWsrKwSlnoEprRsYTdsa7E8JG5nZ5eSkrJ79+7g4GAi2rVr18SJEx0dHW/cuMF0NACQYE0UvDLmrYSW1coK+6PU1dXdunVr06ZNvHYnonHjxpWVlc2ePZvRXAAg8RoXvMb7FnV2cnI6ffq00PKAhCsooNu3336lplJJCVVVkaoqsdmkpvafy8XZbOrS5e1TDb/at2f6ZxCplJQUIpo6dWrDwVmzZs2ePTslJaVHjx4M5QIAide44FesWMF7wOFwwsPDNTQ0hg8frqOj8/jx46NHj75582b+/PkiDwniisul3Fy6fZtu3Xpb6jIyZGtLPXtSYCB9/z1para2wv5YvCWYX7582fDIWUVFBf8pAIBP07jgQ0ND+Q8cHR3PnTunoKDAG1m5cqWfn9/+/fvd3NxEmhHER20tZWf/O/N5QgIpKJCdHdnZ0ZQpZGlJ76wVC83T09OTlZUdOnTo+fPn+YNDhw6Vk5PT1tZmMBgASDqB98HHxMRs2rSJ3+5EJCcnFxoaOnny5B07dogkG4iBmhq6d49u3nz7xVs7nNfokybRL78Q9jI/2+LFiyMiIrp06TJr1iwiioyM/Oeff1atWsV0LgCQbAILnsvlFhUVNRrMz8+Xk/ug2W1BUlVWUmrqv/voKSlkYPC20QMCsHa4MCxatMja2nrMmDGTJ08mIhUVlZMnTw4aNIjpXAAg2QS2ta+vb0REhKmp6YABA1gsFpfLPXPmTERExIgRI0SZD4SuvJzu3Pl3H5233iiv0UePJjs7wg2Twjd48ODy8nKmUwCAVBFY8Bs3bszJyfHz81NTU9PW1i4uLi4rK/Pw8IiMjBRlPmh5xcX/7qDfvElFRWRlRXZ25OVF06ZR9+7U4LwMAABIqKYL/s2bN6WlpefOnYuPj09OTn78+LGOjo6Dg4OHhweWhJc8vEbnfSUlUXU1WVi8bfS5c8nCgvB3CtAKcDic3377LT09XVVV1dvb29XVlelEIFxNF3xeXp6FhcW+fftGjx7t5eUl4kzwWZq50H3SJNq0CRe6A7RCd+7cGTx4cI8ePdzd3SsqKoKDg93c3H7++WdZWVmmo4GwNF3w5ubmffr0iY6ODgwMxC67uGt0ofutW8Rm/9vou3YRVi4BaN24XO6oUaOWLFkyevRo3sjcuXN9fHy2b9/+7bffMpsNhKfpgmexWOPHj1+2bJmDg0P//v07dOgg02A2UCw2w7DmL3S3taV27ZiOCABiJD09vaamht/uRKSoqBgREbFo0SIUvBQTeJHdd999R0QvX77cu3dvo6daW8G/fPny0KFDeXl5gYGB5ubmDCRo/kJ3e3tq04aBVAAgIR49emRoaNho0MjI6N17oUGaCCz4wsJCUeYQW0FBQb/++isRycjIrFy5snPnzqmpqUJfrpt/WRxvH73hhe5z55K5OeG0GQB8MF1d3dzcXC6X2/CUa05Ojr6+PoOpQNg+bhmuzMzMqKgoIUURQ3PmzNm/f//06dPr6+tra2sTExPLysq6dOnS8lsqLqaYGFqyhPz9qVMnsrKitWuprIz8/OjQISoro+vXadMmCgoiS0u0OwB8FEtLSzabvXXrVv5IZWVlRETEmDFjGEwFwtbcTHYHDhy4fft2fX09fzA+Pr6mpqb1HKLfvHmzu7v7hg0beN86OTndvXvXxMTk0qVLffr0+fTP5V/ozttHv3WLFBX/vSxu507q3LllfgAAACIiOnDgwBdffHH8+PH//e9/5eXlBw8eHDZs2NixY5nOBUIksOCXL1++ZMmS7t27Z2ZmqqmpGRgYZGVl1dfXx8bGijIfs6qrq/mr7/AYGxvLy8sfPHjw4wq+0YXuqamkqfn2ZvSpU8nZGRe6A4BQmZqa3r59+8SJE6mpqZ07dz59+rS1tTXTocRdbW3tjh07Lly4UF5ebmNjM2vWLMk6qSGw4Hfv3j1z5szIyMg9e/YcOHDg3LlzlZWV/fr1a20Tar57EUpdXZ2qqup73vbiBaWl/buPnpWFC90BgFlycnLDhw8fPnw400EkQ2VlZd++fbW0tMaPH89ms69cueLg4LBv377+/fszHe1DCSz44uJid3d3Iurbt++MGTPq6+vbt2+/YMGC1atXDxkyRIQJmaShobF+/frZs2fzRzZu3FhfXz9nzpzGLy0rezvzK++rsJC6dHm7jz5pEi50BwCQLOvWrbO0tNy9ezfvWw8Pj759+3711Vd5eXmSsuiawJTq6urZ2dlEpKurKysrm5SU5OTkpKGhkZ6eLsJ4DDtx4oSbm5u6uvrs2bONjIy2bNkSHx/v7e3dsWPH/8z/mpFBz5/jQncAAKkRGxu7c+fOhiNubm4aGhq3b992cHBgKtVHEVjwfn5+kZGRnTp1CgoKsrOzi4qKWr9+/f79+7W1tUWZj1kuLi4ZGRleXl6LFi3qyuU6yslFurr2qqsjdXVq25ZsbcnWlsaPJ1tbMjBgOiwAALSYiooKDQ2NRoMaGhoSdJ5aYMFHRkZWV1cfP348KChozZo1Li4uBw4ckJWV3bNnjwjjMc/c3LywsJDS0mjECOre/W2p29pSp05MRwMAAGExMzO7efOmkZERf6S6ujotLa1r164MpvooAgteVVWVP4edra3to0ePkpKSTExMjFvnUiXW1pSZyXQIAAAQkbCwsGnTpllbW5uZmRERh8OZMWOGi4uLBF1I/6FXCrDZ7H79+gk1CgAAgJgYOHDg48eP3d3du3fvrqamlpCQ4Orqyr/mTiI0Lvh3Tzk04uTkdPr0aaHlAQAAEAvjx48fPnz4jRs3ysvLV65cKZRpTIWpccGvWLGC94DD4YSHh2toaAwfPlxHR+fx48dHjx598+bN/PnzRR4SAACAAaqqqpJ79LpxwfMnbgsNDXV0dDx37pyCggJvZOXKlX5+fvv373dzcxNBsvLy8mfPnhkbG8vifjMAAICPJHCxmZiYmMmTJ/PbnYjk5ORCQ0NjYmKEkWPChAlXrlzhPX727NmgQYPU1NTMzMzU1dW3bt3K5XKFsVEAAABpJbDguVzuu7O05ufnC2kGn19++SUrK4v3eNKkSX/99dfGjRvPnj07derU6dOnHzlyRBgbBQAAkFYC29rX1zciIsLU1HTAgAEsFovL5Z45cyYiImLEiBFCDVRaWnrixIlTp075+/sTkY+PT11dXWRkZEBAgFC3CwAAIE0E7sFv3LjRxsbGz8+vQ4cOVlZWHTp08PPz69GjR2RkpFAD3b9/n4g8PT35I25ubnfv3hXqRgEAAKSMwD14FRWVq1evXrx4MTk5+fHjxzo6Og4ODh4eHiwWS6iBDA0N5eTkCgsLzc3NeSPp6envX70NAAAAGmjuhDqLxfLy8vLy8hJNlMWLFx87dszMzMzY2Hj69Onnzp0jori4uM2bN/v4+IgmA7SI2tra3bt3Jycnm5mZjR07Vl1dnelEAACtjsCCLykpWbRoUWpqan19fcPxnj17bt26tcVzHDlyJDc3Nzc39+7du69fv7569Spv3M/Pz8rKatWqVS2+RRCS06dPjxgxor6+Xltb++DBg/PmzVuxYsXcuXOZzgUA0LoILPhJkybFx8cPGTJEWVm54Xjbtm2FkWPYsGENv+VwOLwHsbGxnp6euBVeUpSWlg4ZMiQoKOiXX37hjWzevHn69OnOzs4eHh7MZgMAaFUEFnxcXBzvEKso0/C1adOG98DLy4vD4VRVVTX6PQPE04oVK9hsNr/diSgsLCw6OnrBggV//fUXg8EAAFobgQWvp6dnIB5rnAcHBx88ePBD5rrhcDjNXG9fX1+PCXOE7e7du+9O1+zi4nLw4EFG8gAAtFoCC37kyJE7d+4MCwsTZZomeXp6fuDue0pKypQpUwQ9++bNm4qKipbLBU1QV1d/8OBBo8HCwkIcgAEAEDGBBW9qavr999+fPHnSzc2NzWY3fGr69OnCD/avkJCQkJCQD3mls7NzcnKyoGfbt2/f6AeBFvftt996eHgkJSU5ODjwRkpLS0+ePDlv3jxmgwEAtDYCC37OnDlE9PjxY/4MsnzCLvja2trS0lI1NTV5eXmhbghanJub29ChQ3v16jVw4EB3d/e7d+/+/vvvRkZGixcvZjoaAEDrInAmu0LBhBQlJSVlwoQJWlpaCgoKnTp1UlRU1NLSGjduXGpqqpC2CMJw+PDh6Ojo3NzclStXXrt2LTw8PDMzk+lQIIXOnz/fq1cvQ0PDvn37vntiCAA+buWYzMzMP/74Qxh78H/++ae3t7eBgUFwcLChoaGqqmpFRUVhYWFMTIyzs/PFixddXFxafKMgJAEBAVg7AITK1dU1Pj5eVVVVTU0tPj7eyMho6tSpmzZtYjoXgBgRWPBcLvfAgQO3b99uONFNfHx8TU2NMAp+wYIFvr6+hw4darRa3bJly0JCQubPn8+f+gYAWrlVq1bFx8f//vvvo0aN4o2EhIT88MMPISEhFhYWzGYDEB8CC3758uVLlizp3r17ZmammpqagYFBVlZWfX19bGysMHKkpqbu2bPn3bVoZWRkgoODBwwYIIyNAoAkioqKMjMz47c7Ee3YsWPfvn3Tpk27cOECg8EAxIrAc/C7d++eOXNmamrqzp07bW1tb9y4UVRUZGVlVV5eLowcJiYmCQkJTT4VHx9vamoqjI0CgCR6+fIlfzEqvg4dOjx8+JCRPADiSeAefHFxsbu7OxH17dt3xowZ9fX17du3X7BgwerVq4cMGdLiOcLCwsaPH19UVDRs2DBDQ0MVFZXKysqCgoKYmJhdu3bt3r27xbcIABJKSUkpLy+v0WB5ebmJiQkjeQDEk8CCV1dXz87OJiJdXV1ZWdmkpCQnJycNDY309HRh5Bg3bpy8vPy6deuio6MbjltZWe3bty8wMFAYGwUASTRmzJioqKi///7bycmJN7J27dpXr16tXr2a2WAAYkVgwfv5+UVGRnbq1CkoKMjOzi4qKmr9+vX79+/X1tYWUpTRo0cHBgYWFBQUFxeXlpaqq6tra2vr6ekJewV6AJAsGzZsOHXqlLOzs66urp6e3r17954/fz506FDcawPQkMCCj4yMrK6uPn78eFBQ0Jo1a1xcXA4cOCArK7tnzx7hpWGxWPr6+vr6+sLbBABIgdzc3M2bN2/YsCEjI0NbW/vo0aNYrhCgkeZuk9u9e7eMjAwR2draPnr0KCkpydDQUFNTU4TxAACaFhYWJg6LZQCILYFX0aupqT19+pT/LZvN7tevX3l5uZ6enkiCAQAAwKdrYg+ef5F8cHAwf112nnv37mloaIgiFwAAAHyGJgqev7Jnu3btlJSUGj7l5OQ0ZswYUeQCAACAz9BEwf/6669ElJubu3Pnzg4dOog8EgAAAHwugRfZJSYmijIHAAAAtKAmLrLj3R3Hn3O+oqJiwoQJJiYm3t7ef//9t2jjAQAAwKdoXPBlZWW9e/ceOnToxSIUhAAAGblJREFUmTNneCPBwcEHDhxwd3cvKSlxd3e/c+eOyEMCAADAx2l8iH7dunX3799PSkqyt7cnotzc3BMnTvz444+TJ0+uqalxd3dfvXr1b7/9xkRUAAAA+FCN9+BPnjw5Y8YMXrsT0dmzZ4lo6NChRCQvLz9y5MgbN26IOCIAAAB8rMYF/+DBA2tra/63Fy5csLW11dLS4n3bqVOnwsJC0aUDAACAT9K44NXU1Pgrvr9+/frKlSteXl78Z4uKijp27Ci6dAAAAPBJGhe8vb39jh07qquriWjfvn2VlZXe3t68p7hc7oEDBxru3wMAAIB4anyR3ZIlS5ydnW1tbW1sbE6ePGlqaurp6VlbW3vp0qWtW7cmJyf/+eefjAQFAACAD9d4D97W1jY+Pr5Lly7Jycl9+vSJjY2VlZUtLy/v37//9evXDx8+7ObmxkhQAAAA+HCN9+Dv37+vpqa2YcMG/kheXl5tbe2VK1d0dXVZLFZeXh6bzVZXVxdtTgAAAPgI/yn4wsJCHx+f2tra5t8zcODAH374QZipAAAA4LP8p+B1dXWzs7OZigIAAAAtRYaIWCzWL7/8curUqYcPHzKdBwAAAFqAHBFt2rTp3Llz27dvT0tLe/XqVY8ePaytrW1sbGxsbKysrBQVFZkOCQAAAB9Hjog8PDw8PDx435eXl9+5c+fmzZt///33rl27UlNTNTU1LSws7Ozs7Ozs7O3tO3fuzGhgAAAAeL/GV9Gz2Ww3Nzf+vXDV1dUZGRmpqampqalRUVEJCQmqqqpjx45dtWqVyKMCAADAh2pc8I0oKCj06NFDV1e3srIyISGBzWaPHDkyMDBQNOEAAADg07yn4J8/fz527Njr16/7+fktWbLEy8tLVlZWNMkAAADgk72n4BUVFW/cuHHq1ClMYAcAACBBGk9V24iysvLSpUuXLVsmmjQAAADQIt6zB09E48eP53A4IogCAAAALeU9e/BEJCcnN23aNBFEAQAAgJby/oIHAAAAiYOCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCckwHkADZ2dmXL18uLy+3sbHp37+/jAx+KwIAAHGHrmoOl8udP3++p6dnampqRUXF8uXLnZ2d8/Pzmc4FAADwHtiDb87evXvj4uIyMjLYbDZvZO3atV999dX169eZDQYAANA87ME3Z/fu3StWrOC3OxF99913RUVF2dnZDKYCAAB4LxR8c/Lz87t27dpwhMVimZmZPXz4kKlIAAAAHwIF3xxNTc2CgoJGg/n5+VpaWozkAQAA+EAo+OaMHDly5cqVtbW1/JEDBw7Iysp2796dwVQAAADvhYvsmjNlypTExMQePXqMGTNGXV390qVLf/75Z0xMDIvFYjoaAABAc7AH3xx5efmDBw9u3Ljx6dOnCQkJLi4umZmZNjY2TOcCAAB4D+zBv1+/fv369evHdAoAAICPgD14AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKYSCBwAAkEIoeAAAACmEggcAAJBCKHgAAAAphIIHAACQQih4AAAAKSSOBV9bW/v06dOamhqmgwAAAEgqMSr4lJSUCRMmaGlpKSgodOrUSVFRUUtLa9y4campqUxHAwAAkDDish78n3/+6e3tbWBgEBwcbGhoqKqqWlFRUVhYGBMT4+zsfPHiRRcXF6YzAgAASAxxKfgFCxb4+voeOnRITu4/kZYtWxYSEjJ//vyrV68ylQ0AAEDiiMsh+tTU1MDAwEbtTkQyMjLBwcEpKSmMpAIAAJBQ4lLwJiYmCQkJTT4VHx9vamoq4jwAAAASTVwO0YeFhY0fP76oqGjYsGGGhoYqKiqVlZUFBQUxMTG7du3avXs30wEBAAAkibgU/Lhx4+Tl5detWxcdHd1w3MrKat++fYGBgUwFAwAAkETiUvBENHr06MDAwIKCguLi4tLSUnV1dW1tbT09PRaLxXQ0AAAACSNGBU9ELBZLX19fX1+/4SCHw6mtrVVWVmYqFQAAgMQRr4JvUnBw8MGDB7lc7ntfmZaWtmLFCkHPvnnz5uXLly0aDQAAQExJQMF7enp+4O67rq5uQECAoGf//vtvGxublssFAAAgvlgfsmcsHXr16rVx40ZnZ2emgwAAQGthY2Ozf/9+a2tr0W9aXO6DbwiLzQAAAHwmMSp4LDYDAADQUsTlHDwWmwEAAGhB4lLwWGwGAACgBYnLIXosNgMAANCCxKXgsdgMAABACxKXQ/RYbAYAAKAFiUvBY7EZAACAFiQuBU8iWWzm5MmTaWlpvMcPHjxIS0tTU1NrqQ9vDbhcbnFxsY6ODtNBJMyzZ89UVVUVFBSYDiJJqqqqqqqq1NXVmQ4iYYqKivD/0I/17NkzNzc3TU1NYXz48+fPhfGxH6IVzWT3888/JyUl8b9NTk6+f/++lpYWg5EkTl1dXW5ubteuXZkOImEKCgrU1dXbtWvHdBBJUlFR8erVK21tbaaDSBIul5udnW1ubs50EAlTVFTUo0cPIf25tWnTZvny5SoqKsL48Oa1ooJvZPv27ampqdu2bWM6iCQpLy83MjIqKytjOoiE8fX1nTZtmo+PD9NBJMmePXuuXr2K628+yps3b1RVVTmc/2vv3qNizv8/gL8/TZdpKtOFXGq6SlGRLSQkTikKXa1LkRTb5pLcrW9Y0iJrndU2Zy210emgm0IsG2dbFaVYR6XlsIZcK5VUamZ+f3zOb37zq2b0baN85vn46/N+fWben/e8jno2n/n4TEtfL+Qz4+vrGxQU5OPj09cL6WX95Sp6AAAA6EUIeAAAAAZCwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyluwKuoqHT+8jqQT1lZWUVFpa9X8flB33oAP6E9oKSkhBsm9gBTf0IV90Y379+/b21t1dLS6uuFfGZqamr09PT6ehWfmbq6Oi6Xq6SkuH9P90BbW1tzc3Of3P/rs4af0B548+aNlpYWi8Xq64X0MsUNeAAAAAbDWwoAAAAGQsADAAAwEAIeAACAgRDwAAAADISABwAAYCAEPAAAAAMh4AEAABgIAQ8AAMBACHgAAAAGQsADAAAwEAIeAACAgRDwAAAADKSgAZ+YmGhvb8/lcqdOnfrnn3/29XL6o9bW1piYGCsrKw6HM2rUqH379r1//16yFw38oLa2tkmTJgUGBkoX0TdZiouLPTw8dHR0TE1N9+3bJ/0lWGhal9ra2vbs2WNhYaGhoTF27NjU1FTpvWhaZ5GRkevXr+9QlNMoJvRQrHiOHz9OCFm7dm16erq3tzebzS4rK+vrRfU7mzZtUlNTi4mJyc3N3b59u6qq6qpVq+hdaGB3rFu3jhCyaNEiSQV9k+XGjRsaGhoLFixIS0ujfwUfPHiQ3oWmybJx40Y1NbXvvvsuJycnIiKCEJKZmUnvQtM6EAqF58+f53A469atk67LaRQzeqhwAS8SiWxtbefPn08P29raLC0tQ0ND+3ZV/U17ezubzd6yZYuksnPnTmVl5Xfv3qGB3ZGVlaWtrW1kZCQJePRNDm9v7xkzZohEInoYFRXl4+MjRtNkE4lE+vr6UVFRksrEiRPnzp0rRtM6ycjI0NTUpN/QSge8nEYxpocKd4peIBDcuXPH39+fHiorK3t7e589e7ZvV9XfVFdXGxsbz5o1S1IxNTVtb29/+fIlGvhBDx8+DA4OPnbs2NChQyVF9E2W+vr67OzspUuXUhRFVw4cOJCRkUHQNLlYLBabzZYM2Ww2i8UiaFonLi4uhYWFd+7cMTQ0lK7LaRRjeqhwAV9dXU0IMTY2llRMTExevHjR3t7ed4vqd3g8XmVl5eTJk+nhu3fvjhw5Ym5ubmhoiAbK19raGhAQEBwc7OPjI11H32R5/PixSCRSUlLy8PDQ1tYePnz49u3bW1tbCZomG0VR69at+/HHH48dO3b9+vWdO3cWFRV99dVXBE3rREdHx8bGxsbGRk1NTboup1GM6aHCBXxtbS0hhMvlSipcLlcsFtN16Ky0tNTZ2bmsrCw5OZnFYqGB8kVFRbFYrL1793aoo2+y0L9Mw8PDHR0dU1JSQkJC4uLi1q5dS9A0ucLDw62srJYtW+bo6Lhjx47Fixe7uroSNK3b5DSKMT1UuIDX1dUlhDQ2NkoqDQ0NFEVpa2v33aL6qZqamqVLlzo4OPB4vDt37jg5ORE0UK709PTU1NRTp06pqqp22IW+yUK/r9qyZcuOHTs8PT23bt26bds2Pp/f1NSEpsnS2trq4OCgra199+7dpqam/Pz8vLy8L7/8UiwWo2ndJKdRjOmhwgX8sGHDCCECgUBSEQgE+vr6nX8jK7iqqipbW9vi4uKioqLMzEwTExO6jgbKkZ+fX1dXZ2JiQlEURVHXr19PSUmhKCorKwt9k8XAwIAQ4uDgIKmMHTtWLBYLBAI0TZbffvutoqLiyJEjo0aN4nA4kydP3rNnz+nTpx88eICmdZOcRjGmhwoX8Dwez9raOicnhx6KRKKcnBzpq8mAECISiXx8fEaPHl1cXDx+/HjpXWigHOHh4ZekWFlZTZ8+/dKlS05OTuibLKampjweLz8/X1IpLCxUU1MzMzND02RRUVEhhLx+/VpSefXqFSFET08PTesmOY1iTA+V+3oBnxpFUZs2bVqyZMnIkSMnTpyYlJRUVVWVnJzc1+vqXwoKCsrLy2fMmJGSkiJdX7hwIYfDQQNlsbS0tLS0lAy5XO7QoUPpT0YJIehbl5SVlTdu3BgVFSUUCh0dHYuKimJjY6Ojo+l3S2hal1xcXL744ouAgIBvvvnGyMiotLQ0JiYmJCRER0eHoGndIycLmBMTffYf9PpUYmKinZ2dlpbWlClTrl271tfL6Xf4fH6X/1qePXtGPwAN7I4JEyZI3+hGjL7JlpCQYGdnx+FwbGxsEhIShEKhZBea1qVXr16tXLnS1NSUzWbT95psaWmR7EXTOjM3N+9woxux3EYxoIeUWOqWkAAAAMAMCvcZPAAAgCJAwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAMAADAQAh4AAICBEPAAAAAMhIAHAABgIAQ8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAP0d5GRkVRXLl++3OvH4vP5AwcO7JWpBg4cyOfze2UqAOgB5b5eAAB8mKamZmJiYoeira3t27dvtbS0lixZkpSURAgJDAx8+/ZtVlZWh20AUEAIeIDPgJqamr+/f+d6c3Ozn5+fg4PDp18SAPRzOEUP8BlTV1cvKipqb28nhDg6OqakpJw5c4aiKCsrK8n2mzdvCCGJiYn29vYaGho2NjbSJwMaGxvDwsIMDQ0NDQ3Dw8NbWlo6H8XHx2fs2LGSoUgkMjAw+Prrrwkh9fX1y5cvHzZsmJqampmZ2a5du8RicecZ2Gz2iRMnJMPg4GAvLy/JUNbabt++7eHhoaOjo6en5+vrKxAI/k2vABQNAh6AIc6dO+fr6+vu7v7s2bP8/HzJ9oABAw4dOrR8+fKZM2eeOnXKxcUlJCQkISGBECIWiz09PU+dOhUVFXXw4MHKyspt27Z1ntnf3//WrVuPHz+mh0VFRdXV1YGBgYSQ1atXZ2VlRUZGnj592sfHJzo6+uTJk//VsmWt7d27d+7u7i0tLfHx8TExMQUFBaGhof+2RwCKBKfoAT4DNTU1FEVJV3bt2tUhjPX09NTV1YVC4ZAhQwghku23b9/u3LnzP//5T3R0NCHE09Pz/fv3e/bsCQ8P//333/Pz83Nycuj303PmzBkxYkRTU1OHo3t5eamqqubk5ERERBBC0tPTTUxMJk6cSAhpaGiIi4tbvHgx/fTLly+XlpbOnz+/m69LztrKy8tfvHiRkZHh5OREv7o//vijJ70DUFQIeIDPQOeL7Kytrbv53Lt379bV1bm6ur5+/ZquTJ8+/ciRIy0tLSUlJTo6Op6ennRdTU0tICCAvl5PGpfLdXNzy87OjoiIEIvFGRkZgYGB9B8cmZmZhBCxWPzkyZO8vLyKigo3N7fuvy45azMyMmKz2atXr960aZObm1tAQEBAQED3ZwYABDzAZ0DWRXbd8ejRI0LIpEmTOtQFAsHz58+HDRsmfW7AwMCgy0n8/PxWrFhRX19///79R48eLVq0iK6XlZVt3ry5tLS0ra1twoQJOjo63VmS5HN6OWuzsLDIzc3dunXrvHnzKIqaPHny1q1bPTw8ujM/ABB8Bg/AeIMGDSKECAQC8f9nYWFhYGDw7Nkz6cvinj9/3uUkc+fOFYvFFy9eTE9Pt7e3t7KyIoS8efPGyclJR0cnOzu7trb24sWLxsbG3VnSy5cvP7g2QoiLi0tBQcHz589PnDghFAq9vLzu37//75oBoEAQ8AAMZ2Njo6Kicu7cOUnlhx9+mDdvnlgsHjduXG1tbW5uLl1vb29PT0/vchJdXd1p06adOXMmPT1d8va9uLi4paUlLi5u4sSJSkpKjY2N9+7d6/LpSkpKNTU19HZ9fX1hYeEH15aWlmZra9vQ0DB48OCFCxcmJycLhUJZ8wNAZzhFD8AcKioqf//9d0lJiZ2dnfR2ZGTk+vXrGxsbR48eXVBQsHv37tjYWIqipk6d6uzsvHDhwt27d/N4PD6f39raKmtyf3//VatWtbe3S66hMzc3V1JS2rNnT2BgYG1tbUxMjFAoLCwsrKyspN/iS4wZM+bAgQMjRoxQV1fftWuXvr4+XdfX15e1Nmtr64qKiqCgoKCgIKFQePz4cS6XO378+I/XPQCmEQNA/7ZmzRo9PT1Zew0MDA4ePEhv5+XlmZiYaGpqvnnzRnpbKBTu379/1KhR6urqVlZW8fHxIpGIfkpDQ0NoaKiBgcGQIUPCwsLy8/NnzZrV5YFevHihpKTk5uYmXUxOTjY3N+dwOBMmTMjIyMjLy7O3t6fXo6enl5CQQD+svLzc2dmZw+HY2Njw+fyjR49u2LCB3iVnbWlpaWPGjFFXV9fV1Z0xY8aNGzd63kQAxUOJu7orBQAAAHzW8Bk8AAAAAyHgAQAAGAgBDwAAwEAIeAAAAAZCwAMAADAQAh4AAICBEPAAjDJw4EA+n9+DJ0ZFRR09erTX19NZRUWFk5OTUCj8BMcCUGQIeAAgZWVlFy5cWLJkySc41siRIy0tLX/66adPcCwARYaABwASGxsbFhamrPyJ7l0dERGxd+/etra2T3M4AMWEgAdgLBMTk+Tk5IiICCMjIzMzs/j4+KdPn3p6eg4cONDY2Dg1NZV+WHV1dVZW1oIFCwghmzdv5vF4khtctrS0cLncffv29eDoly5doijqypUr9DAtLU1FReX27duEEHt7e01NzaysrF54kQAgAwIegMmio6M9PT0fPnwYHBy8cuVKV1fX7du3P3/+3NnZOSQkpLm5mRBy8eJFCwuLIUOGEEL8/PyePHly8+ZN+ukXLlxoaGiQfLvMf8XNzS0oKCg8PLy1tbWhoWHNmjWbNm0aM2YMIYSiKBcXl/Pnz/feCwWAjhDwAEw2bdq0WbNmsVisZcuWEUICAwPHjx+vrKy8dOnSlpYWgUBACLly5cq4cePoxzs4OPB4vMzMTHp48uTJKVOmGBkZdZj21q1bPB7Pzs4uLy+Prhw6dKiurq7Dww4cOPD69eu9e/dGR0draWlt27ZNsmvcuHFXr17t/RcMAP8LAQ/AZJKvbR08eLD0kP7CVpFIRAj5559/hg4dStcpivL19aVPnjc1NWVnZ0u+/V1aXFzcypUrV6xYMWfOnBUrVhQWFsbGxrLZ7A4PGzRo0Pfffx8TExMfH//LL79IP2DYsGGPHz/GtfQAHw8CHoDJWCyWnCGttrZ2wIABkqGvr295eXlVVdW5c+fev3/v7+/f+SnW1tahoaHh4eE3btx4+vSph4dHeHi4urp650fOnz9fU1Nz+PDhTk5O0nUulysSierr63v4wgDgQz7RRbMA0G/p6uo2NjZKhpMmTdLX18/KyioqKvLw8NDT0+v8lC1bttAbo0aNOnv2rJzJ9+/fr6qq+uDBg6SkpJCQEEm9oaFBSUmJy+X20osAgI4Q8ACKztjYuLq6WjJksVje3t4pKSn37t1LSkr6NzNXVlZ+++23qamppaWlGzZsmD179qBBg+hd1dXVhoaGXZ5RAIBegVP0AIrOxcWlpKREuuLn5/fXX38pKyvPnj27x9OKRKKwsDBXV1cfH5+tW7cOGDBg/fr1kr0lJSUuLi49nhwAPggBD6CgbG1ttbS0CCHu7u5VVVUvXryQ7HJxcdHU1PTx8dHQ0Ojx/Hw+v6Sk5PDhwxRFcTicw4cPJycn0/8tXiwWX7161dPT89+/CgCQhZLc0QIAFJa/v//kyZMjIyPp4cOHD83MzM6fPz9z5syPcbjS0lIvL69Hjx6pqqp+jPkBgCDgAYAQcvPmzcWLF9O3mWtubo6KisrNzX348KGKisrHONyyZctGjx69Zs2ajzE5ANBwkR0AEHt7e3d3919//XXq1KkWFhYURWVnZ3+kdK+srCwvL//5558/xuQAIIF38ADwf9ra2q5du2ZlZUXfuRYAPl8IeAAAAAbCVfQAAAAM9D8deZGBo3YifAAAAABJRU5ErkJggg==" /><!-- --></p> -<p>Therefore, in Example 8 in <span class="citation">Massart et al. (1997)</span>, weighted regression is proposed which can be reproduced by</p> -<pre class="r"><code>attach(massart97ex3) -yx <- split(y, x) -ybar <- sapply(yx, mean) -s <- round(sapply(yx, sd), digits = 2) -w <- round(1 / (s^2), digits = 3) -weights <- w[factor(x)] -m <- lm(y ~ x, w = weights)</code></pre> +<p>Therefore, in Example 8 in <span class="citation">Massart et al. (1997)</span>, weighted regression is proposed which can be reproduced by the following code. Note that we are building the model on the mean values for each standard in order to be able to reproduce the results given in the book with the current version of chemCal.</p> +<pre class="r"><code>weights <- with(massart97ex3, { + yx <- split(y, x) + ybar <- sapply(yx, mean) + s <- round(sapply(yx, sd), digits = 2) + w <- round(1 / (s^2), digits = 3) +}) +massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean) + +m <- lm(y ~ x, w = weights, data = massart97ex3.means)</code></pre> <p>If we now want to predict a new x value from measured y values, we use the <code>inverse.predict</code> function:</p> <pre class="r"><code>inverse.predict(m, 15, ws=1.67)</code></pre> <pre><code>## $Prediction @@ -291,8 +297,8 @@ m <- lm(y ~ x, w = weights)</code></pre> ## [1] 36.20523 51.91526</code></pre> <p>The weight <code>ws</code> assigned to the measured y value has to be given by the user in the case of weighted regression, or alternatively, the approximate variance <code>var.s</code> at this location.</p> </div> -<div id="some-theory-for-inverse.predict" class="section level1"> -<h1>Some theory for <code>inverse.predict</code></h1> +<div id="background-for-inverse.predict" class="section level1"> +<h1>Background for <code>inverse.predict</code></h1> <p>Equation 8.28 in <span class="citation">Massart et al. (1997)</span> gives a general equation for predicting the standard error <span class="math inline">\(s_{\hat{x_s}}\)</span> for an <span class="math inline">\(x\)</span> value predicted from measurements of <span class="math inline">\(y\)</span> according to the linear calibration function <span class="math inline">\(y = b_0 + b_1 \cdot x\)</span>:</p> <span class="math display">\[\begin{equation} s_{\hat{x_s}} = \frac{s_e}{b_1} \sqrt{\frac{1}{w_s m} + \frac{1}{\sum{w_i}} + @@ -304,9 +310,11 @@ s_{\hat{x_s}} = \frac{s_e}{b_1} \sqrt{\frac{1}{w_s m} + \frac{1}{\sum{w_i}} + <span class="math display">\[\begin{equation} s_e = \sqrt{ \frac{\sum w_i (y_i - \hat{y_i})^2}{n - 2}} \end{equation}\]</span> -<p>where <span class="math inline">\(w_i\)</span> is the weight for calibration standard <span class="math inline">\(i\)</span>, <span class="math inline">\(y_i\)</span> is the mean <span class="math inline">\(y\)</span> value (!) observed for standard <span class="math inline">\(i\)</span>, <span class="math inline">\(\hat{y_i}\)</span> is the estimated value for standard <span class="math inline">\(i\)</span>, <span class="math inline">\(n\)</span> is the number calibration standards, <span class="math inline">\(w_s\)</span> is the weight attributed to the sample <span class="math inline">\(s\)</span>, <span class="math inline">\(m\)</span> is the number of replicate measurements of sample <span class="math inline">\(s\)</span>, <span class="math inline">\(\bar{y_s}\)</span> is the mean response for the sample, <span class="math inline">\(\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}\)</span> is the weighted mean of responses <span class="math inline">\(y_i\)</span>, and <span class="math inline">\(x_i\)</span> is the given <span class="math inline">\(x\)</span> value for standard <span class="math inline">\(i\)</span>.</p> +<p>In chemCal version before 0.2, I interpreted <span class="math inline">\(w_i\)</span> to be the weight for calibration standard <span class="math inline">\(i\)</span>, <span class="math inline">\(y_i\)</span> to be the mean value observed for standard <span class="math inline">\(i\)</span>, and <span class="math inline">\(n\)</span> to be the number of calibration standards. With this implementation I was able to reproduce the examples given in the book. However, as noted above, I was made aware of the fact that this way of calculation does not take the variation of the y values about the means into account. Furthermore, I noticed that for the case of unweighted linear calibration with replicate standards, <code>inverse.predict</code> produced different results than <code>calibrate</code> from the <code>investr</code> package when using the Wald method.</p> +<p>Both issues are now addressed in chemCal starting from version 0.2.1. Here, <span class="math inline">\(y_i\)</span> is calibration measurement <span class="math inline">\(i\)</span>, <span class="math inline">\(\hat{y_i}\)</span> is the estimated value for calibration measurement <span class="math inline">\(i\)</span> and <span class="math inline">\(n\)</span> is the total number of calibration measurements.</p> +<p><span class="math inline">\(w_s\)</span> is the weight attributed to the sample <span class="math inline">\(s\)</span>, <span class="math inline">\(m\)</span> is the number of replicate measurements of sample <span class="math inline">\(s\)</span>, <span class="math inline">\(\bar{y_s}\)</span> is the mean response for the sample, <span class="math inline">\(\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}\)</span> is the weighted mean of responses <span class="math inline">\(y_i\)</span>, and <span class="math inline">\(x_i\)</span> is the given <span class="math inline">\(x\)</span> value for standard <span class="math inline">\(i\)</span>.</p> <p>The weight <span class="math inline">\(w_s\)</span> for the sample should be estimated or calculated in accordance to the weights used in the linear regression.</p> -<p>I adjusted the above equation in order to be able to take a different precisions in standards and samples into account. In analogy to Equation 8.26 from we get</p> +<p>I had also adjusted the above equation in order to be able to take a different precisions in standards and samples into account. In analogy to Equation 8.26 from I am using</p> <span class="math display">\[\begin{equation} s_{\hat{x_s}} = \frac{1}{b_1} \sqrt{\frac{{s_s}^2}{w_s m} + {s_e}^2 \left( \frac{1}{\sum{w_i}} + @@ -315,6 +323,15 @@ s_{\hat{x_s}} = \frac{1}{b_1} \sqrt{\frac{{s_s}^2}{w_s m} + \end{equation}\]</span> <p>where I interpret <span class="math inline">\(\frac{{s_s}^2}{w_s}\)</span> as an estimator of the variance at location <span class="math inline">\(\hat{x_s}\)</span>, which can be replaced by a user-specified value using the argument <code>var.s</code> of the function <code>inverse.predict</code>.</p> <div id="refs" class="references"> +<div id="ref-amc89"> +<p>Analytical Methods Committee. 1989. “Robust Statistics — How Not to Reject Outliers. Part 1. Basic Concepts.” <em>The Analyst</em> 114: 1693–7.</p> +</div> +<div id="ref-castells00"> +<p>Castells, Reynaldo César, and Marcela Alejandra Castillo. 2000. “Systematic Errors: Detection and Correction by Means of Standard Calibration, Youden Calibration and Standard Additions Method in Conjunction with a Method Response Model.” <em>Analytica Chimica Acta</em> 423: 179–85.</p> +</div> +<div id="ref-currie97"> +<p>Currie, L. A. 1997. “Nomenclature in Evaluation of Analytical Methods Including Detection and Quantification Capabilities (IUPAC Recommendations 1995).” <em>Analytica Chimica Acta</em> 391: 105–26.</p> +</div> <div id="ref-massart97"> <p>Massart, D. L, B. G. M. Vandeginste, L. M. C. Buydens, S. De Jong, P. J. Lewi, and J Smeyers-Verbeke. 1997. <em>Handbook of Chemometrics and Qualimetrics: Part A</em>. Amsterdam: Elsevier.</p> </div> diff --git a/vignettes/figure/unnamed-chunk-1-1.pdf b/vignettes/figure/unnamed-chunk-1-1.pdf Binary files differdeleted file mode 100644 index c70b645..0000000 --- a/vignettes/figure/unnamed-chunk-1-1.pdf +++ /dev/null diff --git a/vignettes/figure/unnamed-chunk-2-1.pdf b/vignettes/figure/unnamed-chunk-2-1.pdf Binary files differdeleted file mode 100644 index c1b934d..0000000 --- a/vignettes/figure/unnamed-chunk-2-1.pdf +++ /dev/null diff --git a/vignettes/references.bib b/vignettes/references.bib new file mode 100644 index 0000000..a710662 --- /dev/null +++ b/vignettes/references.bib @@ -0,0 +1,35 @@ +% Encoding: UTF-8 +@book{massart97, + author = "Massart, D. L and Vandeginste, B. G. M. and Buydens, L. M. C. and De Jong, S. and Lewi, P. J. and Smeyers-Verbeke, J", + title = "Handbook of {C}hemometrics and {Q}ualimetrics: Part {A}", + publisher = "Elsevier", + address = "Amsterdam", + year = "1997" +} +@article{currie97, + author = "Currie, L. A.", + year = "1997", + title = "Nomenclature in evaluation of analytical methods including + detection and quantification capabilities ({IUPAC Recommendations 1995})", + journal = "Analytica Chimica Acta", + volume = "391", + pages = "105 - 126" +} +@article{amc89, + Title = {Robust statistics --- how not to reject outliers. {Part} 1. {Basic} concepts}, + Author = {{Analytical Methods Committee}}, + Journal = {The Analyst}, + Year = {1989}, + Pages = {1693--1697}, + Volume = {114} +} +@article{castells00, + author = {Reynaldo César Castells and Marcela Alejandra Castillo}, + title = {Systematic errors: detection and correction by means of standard calibration, Youden calibration and standard additions method in conjunction with a method response model}, + journal = "Analytica Chimica Acta", + volume = "423", + year = "2000", + pages = "179-185" +} + +@Comment{jabref-meta: databaseType:bibtex;} diff --git a/vignettes/refs.bib b/vignettes/refs.bib deleted file mode 100644 index 514d76b..0000000 --- a/vignettes/refs.bib +++ /dev/null @@ -1,7 +0,0 @@ -@book{massart97, - author = "Massart, D. L and Vandeginste, B. G. M. and Buydens, L. M. C. and De Jong, S. and Lewi, P. J. and Smeyers-Verbeke, J", - title = "Handbook of {C}hemometrics and {Q}ualimetrics: Part {A}", - publisher = "Elsevier", - address = "Amsterdam", - year = "1997" -} |
