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

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
diff --git a/data/pahMeasurements.rda b/data/pahMeasurements.rda
deleted file mode 100644
index 657e552..0000000
--- a/data/pahMeasurements.rda
+++ /dev/null
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}