From 69504b635d388507bce650c44b3bfe17cad3383e Mon Sep 17 00:00:00 2001 From: ranke Date: Fri, 12 May 2006 21:59:33 +0000 Subject: - Fixed the inverse prediction - Now I have a working approach for the calculation of LOD and LOQ, but it seems to be different from what everybody else is doing (e.g. Massart chaper 13). I like it, however. Maybe it even yields a paper. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@8 5fad18fb-23f0-0310-ab10-e59a3bee62b4 --- DESCRIPTION | 9 ++++---- INDEX | 2 ++ R/calplot.R | 2 +- R/inverse.predict.lm.R | 37 ++++++++++++++++-------------- R/lod.R | 25 ++++++++++++++++++++ demo/00Index | 1 - demo/massart97ex3.R | 15 ------------ demo/massart97ex8.R | 6 ++--- inst/doc/chemCal-001.pdf | 4 ++-- inst/doc/chemCal.log | 4 ++-- inst/doc/chemCal.pdf | Bin 107614 -> 107514 bytes inst/doc/chemCal.tex | 6 ++--- man/calplot.lm.Rd | 2 +- man/din32645.Rd | 12 +++++++++- man/inverse.predict.Rd | 13 +++++++---- man/lod.Rd | 58 +++++++++++++++++++++++++++++++++++++++++++++++ man/massart97ex3.Rd | 26 +++++++++++++++++++++ 17 files changed, 167 insertions(+), 55 deletions(-) create mode 100644 R/lod.R delete mode 100644 demo/massart97ex3.R create mode 100644 man/lod.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 963c663..9f45d05 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,15 +1,16 @@ Package: chemCal -Version: 0.05-7 -Date: 2006-05-11 +Version: 0.05-8 +Date: 2006-05-12 Title: Calibration functions for analytical chemistry Author: Johannes Ranke Maintainer: Johannes Ranke Depends: R Suggests: MASS Description: chemCal provides simple functions for plotting linear - calibration functions and for estimating standard errors for measurements + calibration functions and estimating standard errors for measurements according to the Handbook of Chemometrics and Qualimetrics: Part A - by Massart et al. + by Massart et al. There are also functions estimating the limit + of detection (LOQ) and limit of quantification (LOD). The functions work on model objects from - optionally weighted - linear regression (lm) or robust linear regression (rlm from the MASS package). License: GPL version 2 or newer diff --git a/INDEX b/INDEX index 6e85292..6ad7c6d 100644 --- a/INDEX +++ b/INDEX @@ -3,5 +3,7 @@ calplot Plot calibration graphs from univariate linear chemCal-package Calibration functions for analytical chemistry din32645 Calibration data from DIN 32645 inverse.predict Predict x from y for a linear calibration +lod Estimate a limit of detection (LOD) or + quantification (LOQ) massart97ex3 Calibration data from Massart et al. (1997), example 3 diff --git a/R/calplot.R b/R/calplot.R index cea1149..feb9727 100644 --- a/R/calplot.R +++ b/R/calplot.R @@ -27,7 +27,7 @@ calplot.lm <- function(object, xlim = "auto", ylim = "auto", pred.lim <- predict(m, newdata, interval = "prediction",level=level) conf.lim <- predict(m, newdata, interval = "confidence",level=level) if (xlim == "auto") xlim = c(0,max(x)) - if (ylim == "auto") ylim = range(c(pred.lim,y)) + if (ylim == "auto") ylim = range(c(pred.lim,y,0)) plot(1, type = "n", xlab = xlab, diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R index f09dc3e..b48c967 100644 --- a/R/inverse.predict.lm.R +++ b/R/inverse.predict.lm.R @@ -2,44 +2,47 @@ # and Qualimetrics, Part A, Massart et al (1997), page 200, validated with # Example 8 on the same page -inverse.predict <- function(object, newdata, - ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), - alpha=0.05, ss = "auto") +inverse.predict <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, ss = "auto") { UseMethod("inverse.predict") } -inverse.predict.default <- function(object, newdata, - ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), - alpha=0.05, ss = "auto") +inverse.predict.default <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, ss = "auto") { stop("Inverse prediction only implemented for univariate lm objects.") } -inverse.predict.lm <- function(object, newdata, - ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), - alpha=0.05, ss = "auto") +inverse.predict.lm <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, ss = "auto") { + if (ws == "auto") { + ws <- ifelse(length(object$weights) > 0, mean(object$weights), 1) + } if (length(object$weights) > 0) { - wx <- split(object$model$y,object$model$x) + wx <- split(object$weights,object$model$x) w <- sapply(wx,mean) } else { w <- rep(1,length(split(object$model$y,object$model$x))) } - .inverse.predict(object, newdata, ws, alpha, ss, w) + .inverse.predict(object = object, newdata = newdata, + ws = ws, alpha = alpha, ss = ss, w = w) } -inverse.predict.rlm <- function(object, newdata, - ws = mean(object$w), alpha=0.05, ss = "auto") +inverse.predict.rlm <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, ss = "auto") { + if (ws == "auto") { + ws <- mean(object$w) + } wx <- split(object$weights,object$model$x) w <- sapply(wx,mean) - .inverse.predict(object, newdata, ws, alpha, ss, w) + .inverse.predict(object = object, newdata = newdata, + ws = ws, alpha = alpha, ss = ss, w = w) } -.inverse.predict <- function(object, newdata, - ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), - alpha=0.05, ss = "auto", w) +.inverse.predict <- function(object, newdata, ws, alpha, ss, w) { if (length(object$coef) > 2) stop("More than one independent variable in your model - not implemented") diff --git a/R/lod.R b/R/lod.R new file mode 100644 index 0000000..9f90d48 --- /dev/null +++ b/R/lod.R @@ -0,0 +1,25 @@ +lod <- function(object, ..., alpha = 0.05, k = 1, n = 1, w = "auto") +{ + UseMethod("lod") +} + +loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w = "auto") +{ + lod(object = object, alpha = alpha, k = k, n = n, w = w) +} + +lod.default <- function(object, ..., alpha = 0.05, k = 1, n = 1, w = "auto") +{ + stop("lod is only implemented for univariate lm objects.") +} + +lod.lm <- function(object, ..., alpha = 0.05, k = 1, n = 1, w = "auto") +{ + f <- function(x) { + y <- predict(object, data.frame(x = x)) + p <- inverse.predict(object, rep(y, n), ws = w, alpha = alpha) + (p[["Prediction"]] - k * p[["Confidence"]])^2 + } + tmp <- optimize(f,interval=c(0,max(object$model$x))) + return(tmp$minimum) +} diff --git a/demo/00Index b/demo/00Index index 8749adf..a548abc 100644 --- a/demo/00Index +++ b/demo/00Index @@ -1,2 +1 @@ -massart97ex3 Analysis of example 3 in Massart (1997) massart97ex8 Analysis of example 8 in Massart (1997) diff --git a/demo/massart97ex3.R b/demo/massart97ex3.R deleted file mode 100644 index 731aba6..0000000 --- a/demo/massart97ex3.R +++ /dev/null @@ -1,15 +0,0 @@ -library(chemCal) -data(massart97ex3) -attach(massart97ex3) -yx <- split(y,factor(x)) -ybar <- sapply(yx,mean) -s <- round(sapply(yx,sd),digits=2) -w <- round(1/(si^2),digits=3) -data.frame(x=levels(factor(x)),ybar,s,w) - -weights <- w[factor(x)] -m <- lm(y ~ x,w=weights) -inverse.predict(m,15,ws=1.67) -inverse.predict(m,90,ws=0.145) - -calplot(m) diff --git a/demo/massart97ex8.R b/demo/massart97ex8.R index 332bd1d..dca065f 100644 --- a/demo/massart97ex8.R +++ b/demo/massart97ex8.R @@ -1,4 +1,3 @@ -library(chemCal) data(massart97ex3) attach(massart97ex3) xi <- levels(factor(x)) @@ -8,5 +7,6 @@ si <- round(sapply(yx,sd),digits=2) wi <- round(1/(si^2),digits=3) weights <- wi[factor(x)] m <- lm(y ~ x,w=weights) -inverse.predict(m,15) -inverse.predict(m,90) +inverse.predict(m, 15, ws = 1.67) +inverse.predict(m, 90, ws = 0.145) +calplot(m) diff --git a/inst/doc/chemCal-001.pdf b/inst/doc/chemCal-001.pdf index 6e2f9b8..ca1ecd6 100644 --- a/inst/doc/chemCal-001.pdf +++ b/inst/doc/chemCal-001.pdf @@ -2,8 +2,8 @@ %âãÏÓ\r 1 0 obj << -/CreationDate (D:20060511175039) -/ModDate (D:20060511175039) +/CreationDate (D:20060512230530) +/ModDate (D:20060512230530) /Title (R Graphics Output) /Producer (R 2.3.0) /Creator (R) diff --git a/inst/doc/chemCal.log b/inst/doc/chemCal.log index 805f382..ae10f7f 100644 --- a/inst/doc/chemCal.log +++ b/inst/doc/chemCal.log @@ -1,4 +1,4 @@ -This is pdfeTeX, Version 3.141592-1.21a-2.2 (Web2C 7.5.4) (format=pdflatex 2006.4.5) 11 MAY 2006 17:50 +This is pdfeTeX, Version 3.141592-1.21a-2.2 (Web2C 7.5.4) (format=pdflatex 2006.4.5) 12 MAY 2006 23:05 entering extended mode **chemCal.tex (./chemCal.tex @@ -346,4 +346,4 @@ type1/bluesky/cm/cmbx12.pfb> -Output written on chemCal.pdf (3 pages, 107614 bytes). +Output written on chemCal.pdf (3 pages, 107514 bytes). diff --git a/inst/doc/chemCal.pdf b/inst/doc/chemCal.pdf index e3c6a8c..9cec612 100644 Binary files a/inst/doc/chemCal.pdf and b/inst/doc/chemCal.pdf differ diff --git a/inst/doc/chemCal.tex b/inst/doc/chemCal.tex index a4b4459..4dbdce3 100644 --- a/inst/doc/chemCal.tex +++ b/inst/doc/chemCal.tex @@ -63,13 +63,13 @@ $Prediction [1] 5.865367 $`Standard Error` -[1] 10.66054 +[1] 0.892611 $Confidence -[1] 29.59840 +[1] 2.478285 $`Confidence Limits` -[1] -23.73303 35.46376 +[1] 3.387082 8.343652 \end{Soutput} \end{Schunk} diff --git a/man/calplot.lm.Rd b/man/calplot.lm.Rd index efdea3d..6d3f52d 100644 --- a/man/calplot.lm.Rd +++ b/man/calplot.lm.Rd @@ -30,7 +30,7 @@ The label of the y axis. } \item{alpha}{ - The confidence level for the confidence and prediction bands. + The error tolerance level for the confidence and prediction bands. } } \value{ diff --git a/man/din32645.Rd b/man/din32645.Rd index 1a3c046..d251b7c 100644 --- a/man/din32645.Rd +++ b/man/din32645.Rd @@ -13,9 +13,19 @@ data(din32645) m <- lm(y ~ x, data=din32645) calplot(m) -prediction <- inverse.predict(m,3500,alpha=0.01) +(prediction <- inverse.predict(m,3500,alpha=0.01)) # This should give 0.074 according to DIN (cited from the Dintest test data) round(prediction$Confidence,3) + +# According to Dintest, we should get 0.07, but we get 0.0759 +lod(m, alpha = 0.01) + +# In German, there is the "Erfassungsgrenze", with k = 2, +# and we should get 0.14 according to Dintest +lod(m, k = 2, alpha = 0.01) + +# According to Dintest, we should get 0.21, we get 0.212 +loq(m, alpha = 0.01) } \references{ DIN 32645 (equivalent to ISO 11843) diff --git a/man/inverse.predict.Rd b/man/inverse.predict.Rd index d773e58..8c2be9c 100644 --- a/man/inverse.predict.Rd +++ b/man/inverse.predict.Rd @@ -4,9 +4,8 @@ \alias{inverse.predict.rlm} \alias{inverse.predict.default} \title{Predict x from y for a linear calibration} -\usage{inverse.predict(object, newdata, - ws = ifelse(length(object$weights) > 0, mean(object$weights), 1), - alpha=0.05, ss = "auto") +\usage{inverse.predict(object, newdata, \dots, + ws, alpha=0.05, ss = "auto") } \arguments{ \item{object}{ @@ -17,12 +16,16 @@ \item{newdata}{ A vector of observed y values for one sample. } + \item{\dots}{ + Placeholder for further arguments that might be needed by + future implementations. + } \item{ws}{ The weight attributed to the sample. The default is to take the mean of the weights in the model, if there are any. } \item{alpha}{ - The confidence level for the confidence interval to be reported. + The error tolerance level for the confidence interval to be reported. } \item{ss}{ The estimated standard error of the sample measurements. The @@ -62,6 +65,6 @@ w <- round(1/(s^2),digits=3) weights <- w[factor(x)] m <- lm(y ~ x,w=weights) -inverse.predict(m,c(15)) +inverse.predict(m,15,ws = 1.67) } \keyword{manip} diff --git a/man/lod.Rd b/man/lod.Rd new file mode 100644 index 0000000..e6ce345 --- /dev/null +++ b/man/lod.Rd @@ -0,0 +1,58 @@ +\name{lod} +\alias{lod} +\alias{lod.lm} +\alias{lod.rlm} +\alias{lod.default} +\alias{loq} +\alias{loq.lm} +\alias{loq.rlm} +\alias{loq.default} +\title{Estimate a limit of detection (LOD) or quantification (LOQ)} +\usage{ + lod(object, \dots, alpha = 0.05, k = 1, n = 1, w = "auto") + loq(object, \dots, alpha = 0.05, k = 3, n = 1, w = "auto") +} +\arguments{ + \item{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. + } + \item{alpha}{ + The error tolerance for the prediction of x values in the calculation. + } + \item{\dots}{ + Placeholder for further arguments that might be needed by + future implementations. + } + \item{k}{ + The inverse of the maximum relative error tolerated at the + desired LOD/LOQ. + } + \item{n}{ + The number of replicate measurements for which the LOD/LOQ should be + specified. + } + \item{w}{ + The weight that should be attributed to the LOD/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. + } +} +\value{ + The estimated limit of detection for a model used for calibration. +} +\description{ + A useful operationalisation of a lower limit L of a measurement method is + simply the solution of the equation + \deqn{L = k c(L)}{L = k * c(L)} + where c(L) is half of the lenght of the confidence interval at the limit L. +} +\examples{ + data(din32645) + m <- lm(y ~ x, data = din32645) + lod(m) +} +\keyword{manip} diff --git a/man/massart97ex3.Rd b/man/massart97ex3.Rd index 1dabf90..2618709 100644 --- a/man/massart97ex3.Rd +++ b/man/massart97ex3.Rd @@ -10,6 +10,32 @@ 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) +# The following concords with the book +inverse.predict(m, 15, ws = 1.67) +inverse.predict(m, 90, ws = 0.145) + +calplot(m) + +m0 <- lm(y ~ x) +lod(m0) +lod(m) + +# Now we want to take advantage of the lower weights at lower y values +m2 <- lm(y ~ x, w = 1/y) +# To get a reasonable weight for the lod, we need to estimate it and predict +# a y value for it +yhat.lod <- predict(m,data.frame(x = lod(m2))) +lod(m2,w=1/yhat.lod,k=3) +} \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, p. 188 -- cgit v1.2.1