Fix URLs in README, convert to roxygenv0.2.3

- The roxygen conversion was done using Rd2roxygen
- Also edit _pkgdown.yml to group the reference
- Use markdown bullet lists for lod and loq docs

5 files changed, 422 insertions, 4 deletions

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")
 {