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 --- 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 ------------------------- 30 files changed, 1339 deletions(-) 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 (limited to 'branches/0.1') 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} -- cgit v1.2.1