From 08465d77a6ca5a9656ac86047c6008f1e7f3e9c7 Mon Sep 17 00:00:00 2001 From: Johannes Ranke Date: Thu, 31 Mar 2022 19:21:03 +0200 Subject: Fix URLs in README, convert to roxygen - The roxygen conversion was done using Rd2roxygen - Also edit _pkgdown.yml to group the reference - Use markdown bullet lists for lod and loq docs --- R/calplot.R | 44 +++++++++++ R/chemCal-package.R | 194 +++++++++++++++++++++++++++++++++++++++++++++++++ R/inverse.predict.lm.R | 77 +++++++++++++++++++- R/lod.R | 61 ++++++++++++++++ R/loq.R | 50 +++++++++++++ 5 files changed, 422 insertions(+), 4 deletions(-) create mode 100644 R/chemCal-package.R (limited to 'R') diff --git a/R/calplot.R b/R/calplot.R index fd49a54..3c227d4 100644 --- a/R/calplot.R +++ b/R/calplot.R @@ -1,3 +1,45 @@ +#' Plot calibration graphs from univariate linear models +#' +#' Produce graphics of calibration data, the fitted model as well as +#' confidence, and, for unweighted regression, prediction bands. +#' +#' @aliases calplot calplot.default calplot.lm +#' @param object A univariate model object of class \code{\link{lm}} or +#' \code{\link[MASS:rlm]{rlm}} with model formula \code{y ~ x} or \code{y ~ x - +#' 1}. +#' @param xlim The limits of the plot on the x axis. +#' @param ylim The limits of the plot on the y axis. +#' @param xlab The label of the x axis. +#' @param ylab The label of the y axis. +#' @param legend_x An optional numeric value for adjusting the x coordinate of +#' the legend. +#' @param alpha 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{\link{lod}} (see note +#' below). +#' @param varfunc The variance function for generating the weights in the +#' model. Currently, this argument is ignored (see note below). +#' @return A plot of the calibration data, of your fitted model as well as +#' lines showing the confidence limits. Prediction limits are only shown for +#' models from unweighted regression. +#' @note Prediction bands for models from weighted linear regression require +#' weights for the data, for which responses should be predicted. Prediction +#' intervals using weights e.g. from a variance function are currently not +#' supported by the internally used function \code{\link{predict.lm}}, +#' therefore, \code{calplot} does not draw prediction bands for such models. +#' +#' It is possible to compare the \code{\link{calplot}} prediction bands with +#' the \code{\link{lod}} values if the \code{lod()} alpha and beta parameters +#' are half the value of the \code{calplot()} alpha parameter. +#' @author Johannes Ranke +#' @importFrom graphics legend matlines plot points +#' @examples +#' +#' data(massart97ex3) +#' m <- lm(y ~ x, data = massart97ex3) +#' calplot(m) +#' +#' @export calplot calplot <- function(object, xlim = c("auto", "auto"), ylim = c("auto", "auto"), xlab = "Concentration", ylab = "Response", @@ -7,6 +49,7 @@ calplot <- function(object, UseMethod("calplot") } +#' @export calplot.default <- function(object, xlim = c("auto","auto"), ylim = c("auto","auto"), xlab = "Concentration", ylab = "Response", @@ -16,6 +59,7 @@ calplot.default <- function(object, stop("Calibration plots only implemented for univariate lm objects.") } +#' @export calplot.lm <- function(object, xlim = c("auto","auto"), ylim = c("auto","auto"), xlab = "Concentration", ylab = "Response", diff --git a/R/chemCal-package.R b/R/chemCal-package.R new file mode 100644 index 0000000..8cc8c76 --- /dev/null +++ b/R/chemCal-package.R @@ -0,0 +1,194 @@ +#' Calibration data from DIN 32645 +#' +#' Sample dataset to test the package. +#' +#' +#' @name din32645 +#' @docType data +#' @format A dataframe containing 10 rows of x and y values. +#' @references DIN 32645 (equivalent to ISO 11843), Beuth Verlag, Berlin, 1994 +#' +#' Dintest. Plugin for MS Excel for evaluations of calibration data. Written by +#' Georg Schmitt, University of Heidelberg. Formerly available from the Website +#' of the University of Heidelberg. +#' +#' Currie, L. A. (1997) Nomenclature in evaluation of analytical methods +#' including detection and quantification capabilities (IUPAC Recommendations +#' 1995). Analytica Chimica Acta 391, 105 - 126. +#' @keywords datasets +#' @examples +#' +#' m <- lm(y ~ x, data = din32645) +#' calplot(m) +#' +#' ## Prediction of x with confidence interval +#' 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) +#' +#' ## Critical value: +#' 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") +#' round(crit$x, 2) +#' # and according to Dintest test data, we should get 0.0698 from +#' round(crit$x, 4) +#' +#' ## Limit of detection (smallest detectable value given alpha and beta) +#' # In German, the smallest detectable value is the "Erfassungsgrenze", and we +#' # should get 0.14 according to DIN, which we achieve by using the method +#' # described in it: +#' lod.din <- lod(m, alpha = 0.01, beta = 0.01, method = "din") +#' 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) +#' +#' # 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 +#' +NULL + + + + + +#' Calibration data from Massart et al. (1997), example 1 +#' +#' Sample dataset from p. 175 to test the package. +#' +#' +#' @name massart97ex1 +#' @docType data +#' @format A dataframe containing 6 observations of x and y data. +#' @source 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, Chapter 8. +#' @keywords datasets +NULL + + + + + +#' Calibration data from Massart et al. (1997), example 3 +#' +#' Sample dataset from p. 188 to test the package. +#' +#' +#' @name massart97ex3 +#' @docType data +#' @format A dataframe containing 6 levels of x values with 5 observations of y +#' for each level. +#' @source 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, Chapter 8. +#' @keywords datasets +#' @examples +#' +#' # 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(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, data = massart97ex3) +#' lod(m0) +#' +#' # Limit of quantification from unweighted regression +#' loq(m0) +#' +#' # For calculating the limit of quantification from a model from weighted +#' # regression, we need to supply weights, internally used for inverse.predict +#' # 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(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. +#' +NULL + + + + + +#' Cadmium concentrations measured by AAS as reported by Rocke and Lorenzato +#' (1995) +#' +#' Dataset reproduced from Table 1 in Rocke and Lorenzato (1995). +#' +#' +#' @name rl95_cadmium +#' @docType data +#' @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. +#' @keywords datasets +NULL + + + + + +#' Toluene amounts measured by GC/MS as reported by Rocke and Lorenzato (1995) +#' +#' Dataset reproduced from Table 4 in Rocke and Lorenzato (1995). The toluene +#' amount in the calibration samples is given in picograms per 100 µL. +#' Presumably this is the volume that was injected into the instrument. +#' +#' +#' @name rl95_toluene +#' @docType data +#' @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. +#' @keywords datasets +NULL + + + + + +#' Example data for calibration with replicates from University of Toronto +#' +#' Dataset read into R from +#' \url{https://sites.chem.utoronto.ca/chemistry/coursenotes/analsci/stats/files/example14.xls}. +#' +#' +#' @name utstats14 +#' @docType data +#' @format A tibble containing three replicate observations of the response for +#' five calibration concentrations. +#' @source David Stone and Jon Ellis (2011) Statistics in Analytical Chemistry. +#' Tutorial website maintained by the Departments of Chemistry, University of +#' Toronto. +#' \url{https://sites.chem.utoronto.ca/chemistry/coursenotes/analsci/stats/index.html} +#' @keywords datasets +NULL + + + diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R index 77d548f..031f9d4 100644 --- a/R/inverse.predict.lm.R +++ b/R/inverse.predict.lm.R @@ -1,19 +1,87 @@ -# This is an implementation of Equation (8.28) in the Handbook of Chemometrics -# and Qualimetrics, Part A, Massart et al (1997), page 200, validated with -# Example 8 on the same page, extended as specified in the package vignette - +#' Predict x from y for a linear calibration +#' +#' This function predicts x values using a univariate linear model that has +#' been generated for the purpose of calibrating a measurement method. +#' Prediction intervals are given at the specified confidence level. The +#' calculation method was taken from Massart et al. (1997). In particular, +#' Equations 8.26 and 8.28 were combined in order to yield a general treatment +#' of inverse prediction for univariate linear models, taking into account +#' weights that have been used to create the linear model, and at the same time +#' providing the possibility to specify a precision in sample measurements +#' differing from the precision in standard samples used for the calibration. +#' This is elaborated in the package vignette. +#' +#' This is an implementation of Equation (8.28) in the Handbook of Chemometrics +#' and Qualimetrics, Part A, Massart et al (1997), page 200, validated with +#' Example 8 on the same page, extended as specified in the package vignette +#' +#' @aliases inverse.predict inverse.predict.lm inverse.predict.rlm +#' inverse.predict.default +#' @param object A univariate model object of class \code{\link{lm}} or +#' \code{\link[MASS:rlm]{rlm}} with model formula \code{y ~ x} or \code{y ~ x - +#' 1}. +#' @param newdata A vector of observed y values for one sample. +#' @param \dots Placeholder for further arguments that might be needed by +#' future implementations. +#' @param ws The weight attributed to the sample. This argument is obligatory +#' if \code{object} has weights. +#' @param alpha The error tolerance level for the confidence interval to be +#' reported. +#' @param var.s The estimated variance of the sample measurements. The default +#' is to take the residual standard error from the calibration and to adjust it +#' using \code{ws}, if applicable. This means that \code{var.s} overrides +#' \code{ws}. +#' @return A list containing the predicted x value, its standard error and a +#' confidence interval. +#' @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, p. 200 +#' @importFrom stats optimize predict qt +#' @examples +#' +#' # This is example 7 from Chapter 8 in Massart et al. (1997) +#' 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 +#' +#' +#' @export inverse.predict <- function(object, newdata, ..., ws = "auto", alpha = 0.05, var.s = "auto") { UseMethod("inverse.predict") } +#' @export inverse.predict.default <- function(object, newdata, ..., ws = "auto", alpha = 0.05, var.s = "auto") { stop("Inverse prediction only implemented for univariate lm objects.") } +#' @export inverse.predict.lm <- function(object, newdata, ..., ws = "auto", alpha = 0.05, var.s = "auto") { @@ -31,6 +99,7 @@ inverse.predict.lm <- function(object, newdata, ..., ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) } +#' @export inverse.predict.rlm <- function(object, newdata, ..., ws = "auto", alpha = 0.05, var.s = "auto") { diff --git a/R/lod.R b/R/lod.R index bccdf3e..bf6e928 100644 --- a/R/lod.R +++ b/R/lod.R @@ -1,13 +1,74 @@ +#' Estimate a limit of detection (LOD) +#' +#' The decision limit (German: Nachweisgrenze) is defined as the signal or +#' analyte concentration that is significantly different from the blank signal +#' with a first order error alpha (one-sided significance test). The detection +#' limit, or more precise, the minimum detectable value (German: +#' Erfassungsgrenze), is then defined as the signal or analyte concentration +#' where the probability that the signal is not detected although the analyte +#' is present (type II or false negative error), is beta (also a one-sided +#' significance test). +#' +#' @aliases lod lod.lm lod.rlm lod.default +#' @param object A univariate model object of class \code{\link{lm}} or +#' \code{\link[MASS:rlm]{rlm}} with model formula \code{y ~ x} or \code{y ~ x - +#' 1}, optionally from a weighted regression. +#' @param \dots Placeholder for further arguments that might be needed by +#' future implementations. +#' @param alpha The error tolerance for the decision limit (critical value). +#' @param beta The error tolerance beta for the detection limit. +#' @param method The \dQuote{default} method uses a prediction interval at the +#' LOD for the estimation of the LOD, which obviously requires iteration. This +#' is described for example in Massart, p. 432 ff. The \dQuote{din} method +#' uses the prediction interval at x = 0 as an approximation. +#' @param tol When the \dQuote{default} method is used, the default tolerance +#' for the LOD on the x scale is the value of the smallest non-zero standard +#' divided by 1000. Can be set to a numeric value to override this. +#' @return A list containig the corresponding x and y values of the estimated +#' limit of detection of a model used for calibration. +#' @note +#' * The default values for alpha and beta are the ones recommended by IUPAC. +#' * The estimation of the LOD in terms of the analyte amount/concentration xD +#' from the LOD in the signal domain SD is done by simply inverting the +#' calibration function (i.e. assuming a known calibration function). +#' * The calculation of a LOD from weighted calibration models requires a +#' weights argument for the internally used \code{\link{predict.lm}} +#' function, which is currently not supported in R. +#' @seealso Examples for \code{\link{din32645}} +#' @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, Chapter 13.7.8 +#' +#' J. Inczedy, T. Lengyel, and A.M. Ure (2002) International Union of Pure and +#' Applied Chemistry Compendium of Analytical Nomenclature: Definitive Rules. +#' Web edition. +#' +#' Currie, L. A. (1997) Nomenclature in evaluation of analytical methods +#' including detection and quantification capabilities (IUPAC Recommendations +#' 1995). Analytica Chimica Acta 391, 105 - 126. +#' @importFrom stats optimize predict +#' @examples +#' +#' m <- lm(y ~ x, data = din32645) +#' lod(m) +#' +#' # The critical value (decision limit, German Nachweisgrenze) can be obtained +#' # by using beta = 0.5: +#' lod(m, alpha = 0.01, beta = 0.5) +#' +#' @export lod <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") { UseMethod("lod") } +#' @export lod.default <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") { stop("lod is only implemented for univariate lm objects.") } +#' @export lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") { if (length(object$weights) > 0) { diff --git a/R/loq.R b/R/loq.R index 26d6a53..a4dbcfe 100644 --- a/R/loq.R +++ b/R/loq.R @@ -1,15 +1,65 @@ +#' Estimate a limit of quantification (LOQ) +#' +#' The limit of quantification is the x value, where the relative error of the +#' quantification given the calibration model reaches a prespecified value 1/k. +#' Thus, it is the solution of the equation \deqn{L = k c(L)}{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{\link{inverse.predict}}, and L is obtained by iteration. +#' +#' @aliases loq loq.lm loq.rlm loq.default +#' @param object A univariate model object of class \code{\link{lm}} or +#' \code{\link[MASS:rlm]{rlm}} with model formula \code{y ~ x} or \code{y ~ x - +#' 1}, optionally from a weighted regression. If weights are specified in the +#' model, either \code{w.loq} or \code{var.loq} have to be specified. +#' @param alpha The error tolerance for the prediction of x values in the +#' calculation. +#' @param \dots Placeholder for further arguments that might be needed by +#' future implementations. +#' @param k The inverse of the maximum relative error tolerated at the desired +#' LOQ. +#' @param n The number of replicate measurements for which the LOQ should be +#' specified. +#' @param w.loq 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{\link{massart97ex3}} for an example how to take +#' advantage of knowledge about the variance function. +#' @param var.loq The approximate variance at the LOQ. The default value is +#' calculated from the model. +#' @param tol The default tolerance for the LOQ on the x scale is the value of +#' the smallest non-zero standard divided by 1000. Can be set to a numeric +#' value to override this. +#' @return The estimated limit of quantification for a model used for +#' calibration. +#' @note +#' * IUPAC recommends to base the LOQ on the standard deviation of the +#' signal where x = 0. +#' * The calculation of a LOQ based on weighted regression is non-standard and +#' therefore not tested. Feedback is welcome. +#' @seealso Examples for \code{\link{din32645}} +#' @examples +#' +#' m <- lm(y ~ x, data = massart97ex1) +#' loq(m) +#' +#' # We can get better by using replicate measurements +#' loq(m, n = 3) +#' +#' @export loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", var.loq = "auto", tol = "default") { UseMethod("loq") } +#' @export loq.default <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", var.loq = "auto", tol = "default") { stop("loq is only implemented for univariate lm objects.") } +#' @export loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", var.loq = "auto", tol = "default") { -- cgit v1.2.1