aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--DESCRIPTION2
-rw-r--r--R/chemCal.R132
-rw-r--r--data/pahCalibration.rdabin1289 -> 0 bytes
-rw-r--r--data/pahMeasurements.rdabin1348 -> 0 bytes
-rw-r--r--man/calm.Rd43
-rw-r--r--man/calplot.Rd52
-rw-r--r--man/calpredict.Rd64
-rw-r--r--man/din32645.Rd2
-rw-r--r--man/multical.Rd37
-rw-r--r--man/pahCalibration.Rd19
-rw-r--r--man/pahMeasurements.Rd19
-rw-r--r--man/plot.calm.Rd48
-rw-r--r--man/predictx.Rd37
13 files changed, 146 insertions, 309 deletions
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 <jranke@uni-bremen.de>
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
--- a/data/pahCalibration.rda
+++ /dev/null
Binary files differ
diff --git a/data/pahMeasurements.rda b/data/pahMeasurements.rda
deleted file mode 100644
index 657e552..0000000
--- a/data/pahMeasurements.rda
+++ /dev/null
Binary files 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}

Contact - Imprint