- inverse.predict now has a var.s argument instead of the never

  tested ss argument. This is documented in the updated vignette
- loq() now has w.loq and var.loq arguments, and stops with a message
  if neither are specified and the model has weights.
- calplot doesn't stop any more for weighted regression models, but
  only refrains from drawing prediction bands
- Added method = "din" to lod(), now that I actually have it (DIN 32645) and
  was able to see which approximation is used therein.
- removed the demos, as the examples and tests are already partially 
  duplicated
- The vignette is more of a collection of various notes, but should 
  certainly be helpful for the user.
- Version bump to 0.1-xxx



git-svn-id: http://kriemhild.uft.uni-bremen.de/svn/chemCal@16 5fad18fb-23f0-0310-ab10-e59a3bee62b4

6 files changed, 118 insertions, 101 deletions

diff --git a/man/calplot.lm.Rd b/man/calplot.lm.Rd
index de63022..bf3f616 100644
--- a/man/calplot.lm.Rd
+++ b/man/calplot.lm.Rd
@@ -5,12 +5,11 @@
 \title{Plot calibration graphs from univariate linear models}
 \description{
 	Produce graphics of calibration data, the fitted model as well
-  as prediction and confidence bands. 
+  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)
+  calplot(object, xlim = c("auto", "auto"), ylim = c("auto", "auto"),
+    xlab = "Concentration", ylab = "Response", alpha=0.05, varfunc = NULL)
 }
 \arguments{
   \item{object}{
@@ -40,7 +39,8 @@
 }
 \value{
   A plot of the calibration data, of your fitted model as well as lines showing
-  the confidence limits as well as the prediction limits.
+  the confidence limits. Prediction limits are only shown for models from
+  unweighted regression.
 } 
 \note{
   Prediction bands for models from weighted linear regression require weights
@@ -51,7 +51,7 @@
 }
 \examples{
 data(massart97ex3)
-m <- lm(y ~ x, data=massart97ex3)
+m <- lm(y ~ x, data = massart97ex3)
 calplot(m)
 }
 \author{
diff --git a/man/din32645.Rd b/man/din32645.Rd
index 94486c4..cacbf07 100644
--- a/man/din32645.Rd
+++ b/man/din32645.Rd
@@ -11,39 +11,44 @@
 }
 \examples{
 data(din32645)
-m <- lm(y ~ x, data=din32645)
+m <- lm(y ~ x, data = din32645)
 calplot(m)
-(prediction <- inverse.predict(m,3500,alpha=0.01))
-# This should give 0.07434 according to Dintest test data, as 
-# collected from Procontrol 3.1 (isomehr GmbH)
+
+## 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)
 
-# According to Dintest test data, we should get 0.0698 for the critical value
+## 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")
-(lod <- lod(m, alpha = 0.01, beta = 0.5))
-round(lod$x,4)
+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.140 according to Dintest test data, but with chemCal, we can't
-# reproduce this, 
-lod(m, alpha = 0.01, beta = 0.01) 
-# except by using an equivalent to the approximation
-# xD = 2 * Sc / A (Currie 1999, p. 118, or Orange Book, Chapter 18.4.3.7)
-lod.approx <- 2 * lod$x
-round(lod.approx, digits=3)
-# which seems to be the pragmatic definition in DIN 32645, as judging from
-# the Dintest test data. 
-
-# This accords to the test data from Dintest again, except for the last digit
-# of the value cited for Procontrol 3.1 (0.2121)
+# 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)
+  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. 
diff --git a/man/inverse.predict.Rd b/man/inverse.predict.Rd
index 5be0250..925f3e9 100644
--- a/man/inverse.predict.Rd
+++ b/man/inverse.predict.Rd
@@ -5,7 +5,7 @@
 \alias{inverse.predict.default}
 \title{Predict x from y for a linear calibration}
 \usage{inverse.predict(object, newdata, \dots,
-  ws, alpha=0.05, ss = "auto")
+  ws, alpha=0.05, var.s = "auto")
 }
 \arguments{
   \item{object}{
@@ -21,15 +21,17 @@
     future implementations.
   }
   \item{ws}{ 
-    The weight attributed to the sample. The default is to take the 
-    mean of the weights in the model, if there are any.
+    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{ss}{
-    The estimated standard error of the sample measurements. The 
-    default is to take the residual standard error from the calibration.
+  \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{
@@ -59,12 +61,13 @@
 \examples{
   data(massart97ex3)
   attach(massart97ex3)
-  yx <- split(y,factor(x))
-  s <- round(sapply(yx,sd),digits=2)
-  w <- round(1/(s^2),digits=3)
+  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)
+  m <- lm(y ~ x, w = weights)
 
-  inverse.predict(m,15,ws = 1.67)
+  inverse.predict(m, 15, ws = 1.67) # 5.9 +- 2.5
 }
 \keyword{manip}
diff --git a/man/lod.Rd b/man/lod.Rd
index c2fe4e9..055074c 100644
--- a/man/lod.Rd
+++ b/man/lod.Rd
@@ -5,7 +5,7 @@
 \alias{lod.default}
 \title{Estimate a limit of detection (LOD)}
 \usage{
-  lod(object, \dots, alpha = 0.05, beta = 0.05)
+  lod(object, \dots, alpha = 0.05, beta = 0.05, method = "default")
 }
 \arguments{
   \item{object}{
@@ -24,10 +24,18 @@
   \item{beta}{
     The error tolerance beta for the detection limit.
   }
+  \item{method}{
+    The 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.
+  }
 }
 \value{
   A list containig the corresponding x and y values of the estimated limit of
-  detection of a model used for calibration.  }
+  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
@@ -39,7 +47,7 @@
   one-sided significance test).
 }
 \note{
-  - The default values for alpha and beta are recommended by IUPAC.
+  - 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).
@@ -68,8 +76,8 @@
   # The critical value (decision limit, German Nachweisgrenze) can be obtained
   # by using beta = 0.5:
   lod(m, alpha = 0.01, beta = 0.5)
-  # or approximated by 
-  2 * lod(m, alpha = 0.01, beta = 0.5)$x
-  # for the case of known, constant variance (homoscedastic data)
+}
+\seealso{
+  Examples for \code{\link{din32645}}  
 }
 \keyword{manip}
diff --git a/man/loq.Rd b/man/loq.Rd
index 4850487..0250098 100644
--- a/man/loq.Rd
+++ b/man/loq.Rd
@@ -5,14 +5,17 @@
 \alias{loq.default}
 \title{Estimate a limit of quantification (LOQ)}
 \usage{
-  loq(object, \dots, alpha = 0.05, k = 3, n = 1, w = "auto")
+  loq(object, \dots, alpha = 0.05, k = 3, n = 1, w.loq = "auto",
+    var.loq = "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}, 
-    optionally from a weighted regression.
+    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.
@@ -29,53 +32,46 @@
     The number of replicate measurements for which the LOQ should be
     specified.
   }
-  \item{w}{
+  \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.
+  }
 }
 \value{
   The estimated limit of quantification for a model used for calibration.
 }
 \description{
-  A useful operationalisation of a limit of quantification is simply the
-  solution of the equation
+  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 as
-  estimated by \code{\link{inverse.predict}}. By virtue of this formula, the 
-  limit of detection is the x value, where the relative error
-  of the quantification with the given calibration model is 1/k.
+  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 approach taken here is to my knowledge
-  original to the chemCal package.
+  - 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)
-  m0 <- lm(y ~ x)
-  loq(m0)
-
-  # Now we use a weighted regression
-  yx <- split(y,factor(x))
-  s <- round(sapply(yx,sd),digits=2)
-  w <- round(1/(s^2),digits=3)
-  weights <- w[factor(x)]
-  mw <- lm(y ~ x,w=weights)
-  loq(mw)
-
-  # In order to define the weight at the loq, we can use
-  # the variance function 1/y for the model
-  mwy <- lm(y ~ x, w = 1/y)
+  m <- lm(y ~ x)
+  loq(m)
 
-  # Let's do this with one iteration only
-  loq(mwy, w = 1 / predict(mwy,list(x = loq(mwy)$x)))
-
-  # We can get better by doing replicate measurements
-  loq(mwy, n = 3, w = 1 / predict(mwy,list(x = loq(mwy)$x)))
+  # We can get better by using replicate measurements
+  loq(m, n = 3)
+}
+\seealso{
+  Examples for \code{\link{din32645}}  
 }
 \keyword{manip}
diff --git a/man/massart97ex3.Rd b/man/massart97ex3.Rd
index eb00e79..e7cd383 100644
--- a/man/massart97ex3.Rd
+++ b/man/massart97ex3.Rd
@@ -3,7 +3,7 @@
 \alias{massart97ex3}
 \title{Calibration data from Massart et al. (1997), example 3}
 \description{
-  Sample dataset to test the package.
+  Sample dataset from p. 188 to test the package.
 }
 \usage{data(massart97ex3)}
 \format{
@@ -11,37 +11,42 @@
   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)
-# The following concords with the book
-inverse.predict(m, 15, ws = 1.67)
-inverse.predict(m, 90, ws = 0.145)
+  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)
 
-# Some of the following examples are commented out, because the require
-# prediction intervals from predict.lm for weighted models, which is not
-# available in R at the moment.
+  # 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
 
-#calplot(m)
+  # The LOD is only calculated for models from unweighted regression
+  # with this version of chemCal
+  m0 <- lm(y ~ x) 
+  lod(m0)
 
-m0 <- lm(y ~ x)
-lod(m0)
-#lod(m)
+  # Limit of quantification from unweighted regression
+  m0 <- lm(y ~ x)
+  loq(m0)
 
-# Now we want to take advantage of the lower weights at lower y values
-#m2 <- lm(y ~ x, w = 1/y)
-# To get a reasonable weight for the lod, we need to estimate it and predict 
-# a y value for it
-#yhat.lod <- predict(m,data.frame(x = lod(m2)))
-#lod(m2,w=1/yhat.lod,k=3)
+  # 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, p. 188
+  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}