From 965af33bfc386b0c96a50c85fbddf98211e266c4 Mon Sep 17 00:00:00 2001 From: ranke Date: Tue, 15 Feb 2005 19:15:54 +0000 Subject: Cleaned up version, only containing very basic stuff. git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@2 5fad18fb-23f0-0310-ab10-e59a3bee62b4 --- DESCRIPTION | 2 +- R/chemCal.R | 132 ++++++----------------------------------------- data/pahCalibration.rda | Bin 1289 -> 0 bytes data/pahMeasurements.rda | Bin 1348 -> 0 bytes man/calm.Rd | 43 +++++++++++++++ man/calplot.Rd | 52 ------------------- man/calpredict.Rd | 64 ----------------------- man/din32645.Rd | 2 +- man/multical.Rd | 37 ------------- man/pahCalibration.Rd | 19 ------- man/pahMeasurements.Rd | 19 ------- man/plot.calm.Rd | 48 +++++++++++++++++ man/predictx.Rd | 37 +++++++++++++ 13 files changed, 146 insertions(+), 309 deletions(-) delete mode 100644 data/pahCalibration.rda delete mode 100644 data/pahMeasurements.rda create mode 100644 man/calm.Rd delete mode 100644 man/calplot.Rd delete mode 100644 man/calpredict.Rd delete mode 100644 man/multical.Rd delete mode 100644 man/pahCalibration.Rd delete mode 100644 man/pahMeasurements.Rd create mode 100644 man/plot.calm.Rd create mode 100644 man/predictx.Rd diff --git a/DESCRIPTION b/DESCRIPTION index 18bcf09..115bd45 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: chemCal -Version: 0.03-1 +Version: 0.04-1 Date: 2005-01-14 Title: Calibration for analytical chemistry Author: Johannes Ranke diff --git a/R/chemCal.R b/R/chemCal.R index 5de33f8..fab7db4 100644 --- a/R/chemCal.R +++ b/R/chemCal.R @@ -20,10 +20,11 @@ calm <- function(data) predict.calm <- predict.lm print.calm <- print.lm summary.calm <- summary.lm -plot.calm <- function(m, +plot.calm <- function(x,..., xunit="",yunit="",measurand="", level=0.95) { + m <- x x <- m$model$x y <- m$model$y newdata <- data.frame(x = seq(0,max(x),length=250)) @@ -59,137 +60,36 @@ plot.calm <- function(m, col = c("black", "green4", "red"), horiz = FALSE, cex = 0.9, bg = "gray95") } -predictx.calm <- function(m,measurements) +predictx <- function(m,yobs,level=0.95) { s <- summary(m) xi <- m$model$x yi <- m$model$y n <- length(yi) - yobs <- newdata[[1]] p <- length(yobs) - if (!m$intercept) - { - varb1 <- summary(m)$coef["x","Std. Error"]^2 - xpred <- mean(yobs)/b1 - varxpred <- (varyobs + xpred^2 * varb1) / b1^2 - sdxpred <- sqrt(varxpred) - } else - { - b0 <- summary(m)$coef["(Intercept)","Estimate"] - b1 <- summary(m)$coef["xi","Estimate"] - xpred <- (mean(yobs) - b0)/b1 - sumxxbar <- sum((xi - mean(xi))^2) - if (!syobs) - { - yybar <- (mean(yobs) - mean(yi))^2 - sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5 - } else - { - sdxpred <- ((varyobs^0.5/b1) + (S/b1)^2 * (1/n + ((xpred - mean(xi))^2)/sumxxbar))^0.5 - } - } - t <- qt((1 + level)/2,ntot - 2) - confxpred <- t * sdxpred - - result <- c(estimate=xpred,sdxpred=sdxpred,syobs=syobs, - confxpred=confxpred) - digits <- max(c(3,round(log10(xpred/confxpred)) + 2)) - roundedresult <- round(result,digits=digits) - confidenceinterval <- paste(roundedresult["estimate"],"+-", - roundedresult["confxpred"],xunit) - roundedresult[["confidenceinterval"]] <- confidenceinterval - invisible(roundedresult) -} -calpredict <- function(yobs,xi,yi,xunit="",level=0.95,intercept=FALSE,syobs=FALSE) -{ - if (length(xi)!=length(yi)) - { - cat("xi and yi have to be of the same length\n") - } - xi <- xi[!is.na(yi)] - yi <- yi[!is.na(yi)] - n <- length(yi) - p <- length(yobs) - if (!intercept) - { - m <- lm(yi ~ xi - 1) - } else - { - m <- lm(yi ~ xi) - } - S <- summary(m)$sigma - b1 <- summary(m)$coef["xi","Estimate"] - - if (syobs) - { - if (is.numeric(syobs)) - { - varyobs <- syobs^2 - ntot <- n - } else - { - if (length(yobs)==1) - { - cat("yobs has to contain more than one number vector, if you use syobs=TRUE\n") - } - varyobs <- var(yobs) - ntot <- n + p - } - } else - { - varyobs <- S^2 - ntot <- n + if (p > 1) { + varyobs <- var(yobs) + } else { + varyobs <- 0 } - - if (!intercept) - { - varb1 <- summary(m)$coef["xi","Std. Error"]^2 + if (!m$intercept) { + b1 <- summary(m)$coef["x","Estimate"] + varb1 <- summary(m)$coef["x","Std. Error"]^2 xpred <- mean(yobs)/b1 varxpred <- (varyobs + xpred^2 * varb1) / b1^2 sdxpred <- sqrt(varxpred) } else { b0 <- summary(m)$coef["(Intercept)","Estimate"] - b1 <- summary(m)$coef["xi","Estimate"] + b1 <- summary(m)$coef["x","Estimate"] + S <- summary(m)$sigma xpred <- (mean(yobs) - b0)/b1 sumxxbar <- sum((xi - mean(xi))^2) - if (!syobs) - { - yybar <- (mean(yobs) - mean(yi))^2 - sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5 - } else - { - sdxpred <- ((varyobs^0.5/b1) + (S/b1)^2 * (1/n + ((xpred - mean(xi))^2)/sumxxbar))^0.5 - } + yybar <- (mean(yobs) - mean(yi))^2 + sdxpred <- (S/b1) * (1/p + 1/n + yybar/(b1^2 * sumxxbar))^0.5 } - t <- qt((1 + level)/2,ntot - 2) + t <- qt((1 + level)/2,n - 2) confxpred <- t * sdxpred - result <- c(estimate=xpred,sdxpred=sdxpred,syobs=syobs, - confxpred=confxpred) - digits <- max(c(3,round(log10(xpred/confxpred)) + 2)) - roundedresult <- round(result,digits=digits) - confidenceinterval <- paste(roundedresult["estimate"],"+-", - roundedresult["confxpred"],xunit) - roundedresult[["confidenceinterval"]] <- confidenceinterval - invisible(roundedresult) -} - -multical <- function(cf,df,intercept=FALSE) -{ - rf <- data.frame(name=levels(df$name)) - substances <- colnames(df)[-1] - for (s in substances) - { - r <- vector() - for (sample in levels(df$name)) - { - tmp <- calpredict(subset(df,name==sample)[[s]], - cf[["conc"]],cf[[s]], - intercept=intercept) - r <- c(r,tmp[["confidenceinterval"]]) - } - rf[[s]] <- r - } - return(rf) + result <- c(estimate=xpred,sdxpred=sdxpred,confxpred=confxpred) } diff --git a/data/pahCalibration.rda b/data/pahCalibration.rda deleted file mode 100644 index a0ae738..0000000 Binary files a/data/pahCalibration.rda and /dev/null differ diff --git a/data/pahMeasurements.rda b/data/pahMeasurements.rda deleted file mode 100644 index 657e552..0000000 Binary files a/data/pahMeasurements.rda and /dev/null differ diff --git a/man/calm.Rd b/man/calm.Rd new file mode 100644 index 0000000..c16f663 --- /dev/null +++ b/man/calm.Rd @@ -0,0 +1,43 @@ +\name{calm} +\alias{calm} +\alias{print.calm} +\alias{predict.calm} +\alias{summary.calm} +\title{Generate a calibration model} +\description{ + This function fits a calibration model to the data + frame. +} +\usage{ + calm(data) +} +\arguments{ + \item{data}{ + A data frame with numeric x data in the first column and + numeric y data in the second column. + } +} +\value{ + An object of class \code{calm}, which is derived from + a linear model \code{lm}, the only difference being that + it contains the additional attributes \code{xname}, + \code{yname} and \code{intercept}, the latter being a + boolean reporting wether the model uses an intercept or not. +} +\note{ + The decision if the returned model contains an intercept is taken based on + the significance of the fitted intercept on a significance level of 0.95. + The methods \code{\link{print.calm}}, \code{\link{predict.calm}} + \code{\link{summary.calm}} are just newly assigned names for the + corresponding methods from the class \code{\link{lm}}. +} +\examples{ + data(din32645) + calm(din32645) +} +\author{ + Johannes Ranke + \email{jranke@uni-bremen.de} + \url{http://www.uft.uni-bremen.de/chemie/ranke} +} +\keyword{regression} diff --git a/man/calplot.Rd b/man/calplot.Rd deleted file mode 100644 index 7500912..0000000 --- a/man/calplot.Rd +++ /dev/null @@ -1,52 +0,0 @@ -\name{calplot} -\alias{calplot} -\title{Plot calibration graphs} -\description{ - Produce graphics of calibration data, the fitted model as well - as prediction and confidence intervals. -} -\usage{ - calplot(x,y,intercept=FALSE,measurand="substance x",xunit="mg/L",yunit="Area",level=0.95) -} -\arguments{ - \item{x}{ - A vector of x values. - } - \item{y}{ - A vector of y values. - } - \item{intercept}{ - A boolean describing if the calibration curve is to be forced - through zero. - } - \item{measurand}{ - The name of what is being measured as a character vector. - } - \item{xunit}{ - The unit of the given values on the x axis as a character vector. - } - \item{yunit}{ - The unit of the y axis as a character vector. Defaults to "Area". - } - \item{level}{ - The confidence level of the confidence and prediction bands. Defaults to - 0.95. - } -} -\value{ - A linear model object for y ~ x. You will also get a plot of the calibration - data, of your fitted model as well as lines showing the confidence limits and - the prediction limits. -} -\examples{ -data(pahCalibration) -attach(pahCalibration) -\dontrun{calplot(conc,phenanthrene,"Phenanthrene","mg/L")} -detach(pahCalibration) -} -\author{ - Johannes Ranke - \email{jranke@uni-bremen.de} - \url{http://www.uft.uni-bremen.de/chemie/ranke} -} -\keyword{regression} diff --git a/man/calpredict.Rd b/man/calpredict.Rd deleted file mode 100644 index a91d018..0000000 --- a/man/calpredict.Rd +++ /dev/null @@ -1,64 +0,0 @@ -\name{calpredict} -\alias{calpredict} -\title{Estimate measurement results including confidence intervals} -\description{ - This function generates estimates for x values from y values, including - a confidence interval for the x values. The formulas in this function used - for prediction of concentrations from (replicate) measurements are taken from - the "Handbook of Chemometrics and Qualimetrics Part A" by D. L. Massart, - Vandeginste, B. G. M., Buydens, L. M. C., De Jong, S., Lewi, P. J. and - Smeyers-Verbeke, J, Elsevier, Amsterdam, 1997 and from the EURACHEM/CITAC - report on "Quantifying uncertainty in analytical measurement", 2000, - pp. 111f. -} -\usage{ - calpredict(yobs,xi,yi,xunit="",level=0.95,intercept=FALSE,syobs=FALSE) -} -\arguments{ - \item{yobs}{ - A numeric vector containing the observed data. - } - \item{xi}{ - A vector of x values of the calibration. - } - \item{yi}{ - A vector of y values of the calibration. - } - \item{xunit}{ - The unit of the given values on the x axis as a character string. - } - \item{level}{ - The desired confidence level for the confidence interval of the - estimates. Defaults to 0.95. - } - \item{intercept}{ - Logical value determining if an intercept is to be fitted or not. - Default is FALSE. - } - \item{syobs}{ - If TRUE, a standard deviation for the given y values is - calculated, and the resulting confidence interval will - include this variability (not validated yet). If FALSE (default), this - standard deviation is not included in the - confidence interval. If a numeric value is given, - it is used for the standard deviation of "real samples", - in addition to the standard deviation of the y values - in the calibration (also not validated yet). - } -} -\value{ - A list containing the estimate, its standard deviation and its - confidence interval. -} -\examples{ -data(pahCalibration) -attach(pahCalibration) -y <- c(51.2,51.4,51.1,51.8) -estimate <- calpredict(y,conc,acenaphthene,xunit="mg/L") -} -\author{ - Johannes Ranke - \email{jranke@uni-bremen.de} - \url{http://www.uft.uni-bremen.de/chemie/ranke} -} -\keyword{regression} diff --git a/man/din32645.Rd b/man/din32645.Rd index 901bdf0..0a2a790 100644 --- a/man/din32645.Rd +++ b/man/din32645.Rd @@ -5,7 +5,7 @@ \description{ Sample dataset to test the package. } -\usage{data(pahCalibration)} +\usage{data(din32645)} \format{ A dataframe containing 10 rows of x and y values. } diff --git a/man/multical.Rd b/man/multical.Rd deleted file mode 100644 index 7d62277..0000000 --- a/man/multical.Rd +++ /dev/null @@ -1,37 +0,0 @@ -\name{multical} -\alias{multical} -\title{Calculation of results in a dataframe} -\description{ - This function provides the possibility to perform the calibration - of multiple components simultaneously, and to provide the results - including confidence intervals in a dataframe. -} -\usage{ - multical(cf,df,intercept=FALSE) -} -\arguments{ - \item{cf}{ - A data frame containig the data for the calibration standards. - } - \item{df}{ - A data frame with the measured sample data. - } - \item{intercept}{ - Logical value determining if an intercept is to be fitted or not. - Default is FALSE. - } -} -\value{ - A data frame containig the measurement results with a confidence interval. -} -\examples{ -data(pahCalibration) -data(pahMeasurements) -result <- multical(pahCalibration,pahMeasurements) -} -\author{ - Johannes Ranke - \email{jranke@uni-bremen.de} - \url{http://www.uft.uni-bremen.de/chemie/ranke} -} -\keyword{regression} diff --git a/man/pahCalibration.Rd b/man/pahCalibration.Rd deleted file mode 100644 index 540af4a..0000000 --- a/man/pahCalibration.Rd +++ /dev/null @@ -1,19 +0,0 @@ -\name{pahCalibration} -\docType{data} -\alias{pahCalibration} -\title{Calibration data for HPLC measurement of 4 PAH} -\description{ - This dataset was produced during a course on trace analysis - of organic contaminants. The data are far from perfect, but - good enough to serve as an example. -} -\usage{data(pahCalibration)} -\format{ - A dataframe containing the areas for the four polycyclic aromatic - hydrocarbons (PAH) Acenaphthene, Phenanthrene, Anthracene and Pyrene. - Each measurement of a standard solution makes up one row. -} -\source{ - \url{http://www.uft.uni-bremen.de/chemie} -} -\keyword{datasets} diff --git a/man/pahMeasurements.Rd b/man/pahMeasurements.Rd deleted file mode 100644 index 77b646b..0000000 --- a/man/pahMeasurements.Rd +++ /dev/null @@ -1,19 +0,0 @@ -\name{pahMeasurements} -\docType{data} -\alias{pahMeasurements} -\title{Measurement data for HPLC measurement of 4 PAH} -\description{ - This dataset was produced during a course on trace analysis - of organic contaminants. The data are far from perfect, but - good enough to serve as an example. -} -\usage{data(pahMeasurements)} -\format{ - A dataframe containing the areas for the four polycyclic aromatic - hydrocarbons (PAH) Acenaphthene, Phenanthrene, Anthracene and Pyrene - in the different samples. -} -\source{ - \url{http://www.uft.uni-bremen.de/chemie} -} -\keyword{datasets} diff --git a/man/plot.calm.Rd b/man/plot.calm.Rd new file mode 100644 index 0000000..bb302c7 --- /dev/null +++ b/man/plot.calm.Rd @@ -0,0 +1,48 @@ +\name{plot.calm} +\alias{plot.calm} +\title{Plot calibration graphs from calibration models} +\description{ + Produce graphics of calibration data, the fitted model as well + as prediction and confidence intervals. +} +\usage{ + plot.calm(x,...,xunit="",yunit="",measurand="",level=0.95) +} +\arguments{ + \item{x}{ + A calibration model of type \code{\link{calm}}. It is named + x here because the generic plot method expects x to be its + first argument. + } + \item{...}{ + I just included this because I wanted to avoid the error messages + from R CMD check that tell me I should read "Writing R extensions" + which I did ... + } + \item{xunit}{ + The unit of the given values on the x axis as a character vector. + } + \item{yunit}{ + The unit of the y axis as a character vector. + } + \item{measurand}{ + The name of what is being measured as a character vector. + } + \item{level}{ + The confidence level of the confidence and prediction bands. Defaults to + 0.95. + } +} +\value{ + A plot of the calibration data, of your fitted model as well as lines showing + the confidence limits and the prediction limits. +} +\examples{ + +} +\author{ + Johannes Ranke + \email{jranke@uni-bremen.de} + \url{http://www.uft.uni-bremen.de/chemie/ranke} +} +\keyword{regression} diff --git a/man/predictx.Rd b/man/predictx.Rd new file mode 100644 index 0000000..a3946b0 --- /dev/null +++ b/man/predictx.Rd @@ -0,0 +1,37 @@ +\name{predictx} +\alias{predictx} +\title{Predict x from y values for calibration models} +\description{ + This function predicts x values from y values, as in classical calibration, + including a confindence interval. +} +\usage{ + predictx(m,yobs,level=0.95) +} +\arguments{ + \item{m}{ + A calibration model of type \code{\link{calm}}. + } + \item{yobs}{ + A vector of observed y values for one sample. + } + \item{level}{ + The confidence level for the confidence interval to be reported. + } +} +\value{ + A vector containing the estimate (\code{estimate}), its estimated standard + deviation (\code{sdxpred}), its estimated confidence (\code{confxpred}). +} +\examples{ + data(din32645) + m <- calm(din32645) + r <- predictx(m,3500,level=0.95) + cat("\nThe confidence interval is",r[["estimate"]],"+-",r[["confxpred"]],"\n") +} +\author{ + Johannes Ranke + \email{jranke@uni-bremen.de} + \url{http://www.uft.uni-bremen.de/chemie/ranke} +} +\keyword{regression} -- cgit v1.2.1