aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--DESCRIPTION9
-rw-r--r--INDEX2
-rw-r--r--R/calplot.R2
-rw-r--r--R/inverse.predict.lm.R37
-rw-r--r--R/lod.R25
-rw-r--r--demo/00Index1
-rw-r--r--demo/massart97ex3.R15
-rw-r--r--demo/massart97ex8.R6
-rw-r--r--inst/doc/chemCal-001.pdf4
-rw-r--r--inst/doc/chemCal.log4
-rw-r--r--inst/doc/chemCal.pdfbin107614 -> 107514 bytes
-rw-r--r--inst/doc/chemCal.tex6
-rw-r--r--man/calplot.lm.Rd2
-rw-r--r--man/din32645.Rd12
-rw-r--r--man/inverse.predict.Rd13
-rw-r--r--man/lod.Rd58
-rw-r--r--man/massart97ex3.Rd26
17 files changed, 167 insertions, 55 deletions
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 <jranke@uni-bremen.de>
Maintainer: Johannes Ranke <jranke@uni-bremen.de>
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> </var/cache/fonts/pk/ljfour/jknappen/tc/tctt1000.6
xmf-tetex/fonts/type1/bluesky/cm/cmtt10.pfb></usr/share/texmf-tetex/fonts/type1
/bluesky/cm/cmr10.pfb></usr/share/texmf-tetex/fonts/type1/bluesky/cm/cmr12.pfb>
</usr/share/texmf-tetex/fonts/type1/bluesky/cm/cmr17.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
--- a/inst/doc/chemCal.pdf
+++ b/inst/doc/chemCal.pdf
Binary files 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

Contact - Imprint