From d8d6012e98fb4c7f158bcc7c173407c2b5f3e42e Mon Sep 17 00:00:00 2001 From: ranke Date: Sat, 22 Aug 2015 09:03:10 +0000 Subject: Get rid of the branched svn layout I never used git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@36 5fad18fb-23f0-0310-ab10-e59a3bee62b4 --- .Rbuildignore | 9 ++ ChangeLog | 11 ++ DESCRIPTION | 18 +++ GNUmakefile | 57 +++++++++ NAMESPACE | 6 + R/calplot.R | 80 ++++++++++++ R/inverse.predict.lm.R | 97 ++++++++++++++ R/lod.R | 55 ++++++++ R/loq.R | 41 ++++++ TODO | 5 + branches/0.1/chemCal/.Rbuildignore | 9 -- branches/0.1/chemCal/ChangeLog | 11 -- branches/0.1/chemCal/DESCRIPTION | 18 --- branches/0.1/chemCal/GNUmakefile | 57 --------- branches/0.1/chemCal/NAMESPACE | 6 - branches/0.1/chemCal/R/calplot.R | 80 ------------ branches/0.1/chemCal/R/inverse.predict.lm.R | 97 -------------- branches/0.1/chemCal/R/lod.R | 55 -------- branches/0.1/chemCal/R/loq.R | 41 ------ branches/0.1/chemCal/TODO | 5 - branches/0.1/chemCal/data/din32645.rda | Bin 453 -> 0 bytes branches/0.1/chemCal/data/massart97ex1.rda | Bin 193 -> 0 bytes branches/0.1/chemCal/data/massart97ex3.rda | Bin 296 -> 0 bytes branches/0.1/chemCal/man/calplot.lm.Rd | 69 ---------- branches/0.1/chemCal/man/chemCal-package.Rd | 16 --- branches/0.1/chemCal/man/din32645.Rd | 61 --------- branches/0.1/chemCal/man/inverse.predict.Rd | 69 ---------- branches/0.1/chemCal/man/lod.Rd | 88 ------------- branches/0.1/chemCal/man/loq.Rd | 82 ------------ branches/0.1/chemCal/man/massart97ex1.Rd | 17 --- branches/0.1/chemCal/man/massart97ex3.Rd | 51 -------- branches/0.1/chemCal/tests/din32645.R | 7 - branches/0.1/chemCal/tests/din32645.Rout.save | 48 ------- branches/0.1/chemCal/tests/massart97.R | 25 ---- branches/0.1/chemCal/tests/massart97.Rout.save | 128 ------------------- branches/0.1/chemCal/vignettes/chemCal-001.pdf | Bin 12038 -> 0 bytes branches/0.1/chemCal/vignettes/chemCal-002.pdf | Bin 5031 -> 0 bytes branches/0.1/chemCal/vignettes/chemCal.Rnw | 130 ------------------- branches/0.1/chemCal/vignettes/chemCal.pdf | Bin 179149 -> 0 bytes branches/0.1/chemCal/vignettes/chemCal.tex | 169 ------------------------- data/din32645.rda | Bin 0 -> 453 bytes data/massart97ex1.rda | Bin 0 -> 193 bytes data/massart97ex3.rda | Bin 0 -> 296 bytes man/calplot.lm.Rd | 69 ++++++++++ man/chemCal-package.Rd | 16 +++ man/din32645.Rd | 61 +++++++++ man/inverse.predict.Rd | 69 ++++++++++ man/lod.Rd | 88 +++++++++++++ man/loq.Rd | 82 ++++++++++++ man/massart97ex1.Rd | 17 +++ man/massart97ex3.Rd | 51 ++++++++ tests/din32645.R | 7 + tests/din32645.Rout.save | 48 +++++++ tests/massart97.R | 25 ++++ tests/massart97.Rout.save | 128 +++++++++++++++++++ vignettes/chemCal-001.pdf | Bin 0 -> 12038 bytes vignettes/chemCal-002.pdf | Bin 0 -> 5031 bytes vignettes/chemCal.Rnw | 130 +++++++++++++++++++ vignettes/chemCal.pdf | Bin 0 -> 179149 bytes vignettes/chemCal.tex | 169 +++++++++++++++++++++++++ 60 files changed, 1339 insertions(+), 1339 deletions(-) create mode 100644 .Rbuildignore create mode 100644 ChangeLog create mode 100644 DESCRIPTION create mode 100644 GNUmakefile create mode 100644 NAMESPACE create mode 100644 R/calplot.R create mode 100644 R/inverse.predict.lm.R create mode 100644 R/lod.R create mode 100644 R/loq.R create mode 100644 TODO delete mode 100644 branches/0.1/chemCal/.Rbuildignore delete mode 100644 branches/0.1/chemCal/ChangeLog delete mode 100644 branches/0.1/chemCal/DESCRIPTION delete mode 100644 branches/0.1/chemCal/GNUmakefile delete mode 100644 branches/0.1/chemCal/NAMESPACE delete mode 100644 branches/0.1/chemCal/R/calplot.R delete mode 100644 branches/0.1/chemCal/R/inverse.predict.lm.R delete mode 100644 branches/0.1/chemCal/R/lod.R delete mode 100644 branches/0.1/chemCal/R/loq.R delete mode 100644 branches/0.1/chemCal/TODO delete mode 100644 branches/0.1/chemCal/data/din32645.rda delete mode 100644 branches/0.1/chemCal/data/massart97ex1.rda delete mode 100644 branches/0.1/chemCal/data/massart97ex3.rda delete mode 100644 branches/0.1/chemCal/man/calplot.lm.Rd delete mode 100644 branches/0.1/chemCal/man/chemCal-package.Rd delete mode 100644 branches/0.1/chemCal/man/din32645.Rd delete mode 100644 branches/0.1/chemCal/man/inverse.predict.Rd delete mode 100644 branches/0.1/chemCal/man/lod.Rd delete mode 100644 branches/0.1/chemCal/man/loq.Rd delete mode 100644 branches/0.1/chemCal/man/massart97ex1.Rd delete mode 100644 branches/0.1/chemCal/man/massart97ex3.Rd delete mode 100644 branches/0.1/chemCal/tests/din32645.R delete mode 100644 branches/0.1/chemCal/tests/din32645.Rout.save delete mode 100644 branches/0.1/chemCal/tests/massart97.R delete mode 100644 branches/0.1/chemCal/tests/massart97.Rout.save delete mode 100644 branches/0.1/chemCal/vignettes/chemCal-001.pdf delete mode 100644 branches/0.1/chemCal/vignettes/chemCal-002.pdf delete mode 100644 branches/0.1/chemCal/vignettes/chemCal.Rnw delete mode 100644 branches/0.1/chemCal/vignettes/chemCal.pdf delete mode 100644 branches/0.1/chemCal/vignettes/chemCal.tex create mode 100644 data/din32645.rda create mode 100644 data/massart97ex1.rda create mode 100644 data/massart97ex3.rda create mode 100644 man/calplot.lm.Rd create mode 100644 man/chemCal-package.Rd create mode 100644 man/din32645.Rd create mode 100644 man/inverse.predict.Rd create mode 100644 man/lod.Rd create mode 100644 man/loq.Rd create mode 100644 man/massart97ex1.Rd create mode 100644 man/massart97ex3.Rd create mode 100644 tests/din32645.R create mode 100644 tests/din32645.Rout.save create mode 100644 tests/massart97.R create mode 100644 tests/massart97.Rout.save create mode 100644 vignettes/chemCal-001.pdf create mode 100644 vignettes/chemCal-002.pdf create mode 100644 vignettes/chemCal.Rnw create mode 100644 vignettes/chemCal.pdf create mode 100644 vignettes/chemCal.tex diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..0b413f5 --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,9 @@ +GNUmakefile +out$ +toc$ +bbl$ +blg$ +aux$ +log$ +vignettes/chemCal.tex +vignettes/chemCal.pdf diff --git a/ChangeLog b/ChangeLog new file mode 100644 index 0000000..c258c26 --- /dev/null +++ b/ChangeLog @@ -0,0 +1,11 @@ +2014-04-24 Johannes Ranke for chemCal (0.1-33) + + * 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. + +Changes performed in earlier versions are documented in the subversion log +files found at http://kriemhild.uft.uni-bremen.de/viewcvs/?root=chemCal diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..2608804 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,18 @@ +Package: chemCal +Version: 0.1-35 +Date: 2015-07-02 +Title: Calibration functions for analytical chemistry +Author: Johannes Ranke +Maintainer: Johannes Ranke +Suggests: MASS +Description: chemCal provides simple functions for plotting linear + calibration functions and estimating standard errors for measurements + according to the Handbook of Chemometrics and Qualimetrics: Part A + by Massart et al. There are also functions estimating the limit + of detection (LOD) and limit of quantification (LOQ). + The functions work on model objects from - optionally weighted - linear + regression (lm) or robust linear regression (rlm from the MASS package). +License: GPL (>= 2) +URL: http://www.r-project.org, + http://www.uft.uni-bremen.de/chemie/ranke, + http://kriemhild.uft.uni-bremen.de/viewcvs/?root=chemCal diff --git a/GNUmakefile b/GNUmakefile new file mode 100644 index 0000000..803c503 --- /dev/null +++ b/GNUmakefile @@ -0,0 +1,57 @@ +PKGNAME := $(shell sed -n "s/Package: *\([^ ]*\)/\1/p" DESCRIPTION) +PKGVERS := $(shell sed -n "s/Version: *\([^ ]*\)/\1/p" DESCRIPTION) +PKGSRC := $(shell basename $(PWD)) + +# Specify the directory holding R binaries. To use an alternate R build (say a +# pre-prelease version) use `make RBIN=/path/to/other/R/` or `export RBIN=...` +# 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`") + +.PHONY: help + +help: + @echo "\nExecute development tasks for $(PKGNAME)\n" + @echo "Usage: \`make \` where is one of:" + @echo "" + @echo "Development Tasks" + @echo "-----------------" + @echo " build Create the package" + @echo " build-no-vignettes Create the package without rebuilding vignettes" + @echo " check Invoke build and then check the package" + @echo " check-no-vignettes Invoke build without rebuilding vignettes, and then check" + @echo " install Invoke build and then install the result" + @echo " install-no-vignettes Invoke build without rebuilding vignettes and then install the result" + @echo "" + @echo "Using R in: $(RBIN)" + @echo "Set the RBIN environment variable to change this." + @echo "" + + +#------------------------------------------------------------------------------ +# Development Tasks +#------------------------------------------------------------------------------ + +build: + cd ..;\ + "$(RBIN)/R" CMD build $(PKGSRC) + +build-no-vignettes: + cd ..;\ + "$(RBIN)/R" CMD build $(PKGSRC) --no-build-vignettes + +install: build + cd ..;\ + "$(RBIN)/R" CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz + +install-no-vignettes: build-no-vignettes + cd ..;\ + "$(RBIN)/R" CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz + +check: build + cd ..;\ + "$(RBIN)/R" CMD check --as-cran $(PKGNAME)_$(PKGVERS).tar.gz + +check-no-vignettes: build-no-vignettes + cd ..;\ + "$(RBIN)/R" CMD check --as-cran $(PKGNAME)_$(PKGVERS).tar.gz diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..6314e29 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,6 @@ +# Export all names +exportPattern(".") +S3method(lod, lm) +S3method(loq, lm) +S3method(inverse.predict, lm) +S3method(inverse.predict, rlm) diff --git a/R/calplot.R b/R/calplot.R new file mode 100644 index 0000000..6aed9c0 --- /dev/null +++ b/R/calplot.R @@ -0,0 +1,80 @@ +calplot <- function(object, + xlim = c("auto", "auto"), ylim = c("auto", "auto"), + xlab = "Concentration", ylab = "Response", alpha = 0.05, + varfunc = NULL) +{ + UseMethod("calplot") +} + +calplot.default <- function(object, + xlim = c("auto","auto"), ylim = c("auto","auto"), + xlab = "Concentration", ylab = "Response", + alpha=0.05, varfunc = NULL) +{ + stop("Calibration plots only implemented for univariate lm objects.") +} + +calplot.lm <- function(object, + xlim = c("auto","auto"), ylim = c("auto","auto"), + xlab = "Concentration", ylab = "Response", alpha=0.05, + varfunc = NULL) +{ + if (length(object$coef) > 2) + stop("More than one independent variable in your model - not implemented") + + if (alpha <= 0 | alpha >= 1) + stop("Alpha should be between 0 and 1 (exclusive)") + + m <- object + level <- 1 - alpha + y <- m$model[[1]] + x <- m$model[[2]] + if (xlim[1] == "auto") xlim[1] <- 0 + if (xlim[2] == "auto") xlim[2] <- max(x) + xlim <- as.numeric(xlim) + newdata <- list( + x = seq(from = xlim[[1]], to = xlim[[2]], length=250)) + names(newdata) <- names(m$model)[[2]] + if (is.null(varfunc)) { + varfunc <- if (length(m$weights)) { + function(variable) mean(m$weights) + } else function(variable) rep(1,250) + } + pred.lim <- predict(m, newdata, interval = "prediction", + level=level, weights.newdata = varfunc(m)) + conf.lim <- predict(m, newdata, interval = "confidence", + level=level) + yrange.auto <- range(c(0,pred.lim)) + if (ylim[1] == "auto") ylim[1] <- yrange.auto[1] + if (ylim[2] == "auto") ylim[2] <- yrange.auto[2] + plot(1, + type = "n", + xlab = xlab, + ylab = ylab, + xlim = as.numeric(xlim), + ylim = as.numeric(ylim) + ) + points(x,y, pch = 21, bg = "yellow") + matlines(newdata[[1]], pred.lim, lty = c(1, 4, 4), + col = c("black", "red", "red")) + if (length(object$weights) > 0) { + legend(min(x), + max(pred.lim, na.rm = TRUE), + legend = c("Fitted Line", "Confidence Bands"), + lty = c(1, 3), + lwd = 2, + col = c("black", "green4"), + horiz = FALSE, cex = 0.9, bg = "gray95") + } else { + matlines(newdata[[1]], conf.lim, lty = c(1, 3, 3), + col = c("black", "green4", "green4")) + legend(min(x), + max(pred.lim, na.rm = TRUE), + legend = c("Fitted Line", "Confidence Bands", + "Prediction Bands"), + lty = c(1, 3, 4), + lwd = 2, + col = c("black", "green4", "red"), + horiz = FALSE, cex = 0.9, bg = "gray95") + } +} diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R new file mode 100644 index 0000000..927e672 --- /dev/null +++ b/R/inverse.predict.lm.R @@ -0,0 +1,97 @@ +# 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 + +inverse.predict <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, var.s = "auto") +{ + UseMethod("inverse.predict") +} + +inverse.predict.default <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, var.s = "auto") +{ + stop("Inverse prediction only implemented for univariate lm objects.") +} + +inverse.predict.lm <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, var.s = "auto") +{ + yname = names(object$model)[[1]] + xname = names(object$model)[[2]] + if (ws == "auto") { + ws <- ifelse(length(object$weights) > 0, mean(object$weights), 1) + } + if (length(object$weights) > 0) { + wx <- split(object$weights,object$model[[xname]]) + w <- sapply(wx,mean) + } else { + w <- rep(1,length(split(object$model[[yname]],object$model[[xname]]))) + } + .inverse.predict(object = object, newdata = newdata, + ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) +} + +inverse.predict.rlm <- function(object, newdata, ..., + ws = "auto", alpha = 0.05, var.s = "auto") +{ + yname = names(object$model)[[1]] + xname = names(object$model)[[2]] + if (ws == "auto") { + ws <- mean(object$w) + } + wx <- split(object$weights,object$model[[xname]]) + w <- sapply(wx,mean) + .inverse.predict(object = object, newdata = newdata, + ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) +} + +.inverse.predict <- function(object, newdata, ws, alpha, var.s, w, xname, yname) +{ + if (length(object$coef) > 2) + stop("More than one independent variable in your model - not implemented") + + if (alpha <= 0 | alpha >= 1) + stop("Alpha should be between 0 and 1 (exclusive)") + + ybars <- mean(newdata) + m <- length(newdata) + + yx <- split(object$model[[yname]], object$model[[xname]]) + n <- length(yx) + df <- n - length(object$coef) + x <- as.numeric(names(yx)) + ybar <- sapply(yx,mean) + yhatx <- split(object$fitted.values, object$model[[xname]]) + yhat <- sapply(yhatx, mean) + se <- sqrt(sum(w * (ybar - yhat)^2)/df) + + if (var.s == "auto") { + var.s <- se^2/ws + } + + b1 <- object$coef[[xname]] + + ybarw <- sum(w * ybar)/sum(w) + +# This is the adapted form of equation 8.28 (see package vignette) + sxhats <- 1/b1 * sqrt( + (var.s / m) + + se^2 * (1/sum(w) + + (ybars - ybarw)^2 * sum(w) / + (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2))) + ) + + if (names(object$coef)[1] == "(Intercept)") { + b0 <- object$coef[["(Intercept)"]] + } else { + b0 <- 0 + } + + xs <- (ybars - b0) / b1 + t <- qt(1-0.5*alpha, n - 2) + conf <- t * sxhats + result <- list("Prediction"=xs,"Standard Error"=sxhats, + "Confidence"=conf, "Confidence Limits"=c(xs - conf, xs + conf)) + return(result) +} diff --git a/R/lod.R b/R/lod.R new file mode 100644 index 0000000..5b74418 --- /dev/null +++ b/R/lod.R @@ -0,0 +1,55 @@ +lod <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") +{ + UseMethod("lod") +} + +lod.default <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") +{ + stop("lod is only implemented for univariate lm objects.") +} + +lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") +{ + if (length(object$weights) > 0) { + stop(paste( + "\nThe detemination of a lod from calibration models obtained by", + "weighted linear regression requires confidence intervals for", + "predicted y values taking into account weights for the x values", + "from which the predictions are to be generated.", + "This is not supported by the internally used predict.lm method.", + sep = "\n" + )) + } + xname <- names(object$model)[[2]] + xvalues <- object$model[[2]] + yname <- names(object$model)[[1]] + newdata <- data.frame(0) + names(newdata) <- xname + y0 <- predict(object, newdata, interval = "prediction", + level = 1 - 2 * alpha) + yc <- y0[[1,"upr"]] + if (method == "din") { + y0.d <- predict(object, newdata, interval = "prediction", + level = 1 - 2 * beta) + deltay <- y0.d[[1, "upr"]] - y0.d[[1, "fit"]] + lod.y <- yc + deltay + lod.x <- inverse.predict(object, lod.y)$Prediction + } else { + f <- function(x) { + newdata <- data.frame(x) + names(newdata) <- xname + pi.y <- predict(object, newdata, interval = "prediction", + level = 1 - 2 * beta) + yd <- pi.y[[1,"lwr"]] + (yd - yc)^2 + } + if (tol == "default") tol = min(xvalues[xvalues !=0]) / 1000 + lod.x <- optimize(f, interval = c(0, max(xvalues) * 10), tol = tol)$minimum + newdata <- data.frame(x = lod.x) + names(newdata) <- xname + lod.y <- predict(object, newdata) + } + lod <- list(lod.x, lod.y) + names(lod) <- c(xname, yname) + return(lod) +} diff --git a/R/loq.R b/R/loq.R new file mode 100644 index 0000000..f832265 --- /dev/null +++ b/R/loq.R @@ -0,0 +1,41 @@ +loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto", tol = "default") +{ + UseMethod("loq") +} + +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.") +} + +loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto", tol = "default") +{ + if (length(object$weights) > 0 && var.loq == "auto" && w.loq == "auto") { + stop(paste("If you are using a model from weighted regression,", + "you need to specify a reasonable approximation for the", + "weight (w.loq) or the variance (var.loq) at the", + "limit of quantification")) + } + xname <- names(object$model)[[2]] + xvalues <- object$model[[2]] + yname <- names(object$model)[[1]] + f <- function(x) { + newdata <- data.frame(x = x) + names(newdata) <- xname + y <- predict(object, newdata) + p <- inverse.predict(object, rep(y, n), ws = w.loq, + var.s = var.loq, alpha = alpha) + (p[["Prediction"]] - k * p[["Confidence"]])^2 + } + if (tol == "default") tol = min(xvalues[xvalues !=0]) / 1000 + loq.x <- optimize(f, interval = c(0, max(xvalues) * 10), tol = tol)$minimum + newdata <- data.frame(x = loq.x) + names(newdata) <- xname + loq.y <- predict(object, newdata) + loq <- list(loq.x, loq.y) + names(loq) <- c(xname, yname) + return(loq) +} diff --git a/TODO b/TODO new file mode 100644 index 0000000..1fc44a7 --- /dev/null +++ b/TODO @@ -0,0 +1,5 @@ +- 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/branches/0.1/chemCal/.Rbuildignore b/branches/0.1/chemCal/.Rbuildignore deleted file mode 100644 index 0b413f5..0000000 --- a/branches/0.1/chemCal/.Rbuildignore +++ /dev/null @@ -1,9 +0,0 @@ -GNUmakefile -out$ -toc$ -bbl$ -blg$ -aux$ -log$ -vignettes/chemCal.tex -vignettes/chemCal.pdf diff --git a/branches/0.1/chemCal/ChangeLog b/branches/0.1/chemCal/ChangeLog deleted file mode 100644 index c258c26..0000000 --- a/branches/0.1/chemCal/ChangeLog +++ /dev/null @@ -1,11 +0,0 @@ -2014-04-24 Johannes Ranke for chemCal (0.1-33) - - * 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. - -Changes performed in earlier versions are documented in the subversion log -files found at http://kriemhild.uft.uni-bremen.de/viewcvs/?root=chemCal diff --git a/branches/0.1/chemCal/DESCRIPTION b/branches/0.1/chemCal/DESCRIPTION deleted file mode 100644 index 2608804..0000000 --- a/branches/0.1/chemCal/DESCRIPTION +++ /dev/null @@ -1,18 +0,0 @@ -Package: chemCal -Version: 0.1-35 -Date: 2015-07-02 -Title: Calibration functions for analytical chemistry -Author: Johannes Ranke -Maintainer: Johannes Ranke -Suggests: MASS -Description: chemCal provides simple functions for plotting linear - calibration functions and estimating standard errors for measurements - according to the Handbook of Chemometrics and Qualimetrics: Part A - by Massart et al. There are also functions estimating the limit - of detection (LOD) and limit of quantification (LOQ). - The functions work on model objects from - optionally weighted - linear - regression (lm) or robust linear regression (rlm from the MASS package). -License: GPL (>= 2) -URL: http://www.r-project.org, - http://www.uft.uni-bremen.de/chemie/ranke, - http://kriemhild.uft.uni-bremen.de/viewcvs/?root=chemCal diff --git a/branches/0.1/chemCal/GNUmakefile b/branches/0.1/chemCal/GNUmakefile deleted file mode 100644 index 803c503..0000000 --- a/branches/0.1/chemCal/GNUmakefile +++ /dev/null @@ -1,57 +0,0 @@ -PKGNAME := $(shell sed -n "s/Package: *\([^ ]*\)/\1/p" DESCRIPTION) -PKGVERS := $(shell sed -n "s/Version: *\([^ ]*\)/\1/p" DESCRIPTION) -PKGSRC := $(shell basename $(PWD)) - -# Specify the directory holding R binaries. To use an alternate R build (say a -# pre-prelease version) use `make RBIN=/path/to/other/R/` or `export RBIN=...` -# 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`") - -.PHONY: help - -help: - @echo "\nExecute development tasks for $(PKGNAME)\n" - @echo "Usage: \`make \` where is one of:" - @echo "" - @echo "Development Tasks" - @echo "-----------------" - @echo " build Create the package" - @echo " build-no-vignettes Create the package without rebuilding vignettes" - @echo " check Invoke build and then check the package" - @echo " check-no-vignettes Invoke build without rebuilding vignettes, and then check" - @echo " install Invoke build and then install the result" - @echo " install-no-vignettes Invoke build without rebuilding vignettes and then install the result" - @echo "" - @echo "Using R in: $(RBIN)" - @echo "Set the RBIN environment variable to change this." - @echo "" - - -#------------------------------------------------------------------------------ -# Development Tasks -#------------------------------------------------------------------------------ - -build: - cd ..;\ - "$(RBIN)/R" CMD build $(PKGSRC) - -build-no-vignettes: - cd ..;\ - "$(RBIN)/R" CMD build $(PKGSRC) --no-build-vignettes - -install: build - cd ..;\ - "$(RBIN)/R" CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz - -install-no-vignettes: build-no-vignettes - cd ..;\ - "$(RBIN)/R" CMD INSTALL $(PKGNAME)_$(PKGVERS).tar.gz - -check: build - cd ..;\ - "$(RBIN)/R" CMD check --as-cran $(PKGNAME)_$(PKGVERS).tar.gz - -check-no-vignettes: build-no-vignettes - cd ..;\ - "$(RBIN)/R" CMD check --as-cran $(PKGNAME)_$(PKGVERS).tar.gz diff --git a/branches/0.1/chemCal/NAMESPACE b/branches/0.1/chemCal/NAMESPACE deleted file mode 100644 index 6314e29..0000000 --- a/branches/0.1/chemCal/NAMESPACE +++ /dev/null @@ -1,6 +0,0 @@ -# Export all names -exportPattern(".") -S3method(lod, lm) -S3method(loq, lm) -S3method(inverse.predict, lm) -S3method(inverse.predict, rlm) diff --git a/branches/0.1/chemCal/R/calplot.R b/branches/0.1/chemCal/R/calplot.R deleted file mode 100644 index 6aed9c0..0000000 --- a/branches/0.1/chemCal/R/calplot.R +++ /dev/null @@ -1,80 +0,0 @@ -calplot <- function(object, - xlim = c("auto", "auto"), ylim = c("auto", "auto"), - xlab = "Concentration", ylab = "Response", alpha = 0.05, - varfunc = NULL) -{ - UseMethod("calplot") -} - -calplot.default <- function(object, - xlim = c("auto","auto"), ylim = c("auto","auto"), - xlab = "Concentration", ylab = "Response", - alpha=0.05, varfunc = NULL) -{ - stop("Calibration plots only implemented for univariate lm objects.") -} - -calplot.lm <- function(object, - xlim = c("auto","auto"), ylim = c("auto","auto"), - xlab = "Concentration", ylab = "Response", alpha=0.05, - varfunc = NULL) -{ - if (length(object$coef) > 2) - stop("More than one independent variable in your model - not implemented") - - if (alpha <= 0 | alpha >= 1) - stop("Alpha should be between 0 and 1 (exclusive)") - - m <- object - level <- 1 - alpha - y <- m$model[[1]] - x <- m$model[[2]] - if (xlim[1] == "auto") xlim[1] <- 0 - if (xlim[2] == "auto") xlim[2] <- max(x) - xlim <- as.numeric(xlim) - newdata <- list( - x = seq(from = xlim[[1]], to = xlim[[2]], length=250)) - names(newdata) <- names(m$model)[[2]] - if (is.null(varfunc)) { - varfunc <- if (length(m$weights)) { - function(variable) mean(m$weights) - } else function(variable) rep(1,250) - } - pred.lim <- predict(m, newdata, interval = "prediction", - level=level, weights.newdata = varfunc(m)) - conf.lim <- predict(m, newdata, interval = "confidence", - level=level) - yrange.auto <- range(c(0,pred.lim)) - if (ylim[1] == "auto") ylim[1] <- yrange.auto[1] - if (ylim[2] == "auto") ylim[2] <- yrange.auto[2] - plot(1, - type = "n", - xlab = xlab, - ylab = ylab, - xlim = as.numeric(xlim), - ylim = as.numeric(ylim) - ) - points(x,y, pch = 21, bg = "yellow") - matlines(newdata[[1]], pred.lim, lty = c(1, 4, 4), - col = c("black", "red", "red")) - if (length(object$weights) > 0) { - legend(min(x), - max(pred.lim, na.rm = TRUE), - legend = c("Fitted Line", "Confidence Bands"), - lty = c(1, 3), - lwd = 2, - col = c("black", "green4"), - horiz = FALSE, cex = 0.9, bg = "gray95") - } else { - matlines(newdata[[1]], conf.lim, lty = c(1, 3, 3), - col = c("black", "green4", "green4")) - legend(min(x), - max(pred.lim, na.rm = TRUE), - legend = c("Fitted Line", "Confidence Bands", - "Prediction Bands"), - lty = c(1, 3, 4), - lwd = 2, - col = c("black", "green4", "red"), - horiz = FALSE, cex = 0.9, bg = "gray95") - } -} diff --git a/branches/0.1/chemCal/R/inverse.predict.lm.R b/branches/0.1/chemCal/R/inverse.predict.lm.R deleted file mode 100644 index 927e672..0000000 --- a/branches/0.1/chemCal/R/inverse.predict.lm.R +++ /dev/null @@ -1,97 +0,0 @@ -# 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 - -inverse.predict <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, var.s = "auto") -{ - UseMethod("inverse.predict") -} - -inverse.predict.default <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, var.s = "auto") -{ - stop("Inverse prediction only implemented for univariate lm objects.") -} - -inverse.predict.lm <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, var.s = "auto") -{ - yname = names(object$model)[[1]] - xname = names(object$model)[[2]] - if (ws == "auto") { - ws <- ifelse(length(object$weights) > 0, mean(object$weights), 1) - } - if (length(object$weights) > 0) { - wx <- split(object$weights,object$model[[xname]]) - w <- sapply(wx,mean) - } else { - w <- rep(1,length(split(object$model[[yname]],object$model[[xname]]))) - } - .inverse.predict(object = object, newdata = newdata, - ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) -} - -inverse.predict.rlm <- function(object, newdata, ..., - ws = "auto", alpha = 0.05, var.s = "auto") -{ - yname = names(object$model)[[1]] - xname = names(object$model)[[2]] - if (ws == "auto") { - ws <- mean(object$w) - } - wx <- split(object$weights,object$model[[xname]]) - w <- sapply(wx,mean) - .inverse.predict(object = object, newdata = newdata, - ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname) -} - -.inverse.predict <- function(object, newdata, ws, alpha, var.s, w, xname, yname) -{ - if (length(object$coef) > 2) - stop("More than one independent variable in your model - not implemented") - - if (alpha <= 0 | alpha >= 1) - stop("Alpha should be between 0 and 1 (exclusive)") - - ybars <- mean(newdata) - m <- length(newdata) - - yx <- split(object$model[[yname]], object$model[[xname]]) - n <- length(yx) - df <- n - length(object$coef) - x <- as.numeric(names(yx)) - ybar <- sapply(yx,mean) - yhatx <- split(object$fitted.values, object$model[[xname]]) - yhat <- sapply(yhatx, mean) - se <- sqrt(sum(w * (ybar - yhat)^2)/df) - - if (var.s == "auto") { - var.s <- se^2/ws - } - - b1 <- object$coef[[xname]] - - ybarw <- sum(w * ybar)/sum(w) - -# This is the adapted form of equation 8.28 (see package vignette) - sxhats <- 1/b1 * sqrt( - (var.s / m) + - se^2 * (1/sum(w) + - (ybars - ybarw)^2 * sum(w) / - (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2))) - ) - - if (names(object$coef)[1] == "(Intercept)") { - b0 <- object$coef[["(Intercept)"]] - } else { - b0 <- 0 - } - - xs <- (ybars - b0) / b1 - t <- qt(1-0.5*alpha, n - 2) - conf <- t * sxhats - result <- list("Prediction"=xs,"Standard Error"=sxhats, - "Confidence"=conf, "Confidence Limits"=c(xs - conf, xs + conf)) - return(result) -} diff --git a/branches/0.1/chemCal/R/lod.R b/branches/0.1/chemCal/R/lod.R deleted file mode 100644 index 5b74418..0000000 --- a/branches/0.1/chemCal/R/lod.R +++ /dev/null @@ -1,55 +0,0 @@ -lod <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") -{ - UseMethod("lod") -} - -lod.default <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") -{ - stop("lod is only implemented for univariate lm objects.") -} - -lod.lm <- function(object, ..., alpha = 0.05, beta = 0.05, method = "default", tol = "default") -{ - if (length(object$weights) > 0) { - stop(paste( - "\nThe detemination of a lod from calibration models obtained by", - "weighted linear regression requires confidence intervals for", - "predicted y values taking into account weights for the x values", - "from which the predictions are to be generated.", - "This is not supported by the internally used predict.lm method.", - sep = "\n" - )) - } - xname <- names(object$model)[[2]] - xvalues <- object$model[[2]] - yname <- names(object$model)[[1]] - newdata <- data.frame(0) - names(newdata) <- xname - y0 <- predict(object, newdata, interval = "prediction", - level = 1 - 2 * alpha) - yc <- y0[[1,"upr"]] - if (method == "din") { - y0.d <- predict(object, newdata, interval = "prediction", - level = 1 - 2 * beta) - deltay <- y0.d[[1, "upr"]] - y0.d[[1, "fit"]] - lod.y <- yc + deltay - lod.x <- inverse.predict(object, lod.y)$Prediction - } else { - f <- function(x) { - newdata <- data.frame(x) - names(newdata) <- xname - pi.y <- predict(object, newdata, interval = "prediction", - level = 1 - 2 * beta) - yd <- pi.y[[1,"lwr"]] - (yd - yc)^2 - } - if (tol == "default") tol = min(xvalues[xvalues !=0]) / 1000 - lod.x <- optimize(f, interval = c(0, max(xvalues) * 10), tol = tol)$minimum - newdata <- data.frame(x = lod.x) - names(newdata) <- xname - lod.y <- predict(object, newdata) - } - lod <- list(lod.x, lod.y) - names(lod) <- c(xname, yname) - return(lod) -} diff --git a/branches/0.1/chemCal/R/loq.R b/branches/0.1/chemCal/R/loq.R deleted file mode 100644 index f832265..0000000 --- a/branches/0.1/chemCal/R/loq.R +++ /dev/null @@ -1,41 +0,0 @@ -loq <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", - var.loq = "auto", tol = "default") -{ - UseMethod("loq") -} - -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.") -} - -loq.lm <- function(object, ..., alpha = 0.05, k = 3, n = 1, w.loq = "auto", - var.loq = "auto", tol = "default") -{ - if (length(object$weights) > 0 && var.loq == "auto" && w.loq == "auto") { - stop(paste("If you are using a model from weighted regression,", - "you need to specify a reasonable approximation for the", - "weight (w.loq) or the variance (var.loq) at the", - "limit of quantification")) - } - xname <- names(object$model)[[2]] - xvalues <- object$model[[2]] - yname <- names(object$model)[[1]] - f <- function(x) { - newdata <- data.frame(x = x) - names(newdata) <- xname - y <- predict(object, newdata) - p <- inverse.predict(object, rep(y, n), ws = w.loq, - var.s = var.loq, alpha = alpha) - (p[["Prediction"]] - k * p[["Confidence"]])^2 - } - if (tol == "default") tol = min(xvalues[xvalues !=0]) / 1000 - loq.x <- optimize(f, interval = c(0, max(xvalues) * 10), tol = tol)$minimum - newdata <- data.frame(x = loq.x) - names(newdata) <- xname - loq.y <- predict(object, newdata) - loq <- list(loq.x, loq.y) - names(loq) <- c(xname, yname) - return(loq) -} diff --git a/branches/0.1/chemCal/TODO b/branches/0.1/chemCal/TODO deleted file mode 100644 index 1fc44a7..0000000 --- a/branches/0.1/chemCal/TODO +++ /dev/null @@ -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/branches/0.1/chemCal/data/din32645.rda b/branches/0.1/chemCal/data/din32645.rda deleted file mode 100644 index 4e2913a..0000000 Binary files a/branches/0.1/chemCal/data/din32645.rda and /dev/null differ diff --git a/branches/0.1/chemCal/data/massart97ex1.rda b/branches/0.1/chemCal/data/massart97ex1.rda deleted file mode 100644 index 1a6fd4b..0000000 Binary files a/branches/0.1/chemCal/data/massart97ex1.rda and /dev/null differ diff --git a/branches/0.1/chemCal/data/massart97ex3.rda b/branches/0.1/chemCal/data/massart97ex3.rda deleted file mode 100644 index 7f60f3a..0000000 Binary files a/branches/0.1/chemCal/data/massart97ex3.rda and /dev/null differ diff --git a/branches/0.1/chemCal/man/calplot.lm.Rd b/branches/0.1/chemCal/man/calplot.lm.Rd deleted file mode 100644 index b60a032..0000000 --- a/branches/0.1/chemCal/man/calplot.lm.Rd +++ /dev/null @@ -1,69 +0,0 @@ -\name{calplot} -\alias{calplot} -\alias{calplot.default} -\alias{calplot.lm} -\title{Plot calibration graphs from univariate linear models} -\description{ - Produce graphics of calibration data, the fitted model as well - as confidence, and, for unweighted regression, prediction bands. -} -\usage{ - calplot(object, xlim = c("auto", "auto"), ylim = c("auto", "auto"), - xlab = "Concentration", ylab = "Response", alpha=0.05, varfunc = NULL) -} -\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}. - } - \item{xlim}{ - The limits of the plot on the x axis. - } - \item{ylim}{ - The limits of the plot on the y axis. - } - \item{xlab}{ - The label of the x axis. - } - \item{ylab}{ - The label of the y axis. - } - \item{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). - } - \item{varfunc}{ - The variance function for generating the weights in the model. - Currently, this argument is ignored (see note below). - } -} -\value{ - 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. - -} -\examples{ -data(massart97ex3) -m <- lm(y ~ x, data = massart97ex3) -calplot(m) -} -\author{ - Johannes Ranke - \email{jranke@uni-bremen.de} - \url{http://www.uft.uni-bremen.de/chemie/ranke} -} -\keyword{regression} diff --git a/branches/0.1/chemCal/man/chemCal-package.Rd b/branches/0.1/chemCal/man/chemCal-package.Rd deleted file mode 100644 index fb09cb1..0000000 --- a/branches/0.1/chemCal/man/chemCal-package.Rd +++ /dev/null @@ -1,16 +0,0 @@ -\name{chemCal-package} -\alias{chemCal-package} -\docType{package} -\title{ - Calibration functions for analytical chemistry -} -\description{ - See \url{../DESCRIPTION} -} -\details{ - There is a package vignette located in \url{../doc/chemCal.pdf}. -} -\author{ - Author and Maintainer: Johannes Ranke -} -\keyword{manip} diff --git a/branches/0.1/chemCal/man/din32645.Rd b/branches/0.1/chemCal/man/din32645.Rd deleted file mode 100644 index cacbf07..0000000 --- a/branches/0.1/chemCal/man/din32645.Rd +++ /dev/null @@ -1,61 +0,0 @@ -\name{din32645} -\docType{data} -\alias{din32645} -\title{Calibration data from DIN 32645} -\description{ - Sample dataset to test the package. -} -\usage{data(din32645)} -\format{ - 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)) - -# 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 -} -\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. - \url{http://www.rzuser.uni-heidelberg.de/~df6/download/dintest.htm} - - Currie, L. A. (1997) Nomenclature in evaluation of analytical methods including - detection and quantification capabilities (IUPAC Recommendations 1995). - Analytica Chimica Acta 391, 105 - 126. -} -\keyword{datasets} diff --git a/branches/0.1/chemCal/man/inverse.predict.Rd b/branches/0.1/chemCal/man/inverse.predict.Rd deleted file mode 100644 index 347d670..0000000 --- a/branches/0.1/chemCal/man/inverse.predict.Rd +++ /dev/null @@ -1,69 +0,0 @@ -\name{inverse.predict} -\alias{inverse.predict} -\alias{inverse.predict.lm} -\alias{inverse.predict.rlm} -\alias{inverse.predict.default} -\title{Predict x from y for a linear calibration} -\usage{inverse.predict(object, newdata, \dots, - ws, alpha=0.05, var.s = "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}. - } - \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. This argument is obligatory - if \code{object} has weights. - } - \item{alpha}{ - The error tolerance level for the confidence interval to be reported. - } - \item{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}. - } -} -\value{ - A list containing the predicted x value, its standard error and a - confidence interval. -} -\description{ - 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. -} -\note{ - The function was validated with examples 7 and 8 from Massart et al. (1997). -} -\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 -} -\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 -} -\keyword{manip} diff --git a/branches/0.1/chemCal/man/lod.Rd b/branches/0.1/chemCal/man/lod.Rd deleted file mode 100644 index 04aac1f..0000000 --- a/branches/0.1/chemCal/man/lod.Rd +++ /dev/null @@ -1,88 +0,0 @@ -\name{lod} -\alias{lod} -\alias{lod.lm} -\alias{lod.rlm} -\alias{lod.default} -\title{Estimate a limit of detection (LOD)} -\usage{ - lod(object, \dots, alpha = 0.05, beta = 0.05, method = "default", tol = "default") -} -\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{\dots}{ - Placeholder for further arguments that might be needed by - future implementations. - } - \item{alpha}{ - The error tolerance for the decision limit (critical value). - } - \item{beta}{ - The error tolerance beta for the detection limit. - } - \item{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. - } - \item{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. - } -} -\value{ - A list containig the corresponding x and y values of the estimated limit of - detection of a model used for calibration. -} -\description{ - 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). -} -\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. -} -\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. -} -\examples{ -data(din32645) -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) -} -\seealso{ - Examples for \code{\link{din32645}} -} -\keyword{manip} diff --git a/branches/0.1/chemCal/man/loq.Rd b/branches/0.1/chemCal/man/loq.Rd deleted file mode 100644 index 082cf34..0000000 --- a/branches/0.1/chemCal/man/loq.Rd +++ /dev/null @@ -1,82 +0,0 @@ -\name{loq} -\alias{loq} -\alias{loq.lm} -\alias{loq.rlm} -\alias{loq.default} -\title{Estimate a limit of quantification (LOQ)} -\usage{ - loq(object, \dots, alpha = 0.05, k = 3, n = 1, w.loq = "auto", - var.loq = "auto", tol = "default") -} -\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. If weights are specified - in the model, either \code{w.loq} or \code{var.loq} have to - be specified. - } - \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 LOQ. - } - \item{n}{ - The number of replicate measurements for which the LOQ should be - specified. - } - \item{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. - } - \item{var.loq}{ - The approximate variance at the LOQ. The default value is - calculated from the model. - } - \item{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. - } -} -\value{ - The estimated limit of quantification for a model used for calibration. -} -\description{ - 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. -} -\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. -} -\examples{ -data(massart97ex3) -attach(massart97ex3) -m <- lm(y ~ x) -loq(m) - -# We can get better by using replicate measurements -loq(m, n = 3) -} -\seealso{ - Examples for \code{\link{din32645}} -} -\keyword{manip} diff --git a/branches/0.1/chemCal/man/massart97ex1.Rd b/branches/0.1/chemCal/man/massart97ex1.Rd deleted file mode 100644 index 44e1b85..0000000 --- a/branches/0.1/chemCal/man/massart97ex1.Rd +++ /dev/null @@ -1,17 +0,0 @@ -\name{massart97ex1} -\docType{data} -\alias{massart97ex1} -\title{Calibration data from Massart et al. (1997), example 1} -\description{ - Sample dataset from p. 175 to test the package. -} -\usage{data(massart97ex1)} -\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. -} -\keyword{datasets} diff --git a/branches/0.1/chemCal/man/massart97ex3.Rd b/branches/0.1/chemCal/man/massart97ex3.Rd deleted file mode 100644 index efdcf02..0000000 --- a/branches/0.1/chemCal/man/massart97ex3.Rd +++ /dev/null @@ -1,51 +0,0 @@ -\name{massart97ex3} -\docType{data} -\alias{massart97ex3} -\title{Calibration data from Massart et al. (1997), example 3} -\description{ - Sample dataset from p. 188 to test the package. -} -\usage{data(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) - -# 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 - -# The LOD is only calculated for models from unweighted regression -# with this version of chemCal -m0 <- lm(y ~ x) -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(m, 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. -} -\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. -} -\keyword{datasets} diff --git a/branches/0.1/chemCal/tests/din32645.R b/branches/0.1/chemCal/tests/din32645.R deleted file mode 100644 index e5ffed7..0000000 --- a/branches/0.1/chemCal/tests/din32645.R +++ /dev/null @@ -1,7 +0,0 @@ -require(chemCal) -data(din32645) -m <- lm(y ~ x, data = din32645) -inverse.predict(m, 3500, alpha = 0.01) -lod <- lod(m, alpha = 0.01, beta = 0.5) -lod(m, alpha = 0.01, beta = 0.01) -loq <- loq(m, alpha = 0.01) diff --git a/branches/0.1/chemCal/tests/din32645.Rout.save b/branches/0.1/chemCal/tests/din32645.Rout.save deleted file mode 100644 index 7c9e55d..0000000 --- a/branches/0.1/chemCal/tests/din32645.Rout.save +++ /dev/null @@ -1,48 +0,0 @@ - -R version 3.1.0 (2014-04-10) -- "Spring Dance" -Copyright (C) 2014 The R Foundation for Statistical Computing -Platform: x86_64-pc-linux-gnu (64-bit) - -R is free software and comes with ABSOLUTELY NO WARRANTY. -You are welcome to redistribute it under certain conditions. -Type 'license()' or 'licence()' for distribution details. - -R is a collaborative project with many contributors. -Type 'contributors()' for more information and -'citation()' on how to cite R or R packages in publications. - -Type 'demo()' for some demos, 'help()' for on-line help, or -'help.start()' for an HTML browser interface to help. -Type 'q()' to quit R. - -> require(chemCal) -Loading required package: chemCal -> data(din32645) -> m <- lm(y ~ x, data = din32645) -> inverse.predict(m, 3500, alpha = 0.01) -$Prediction -[1] 0.1054792 - -$`Standard Error` -[1] 0.02215619 - -$Confidence -[1] 0.07434261 - -$`Confidence Limits` -[1] 0.03113656 0.17982178 - -> lod <- lod(m, alpha = 0.01, beta = 0.5) -> lod(m, alpha = 0.01, beta = 0.01) -$x -[1] 0.132909 - -$y - 1 -3765.025 - -> loq <- loq(m, alpha = 0.01) -> -> proc.time() - user system elapsed - 0.472 0.302 0.354 diff --git a/branches/0.1/chemCal/tests/massart97.R b/branches/0.1/chemCal/tests/massart97.R deleted file mode 100644 index 00f837f..0000000 --- a/branches/0.1/chemCal/tests/massart97.R +++ /dev/null @@ -1,25 +0,0 @@ -require(chemCal) -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 - -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) - -inverse.predict(m, 15, ws = 1.67) # 5.9 +- 2.5 -inverse.predict(m, 90, ws = 0.145) # 44.1 +- 7.9 - -m0 <- lm(y ~ x) -lod(m0) - -loq(m0) -loq(m, w.loq = 1.67) diff --git a/branches/0.1/chemCal/tests/massart97.Rout.save b/branches/0.1/chemCal/tests/massart97.Rout.save deleted file mode 100644 index ce99c30..0000000 --- a/branches/0.1/chemCal/tests/massart97.Rout.save +++ /dev/null @@ -1,128 +0,0 @@ - -R version 3.1.0 (2014-04-10) -- "Spring Dance" -Copyright (C) 2014 The R Foundation for Statistical Computing -Platform: x86_64-pc-linux-gnu (64-bit) - -R is free software and comes with ABSOLUTELY NO WARRANTY. -You are welcome to redistribute it under certain conditions. -Type 'license()' or 'licence()' for distribution details. - -R is a collaborative project with many contributors. -Type 'contributors()' for more information and -'citation()' on how to cite R or R packages in publications. - -Type 'demo()' for some demos, 'help()' for on-line help, or -'help.start()' for an HTML browser interface to help. -Type 'q()' to quit R. - -> require(chemCal) -Loading required package: chemCal -> data(massart97ex1) -> m <- lm(y ~ x, data = massart97ex1) -> inverse.predict(m, 15) # 6.1 +- 4.9 -$Prediction -[1] 6.09381 - -$`Standard Error` -[1] 1.767278 - -$Confidence -[1] 4.906751 - -$`Confidence Limits` -[1] 1.187059 11.000561 - -> inverse.predict(m, 90) # 43.9 +- 4.9 -$Prediction -[1] 43.93983 - -$`Standard Error` -[1] 1.767747 - -$Confidence -[1] 4.908053 - -$`Confidence Limits` -[1] 39.03178 48.84788 - -> inverse.predict(m, rep(90,5)) # 43.9 +- 3.2 -$Prediction -[1] 43.93983 - -$`Standard Error` -[1] 1.141204 - -$Confidence -[1] 3.168489 - -$`Confidence Limits` -[1] 40.77134 47.10832 - -> -> 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) -> -> inverse.predict(m, 15, ws = 1.67) # 5.9 +- 2.5 -$Prediction -[1] 5.865367 - -$`Standard Error` -[1] 0.8926109 - -$Confidence -[1] 2.478285 - -$`Confidence Limits` -[1] 3.387082 8.343652 - -> inverse.predict(m, 90, ws = 0.145) # 44.1 +- 7.9 -$Prediction -[1] 44.06025 - -$`Standard Error` -[1] 2.829162 - -$Confidence -[1] 7.855012 - -$`Confidence Limits` -[1] 36.20523 51.91526 - -> -> m0 <- lm(y ~ x) -> lod(m0) -$x -[1] 5.407085 - -$y - 1 -13.63911 - -> -> loq(m0) -$x -[1] 13.97764 - -$y - 1 -30.6235 - -> loq(m, w.loq = 1.67) -$x -[1] 7.346195 - -$y - 1 -17.90777 - -> -> proc.time() - user system elapsed - 0.529 0.327 0.443 diff --git a/branches/0.1/chemCal/vignettes/chemCal-001.pdf b/branches/0.1/chemCal/vignettes/chemCal-001.pdf deleted file mode 100644 index ef9e332..0000000 Binary files a/branches/0.1/chemCal/vignettes/chemCal-001.pdf and /dev/null differ diff --git a/branches/0.1/chemCal/vignettes/chemCal-002.pdf b/branches/0.1/chemCal/vignettes/chemCal-002.pdf deleted file mode 100644 index 2b56f14..0000000 Binary files a/branches/0.1/chemCal/vignettes/chemCal-002.pdf and /dev/null differ diff --git a/branches/0.1/chemCal/vignettes/chemCal.Rnw b/branches/0.1/chemCal/vignettes/chemCal.Rnw deleted file mode 100644 index 30e0331..0000000 --- a/branches/0.1/chemCal/vignettes/chemCal.Rnw +++ /dev/null @@ -1,130 +0,0 @@ -\documentclass[a4paper]{article} -%\VignetteIndexEntry{Short manual for the chemCal package} -\usepackage{hyperref} - -\title{Basic calibration functions for analytical chemistry} -\author{Johannes Ranke} - -\begin{document} -\maketitle - -The \texttt{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 -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 \cite{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 \texttt{lm} or \texttt{rlm}, plus to datasets for -validation. - -A \href{http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877}{bug -report (PR\#8877)} 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. - -When calibrating an analytical method, the first task is to generate a suitable -model. If we want to use the \texttt{chemCal} functions, we will have to restrict -ourselves to univariate, possibly weighted, linear regression so far. - -Once such a model has been created, the calibration can be graphically -shown by using the \texttt{calplot} function: - -<>= -library(chemCal) -data(massart97ex3) -m0 <- lm(y ~ x, data = massart97ex3) -calplot(m0) -@ - -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: - -<>= -plot(m0,which=3) -@ - -Therefore, in Example 8 in \cite{massart97} weighted regression -is proposed which can be reproduced by - -<<>>= -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) -@ - -If we now want to predict a new x value from measured y values, -we use the \texttt{inverse.predict} function: - -<<>>= -inverse.predict(m, 15, ws=1.67) -inverse.predict(m, 90, ws = 0.145) -@ - -The weight \texttt{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 \texttt{var.s} at this location. - -\section*{Theory for \texttt{inverse.predict}} -Equation 8.28 in \cite{massart97} gives a general equation for predicting the -standard error $s_{\hat{x_s}}$ for an x value predicted from measurements of y -according to the linear calibration function $ y = b_0 + b_1 \cdot x$: - -\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} - -with - -\begin{equation} -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 mean $y$ -value (!) observed for standard $i$, $\hat{y_i}$ is the estimated value for -standard $i$, $n$ is the number calibration standards, $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 -precisions in standards and samples into account. In analogy to Equation 8.26 -from \cite{massart97} we get - -\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} - -where I interpret $\frac{{s_s}^2}{w_s}$ as an estimator of the variance at location -$\hat{x_s}$, which can be replaced by a user-specified value using the argument -\texttt{var.s} of the function \texttt{inverse.predict}. - -\begin{thebibliography}{1} -\bibitem{massart97} -Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J., -Smeyers-Verbeke, J. -\newblock Handbook of Chemometrics and Qualimetrics: Part A, -\newblock Elsevier, Amsterdam, 1997 -\end{thebibliography} - -\end{document} diff --git a/branches/0.1/chemCal/vignettes/chemCal.pdf b/branches/0.1/chemCal/vignettes/chemCal.pdf deleted file mode 100644 index 935b160..0000000 Binary files a/branches/0.1/chemCal/vignettes/chemCal.pdf and /dev/null differ diff --git a/branches/0.1/chemCal/vignettes/chemCal.tex b/branches/0.1/chemCal/vignettes/chemCal.tex deleted file mode 100644 index 6326126..0000000 --- a/branches/0.1/chemCal/vignettes/chemCal.tex +++ /dev/null @@ -1,169 +0,0 @@ -\documentclass[a4paper]{article} -%\VignetteIndexEntry{Short manual for the chemCal package} -\usepackage{hyperref} - -\title{Basic calibration functions for analytical chemistry} -\author{Johannes Ranke} - -\usepackage{Sweave} -\begin{document} -\maketitle - -The \texttt{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 -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 \cite{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 \texttt{lm} or \texttt{rlm}, plus to datasets for -validation. - -A \href{http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877}{bug -report (PR\#8877)} 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. - -When calibrating an analytical method, the first task is to generate a suitable -model. If we want to use the \texttt{chemCal} functions, we will have to restrict -ourselves to univariate, possibly weighted, linear regression so far. - -Once such a model has been created, the calibration can be graphically -shown by using the \texttt{calplot} function: - -\begin{Schunk} -\begin{Sinput} -> library(chemCal) -> data(massart97ex3) -> m0 <- lm(y ~ x, data = massart97ex3) -> calplot(m0) -\end{Sinput} -\end{Schunk} -\includegraphics{chemCal-001} - -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: - -\begin{Schunk} -\begin{Sinput} -> plot(m0,which=3) -\end{Sinput} -\end{Schunk} -\includegraphics{chemCal-002} - -Therefore, in Example 8 in \cite{massart97} weighted regression -is proposed which can be reproduced by - -\begin{Schunk} -\begin{Sinput} -> 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) -\end{Sinput} -\end{Schunk} - -If we now want to predict a new x value from measured y values, -we use the \texttt{inverse.predict} function: - -\begin{Schunk} -\begin{Sinput} -> inverse.predict(m, 15, ws=1.67) -\end{Sinput} -\begin{Soutput} -$Prediction -[1] 5.865367 - -$`Standard Error` -[1] 0.8926109 - -$Confidence -[1] 2.478285 - -$`Confidence Limits` -[1] 3.387082 8.343652 -\end{Soutput} -\begin{Sinput} -> inverse.predict(m, 90, ws = 0.145) -\end{Sinput} -\begin{Soutput} -$Prediction -[1] 44.06025 - -$`Standard Error` -[1] 2.829162 - -$Confidence -[1] 7.855012 - -$`Confidence Limits` -[1] 36.20523 51.91526 -\end{Soutput} -\end{Schunk} - -The weight \texttt{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 \texttt{var.s} at this location. - -\section*{Theory for \texttt{inverse.predict}} -Equation 8.28 in \cite{massart97} gives a general equation for predicting the -standard error $s_{\hat{x_s}}$ for an x value predicted from measurements of y -according to the linear calibration function $ y = b_0 + b_1 \cdot x$: - -\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} - -with - -\begin{equation} -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 mean $y$ -value (!) observed for standard $i$, $\hat{y_i}$ is the estimated value for -standard $i$, $n$ is the number calibration standards, $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 -precisions in standards and samples into account. In analogy to Equation 8.26 -from \cite{massart97} we get - -\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} - -where I interpret $\frac{{s_s}^2}{w_s}$ as an estimator of the variance at location -$\hat{x_s}$, which can be replaced by a user-specified value using the argument -\texttt{var.s} of the function \texttt{inverse.predict}. - -\begin{thebibliography}{1} -\bibitem{massart97} -Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J., -Smeyers-Verbeke, J. -\newblock Handbook of Chemometrics and Qualimetrics: Part A, -\newblock Elsevier, Amsterdam, 1997 -\end{thebibliography} - -\end{document} diff --git a/data/din32645.rda b/data/din32645.rda new file mode 100644 index 0000000..4e2913a Binary files /dev/null and b/data/din32645.rda differ diff --git a/data/massart97ex1.rda b/data/massart97ex1.rda new file mode 100644 index 0000000..1a6fd4b Binary files /dev/null and b/data/massart97ex1.rda differ diff --git a/data/massart97ex3.rda b/data/massart97ex3.rda new file mode 100644 index 0000000..7f60f3a Binary files /dev/null and b/data/massart97ex3.rda differ diff --git a/man/calplot.lm.Rd b/man/calplot.lm.Rd new file mode 100644 index 0000000..b60a032 --- /dev/null +++ b/man/calplot.lm.Rd @@ -0,0 +1,69 @@ +\name{calplot} +\alias{calplot} +\alias{calplot.default} +\alias{calplot.lm} +\title{Plot calibration graphs from univariate linear models} +\description{ + Produce graphics of calibration data, the fitted model as well + as confidence, and, for unweighted regression, prediction bands. +} +\usage{ + calplot(object, xlim = c("auto", "auto"), ylim = c("auto", "auto"), + xlab = "Concentration", ylab = "Response", alpha=0.05, varfunc = NULL) +} +\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}. + } + \item{xlim}{ + The limits of the plot on the x axis. + } + \item{ylim}{ + The limits of the plot on the y axis. + } + \item{xlab}{ + The label of the x axis. + } + \item{ylab}{ + The label of the y axis. + } + \item{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). + } + \item{varfunc}{ + The variance function for generating the weights in the model. + Currently, this argument is ignored (see note below). + } +} +\value{ + 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. + +} +\examples{ +data(massart97ex3) +m <- lm(y ~ x, data = massart97ex3) +calplot(m) +} +\author{ + Johannes Ranke + \email{jranke@uni-bremen.de} + \url{http://www.uft.uni-bremen.de/chemie/ranke} +} +\keyword{regression} diff --git a/man/chemCal-package.Rd b/man/chemCal-package.Rd new file mode 100644 index 0000000..fb09cb1 --- /dev/null +++ b/man/chemCal-package.Rd @@ -0,0 +1,16 @@ +\name{chemCal-package} +\alias{chemCal-package} +\docType{package} +\title{ + Calibration functions for analytical chemistry +} +\description{ + See \url{../DESCRIPTION} +} +\details{ + There is a package vignette located in \url{../doc/chemCal.pdf}. +} +\author{ + Author and Maintainer: Johannes Ranke +} +\keyword{manip} diff --git a/man/din32645.Rd b/man/din32645.Rd new file mode 100644 index 0000000..cacbf07 --- /dev/null +++ b/man/din32645.Rd @@ -0,0 +1,61 @@ +\name{din32645} +\docType{data} +\alias{din32645} +\title{Calibration data from DIN 32645} +\description{ + Sample dataset to test the package. +} +\usage{data(din32645)} +\format{ + 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)) + +# 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 +} +\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. + \url{http://www.rzuser.uni-heidelberg.de/~df6/download/dintest.htm} + + Currie, L. A. (1997) Nomenclature in evaluation of analytical methods including + detection and quantification capabilities (IUPAC Recommendations 1995). + Analytica Chimica Acta 391, 105 - 126. +} +\keyword{datasets} diff --git a/man/inverse.predict.Rd b/man/inverse.predict.Rd new file mode 100644 index 0000000..347d670 --- /dev/null +++ b/man/inverse.predict.Rd @@ -0,0 +1,69 @@ +\name{inverse.predict} +\alias{inverse.predict} +\alias{inverse.predict.lm} +\alias{inverse.predict.rlm} +\alias{inverse.predict.default} +\title{Predict x from y for a linear calibration} +\usage{inverse.predict(object, newdata, \dots, + ws, alpha=0.05, var.s = "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}. + } + \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. This argument is obligatory + if \code{object} has weights. + } + \item{alpha}{ + The error tolerance level for the confidence interval to be reported. + } + \item{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}. + } +} +\value{ + A list containing the predicted x value, its standard error and a + confidence interval. +} +\description{ + 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. +} +\note{ + The function was validated with examples 7 and 8 from Massart et al. (1997). +} +\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 +} +\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 +} +\keyword{manip} diff --git a/man/lod.Rd b/man/lod.Rd new file mode 100644 index 0000000..04aac1f --- /dev/null +++ b/man/lod.Rd @@ -0,0 +1,88 @@ +\name{lod} +\alias{lod} +\alias{lod.lm} +\alias{lod.rlm} +\alias{lod.default} +\title{Estimate a limit of detection (LOD)} +\usage{ + lod(object, \dots, alpha = 0.05, beta = 0.05, method = "default", tol = "default") +} +\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{\dots}{ + Placeholder for further arguments that might be needed by + future implementations. + } + \item{alpha}{ + The error tolerance for the decision limit (critical value). + } + \item{beta}{ + The error tolerance beta for the detection limit. + } + \item{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. + } + \item{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. + } +} +\value{ + A list containig the corresponding x and y values of the estimated limit of + detection of a model used for calibration. +} +\description{ + 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). +} +\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. +} +\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. +} +\examples{ +data(din32645) +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) +} +\seealso{ + Examples for \code{\link{din32645}} +} +\keyword{manip} diff --git a/man/loq.Rd b/man/loq.Rd new file mode 100644 index 0000000..082cf34 --- /dev/null +++ b/man/loq.Rd @@ -0,0 +1,82 @@ +\name{loq} +\alias{loq} +\alias{loq.lm} +\alias{loq.rlm} +\alias{loq.default} +\title{Estimate a limit of quantification (LOQ)} +\usage{ + loq(object, \dots, alpha = 0.05, k = 3, n = 1, w.loq = "auto", + var.loq = "auto", tol = "default") +} +\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. If weights are specified + in the model, either \code{w.loq} or \code{var.loq} have to + be specified. + } + \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 LOQ. + } + \item{n}{ + The number of replicate measurements for which the LOQ should be + specified. + } + \item{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. + } + \item{var.loq}{ + The approximate variance at the LOQ. The default value is + calculated from the model. + } + \item{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. + } +} +\value{ + The estimated limit of quantification for a model used for calibration. +} +\description{ + 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. +} +\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. +} +\examples{ +data(massart97ex3) +attach(massart97ex3) +m <- lm(y ~ x) +loq(m) + +# We can get better by using replicate measurements +loq(m, n = 3) +} +\seealso{ + Examples for \code{\link{din32645}} +} +\keyword{manip} diff --git a/man/massart97ex1.Rd b/man/massart97ex1.Rd new file mode 100644 index 0000000..44e1b85 --- /dev/null +++ b/man/massart97ex1.Rd @@ -0,0 +1,17 @@ +\name{massart97ex1} +\docType{data} +\alias{massart97ex1} +\title{Calibration data from Massart et al. (1997), example 1} +\description{ + Sample dataset from p. 175 to test the package. +} +\usage{data(massart97ex1)} +\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. +} +\keyword{datasets} diff --git a/man/massart97ex3.Rd b/man/massart97ex3.Rd new file mode 100644 index 0000000..efdcf02 --- /dev/null +++ b/man/massart97ex3.Rd @@ -0,0 +1,51 @@ +\name{massart97ex3} +\docType{data} +\alias{massart97ex3} +\title{Calibration data from Massart et al. (1997), example 3} +\description{ + Sample dataset from p. 188 to test the package. +} +\usage{data(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) + +# 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 + +# The LOD is only calculated for models from unweighted regression +# with this version of chemCal +m0 <- lm(y ~ x) +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(m, 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. +} +\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. +} +\keyword{datasets} diff --git a/tests/din32645.R b/tests/din32645.R new file mode 100644 index 0000000..e5ffed7 --- /dev/null +++ b/tests/din32645.R @@ -0,0 +1,7 @@ +require(chemCal) +data(din32645) +m <- lm(y ~ x, data = din32645) +inverse.predict(m, 3500, alpha = 0.01) +lod <- lod(m, alpha = 0.01, beta = 0.5) +lod(m, alpha = 0.01, beta = 0.01) +loq <- loq(m, alpha = 0.01) diff --git a/tests/din32645.Rout.save b/tests/din32645.Rout.save new file mode 100644 index 0000000..7c9e55d --- /dev/null +++ b/tests/din32645.Rout.save @@ -0,0 +1,48 @@ + +R version 3.1.0 (2014-04-10) -- "Spring Dance" +Copyright (C) 2014 The R Foundation for Statistical Computing +Platform: x86_64-pc-linux-gnu (64-bit) + +R is free software and comes with ABSOLUTELY NO WARRANTY. +You are welcome to redistribute it under certain conditions. +Type 'license()' or 'licence()' for distribution details. + +R is a collaborative project with many contributors. +Type 'contributors()' for more information and +'citation()' on how to cite R or R packages in publications. + +Type 'demo()' for some demos, 'help()' for on-line help, or +'help.start()' for an HTML browser interface to help. +Type 'q()' to quit R. + +> require(chemCal) +Loading required package: chemCal +> data(din32645) +> m <- lm(y ~ x, data = din32645) +> inverse.predict(m, 3500, alpha = 0.01) +$Prediction +[1] 0.1054792 + +$`Standard Error` +[1] 0.02215619 + +$Confidence +[1] 0.07434261 + +$`Confidence Limits` +[1] 0.03113656 0.17982178 + +> lod <- lod(m, alpha = 0.01, beta = 0.5) +> lod(m, alpha = 0.01, beta = 0.01) +$x +[1] 0.132909 + +$y + 1 +3765.025 + +> loq <- loq(m, alpha = 0.01) +> +> proc.time() + user system elapsed + 0.472 0.302 0.354 diff --git a/tests/massart97.R b/tests/massart97.R new file mode 100644 index 0000000..00f837f --- /dev/null +++ b/tests/massart97.R @@ -0,0 +1,25 @@ +require(chemCal) +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 + +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) + +inverse.predict(m, 15, ws = 1.67) # 5.9 +- 2.5 +inverse.predict(m, 90, ws = 0.145) # 44.1 +- 7.9 + +m0 <- lm(y ~ x) +lod(m0) + +loq(m0) +loq(m, w.loq = 1.67) diff --git a/tests/massart97.Rout.save b/tests/massart97.Rout.save new file mode 100644 index 0000000..ce99c30 --- /dev/null +++ b/tests/massart97.Rout.save @@ -0,0 +1,128 @@ + +R version 3.1.0 (2014-04-10) -- "Spring Dance" +Copyright (C) 2014 The R Foundation for Statistical Computing +Platform: x86_64-pc-linux-gnu (64-bit) + +R is free software and comes with ABSOLUTELY NO WARRANTY. +You are welcome to redistribute it under certain conditions. +Type 'license()' or 'licence()' for distribution details. + +R is a collaborative project with many contributors. +Type 'contributors()' for more information and +'citation()' on how to cite R or R packages in publications. + +Type 'demo()' for some demos, 'help()' for on-line help, or +'help.start()' for an HTML browser interface to help. +Type 'q()' to quit R. + +> require(chemCal) +Loading required package: chemCal +> data(massart97ex1) +> m <- lm(y ~ x, data = massart97ex1) +> inverse.predict(m, 15) # 6.1 +- 4.9 +$Prediction +[1] 6.09381 + +$`Standard Error` +[1] 1.767278 + +$Confidence +[1] 4.906751 + +$`Confidence Limits` +[1] 1.187059 11.000561 + +> inverse.predict(m, 90) # 43.9 +- 4.9 +$Prediction +[1] 43.93983 + +$`Standard Error` +[1] 1.767747 + +$Confidence +[1] 4.908053 + +$`Confidence Limits` +[1] 39.03178 48.84788 + +> inverse.predict(m, rep(90,5)) # 43.9 +- 3.2 +$Prediction +[1] 43.93983 + +$`Standard Error` +[1] 1.141204 + +$Confidence +[1] 3.168489 + +$`Confidence Limits` +[1] 40.77134 47.10832 + +> +> 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) +> +> inverse.predict(m, 15, ws = 1.67) # 5.9 +- 2.5 +$Prediction +[1] 5.865367 + +$`Standard Error` +[1] 0.8926109 + +$Confidence +[1] 2.478285 + +$`Confidence Limits` +[1] 3.387082 8.343652 + +> inverse.predict(m, 90, ws = 0.145) # 44.1 +- 7.9 +$Prediction +[1] 44.06025 + +$`Standard Error` +[1] 2.829162 + +$Confidence +[1] 7.855012 + +$`Confidence Limits` +[1] 36.20523 51.91526 + +> +> m0 <- lm(y ~ x) +> lod(m0) +$x +[1] 5.407085 + +$y + 1 +13.63911 + +> +> loq(m0) +$x +[1] 13.97764 + +$y + 1 +30.6235 + +> loq(m, w.loq = 1.67) +$x +[1] 7.346195 + +$y + 1 +17.90777 + +> +> proc.time() + user system elapsed + 0.529 0.327 0.443 diff --git a/vignettes/chemCal-001.pdf b/vignettes/chemCal-001.pdf new file mode 100644 index 0000000..ef9e332 Binary files /dev/null and b/vignettes/chemCal-001.pdf differ diff --git a/vignettes/chemCal-002.pdf b/vignettes/chemCal-002.pdf new file mode 100644 index 0000000..2b56f14 Binary files /dev/null and b/vignettes/chemCal-002.pdf differ diff --git a/vignettes/chemCal.Rnw b/vignettes/chemCal.Rnw new file mode 100644 index 0000000..30e0331 --- /dev/null +++ b/vignettes/chemCal.Rnw @@ -0,0 +1,130 @@ +\documentclass[a4paper]{article} +%\VignetteIndexEntry{Short manual for the chemCal package} +\usepackage{hyperref} + +\title{Basic calibration functions for analytical chemistry} +\author{Johannes Ranke} + +\begin{document} +\maketitle + +The \texttt{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 +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 \cite{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 \texttt{lm} or \texttt{rlm}, plus to datasets for +validation. + +A \href{http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877}{bug +report (PR\#8877)} 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. + +When calibrating an analytical method, the first task is to generate a suitable +model. If we want to use the \texttt{chemCal} functions, we will have to restrict +ourselves to univariate, possibly weighted, linear regression so far. + +Once such a model has been created, the calibration can be graphically +shown by using the \texttt{calplot} function: + +<>= +library(chemCal) +data(massart97ex3) +m0 <- lm(y ~ x, data = massart97ex3) +calplot(m0) +@ + +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: + +<>= +plot(m0,which=3) +@ + +Therefore, in Example 8 in \cite{massart97} weighted regression +is proposed which can be reproduced by + +<<>>= +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) +@ + +If we now want to predict a new x value from measured y values, +we use the \texttt{inverse.predict} function: + +<<>>= +inverse.predict(m, 15, ws=1.67) +inverse.predict(m, 90, ws = 0.145) +@ + +The weight \texttt{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 \texttt{var.s} at this location. + +\section*{Theory for \texttt{inverse.predict}} +Equation 8.28 in \cite{massart97} gives a general equation for predicting the +standard error $s_{\hat{x_s}}$ for an x value predicted from measurements of y +according to the linear calibration function $ y = b_0 + b_1 \cdot x$: + +\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} + +with + +\begin{equation} +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 mean $y$ +value (!) observed for standard $i$, $\hat{y_i}$ is the estimated value for +standard $i$, $n$ is the number calibration standards, $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 +precisions in standards and samples into account. In analogy to Equation 8.26 +from \cite{massart97} we get + +\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} + +where I interpret $\frac{{s_s}^2}{w_s}$ as an estimator of the variance at location +$\hat{x_s}$, which can be replaced by a user-specified value using the argument +\texttt{var.s} of the function \texttt{inverse.predict}. + +\begin{thebibliography}{1} +\bibitem{massart97} +Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J., +Smeyers-Verbeke, J. +\newblock Handbook of Chemometrics and Qualimetrics: Part A, +\newblock Elsevier, Amsterdam, 1997 +\end{thebibliography} + +\end{document} diff --git a/vignettes/chemCal.pdf b/vignettes/chemCal.pdf new file mode 100644 index 0000000..935b160 Binary files /dev/null and b/vignettes/chemCal.pdf differ diff --git a/vignettes/chemCal.tex b/vignettes/chemCal.tex new file mode 100644 index 0000000..6326126 --- /dev/null +++ b/vignettes/chemCal.tex @@ -0,0 +1,169 @@ +\documentclass[a4paper]{article} +%\VignetteIndexEntry{Short manual for the chemCal package} +\usepackage{hyperref} + +\title{Basic calibration functions for analytical chemistry} +\author{Johannes Ranke} + +\usepackage{Sweave} +\begin{document} +\maketitle + +The \texttt{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 +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 \cite{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 \texttt{lm} or \texttt{rlm}, plus to datasets for +validation. + +A \href{http://bugs.r-project.org/bugzilla3/show_bug.cgi?id=8877}{bug +report (PR\#8877)} 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. + +When calibrating an analytical method, the first task is to generate a suitable +model. If we want to use the \texttt{chemCal} functions, we will have to restrict +ourselves to univariate, possibly weighted, linear regression so far. + +Once such a model has been created, the calibration can be graphically +shown by using the \texttt{calplot} function: + +\begin{Schunk} +\begin{Sinput} +> library(chemCal) +> data(massart97ex3) +> m0 <- lm(y ~ x, data = massart97ex3) +> calplot(m0) +\end{Sinput} +\end{Schunk} +\includegraphics{chemCal-001} + +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: + +\begin{Schunk} +\begin{Sinput} +> plot(m0,which=3) +\end{Sinput} +\end{Schunk} +\includegraphics{chemCal-002} + +Therefore, in Example 8 in \cite{massart97} weighted regression +is proposed which can be reproduced by + +\begin{Schunk} +\begin{Sinput} +> 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) +\end{Sinput} +\end{Schunk} + +If we now want to predict a new x value from measured y values, +we use the \texttt{inverse.predict} function: + +\begin{Schunk} +\begin{Sinput} +> inverse.predict(m, 15, ws=1.67) +\end{Sinput} +\begin{Soutput} +$Prediction +[1] 5.865367 + +$`Standard Error` +[1] 0.8926109 + +$Confidence +[1] 2.478285 + +$`Confidence Limits` +[1] 3.387082 8.343652 +\end{Soutput} +\begin{Sinput} +> inverse.predict(m, 90, ws = 0.145) +\end{Sinput} +\begin{Soutput} +$Prediction +[1] 44.06025 + +$`Standard Error` +[1] 2.829162 + +$Confidence +[1] 7.855012 + +$`Confidence Limits` +[1] 36.20523 51.91526 +\end{Soutput} +\end{Schunk} + +The weight \texttt{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 \texttt{var.s} at this location. + +\section*{Theory for \texttt{inverse.predict}} +Equation 8.28 in \cite{massart97} gives a general equation for predicting the +standard error $s_{\hat{x_s}}$ for an x value predicted from measurements of y +according to the linear calibration function $ y = b_0 + b_1 \cdot x$: + +\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} + +with + +\begin{equation} +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 mean $y$ +value (!) observed for standard $i$, $\hat{y_i}$ is the estimated value for +standard $i$, $n$ is the number calibration standards, $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 +precisions in standards and samples into account. In analogy to Equation 8.26 +from \cite{massart97} we get + +\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} + +where I interpret $\frac{{s_s}^2}{w_s}$ as an estimator of the variance at location +$\hat{x_s}$, which can be replaced by a user-specified value using the argument +\texttt{var.s} of the function \texttt{inverse.predict}. + +\begin{thebibliography}{1} +\bibitem{massart97} +Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J., +Smeyers-Verbeke, J. +\newblock Handbook of Chemometrics and Qualimetrics: Part A, +\newblock Elsevier, Amsterdam, 1997 +\end{thebibliography} + +\end{document} -- cgit v1.2.1