- Fixed the inverse prediction

- Now I have a working approach for the calculation of LOD and LOQ,
  but it seems to be different from what everybody else is doing
  (e.g. Massart chaper 13). I like it, however. Maybe it even 
  yields a paper.



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

5 files changed, 104 insertions, 7 deletions

diff --git a/man/calplot.lm.Rd b/man/calplot.lm.Rd
index efdea3d..6d3f52d 100644
--- a/man/calplot.lm.Rd
+++ b/man/calplot.lm.Rd
@@ -30,7 +30,7 @@
     The label of the y axis.
   }
   \item{alpha}{
-    The confidence level for the confidence and prediction bands.
+    The error tolerance level for the confidence and prediction bands.
   }
 }
 \value{
diff --git a/man/din32645.Rd b/man/din32645.Rd
index 1a3c046..d251b7c 100644
--- a/man/din32645.Rd
+++ b/man/din32645.Rd
@@ -13,9 +13,19 @@
 data(din32645)
 m <- lm(y ~ x, data=din32645)
 calplot(m)
-prediction <- inverse.predict(m,3500,alpha=0.01)
+(prediction <- inverse.predict(m,3500,alpha=0.01))
 # This should give 0.074 according to DIN (cited from the Dintest test data)
 round(prediction$Confidence,3)
+
+# According to Dintest, we should get 0.07, but we get 0.0759
+lod(m, alpha = 0.01)
+
+# In German, there is the "Erfassungsgrenze", with k = 2, 
+# and we should get 0.14 according to Dintest
+lod(m, k = 2, alpha = 0.01)
+
+# According to Dintest, we should get 0.21, we get 0.212
+loq(m, alpha = 0.01)
 }
 \references{
   DIN 32645 (equivalent to ISO 11843)
diff --git a/man/inverse.predict.Rd b/man/inverse.predict.Rd
index d773e58..8c2be9c 100644
--- a/man/inverse.predict.Rd
+++ b/man/inverse.predict.Rd
@@ -4,9 +4,8 @@
 \alias{inverse.predict.rlm}
 \alias{inverse.predict.default}
 \title{Predict x from y for a linear calibration}
-\usage{inverse.predict(object, newdata, 
-  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
-  alpha=0.05, ss = "auto")
+\usage{inverse.predict(object, newdata, \dots,
+  ws, alpha=0.05, ss = "auto")
 }
 \arguments{
   \item{object}{
@@ -17,12 +16,16 @@
   \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. The default is to take the 
     mean of the weights in the model, if there are any.
   }
   \item{alpha}{
-    The confidence level for the confidence interval to be reported.
+    The error tolerance level for the confidence interval to be reported.
   }
   \item{ss}{
     The estimated standard error of the sample measurements. The 
@@ -62,6 +65,6 @@ w <- round(1/(s^2),digits=3)
 weights <- w[factor(x)]
 m <- lm(y ~ x,w=weights)
 
-inverse.predict(m,c(15))
+inverse.predict(m,15,ws = 1.67)
 }
 \keyword{manip}
diff --git a/man/lod.Rd b/man/lod.Rd
new file mode 100644
index 0000000..e6ce345
--- /dev/null
+++ b/man/lod.Rd
@@ -0,0 +1,58 @@
+\name{lod}
+\alias{lod}
+\alias{lod.lm}
+\alias{lod.rlm}
+\alias{lod.default}
+\alias{loq}
+\alias{loq.lm}
+\alias{loq.rlm}
+\alias{loq.default}
+\title{Estimate a limit of detection (LOD) or quantification (LOQ)}
+\usage{
+  lod(object, \dots, alpha = 0.05, k = 1, n = 1, w = "auto")
+  loq(object, \dots, alpha = 0.05, k = 3, n = 1, w = "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.
+  }
+  \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 LOD/LOQ.
+  }
+  \item{n}{
+    The number of replicate measurements for which the LOD/LOQ should be
+    specified.
+  }
+  \item{w}{
+    The weight that should be attributed to the LOD/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.
+  }
+}
+\value{
+  The estimated limit of detection for a model used for calibration.
+}
+\description{
+  A useful operationalisation of a lower limit L of a measurement method is
+  simply the solution of the equation
+    \deqn{L = k c(L)}{L = k * c(L)}
+  where c(L) is half of the lenght of the confidence interval at the limit L.
+}
+\examples{
+  data(din32645)
+  m <- lm(y ~ x, data = din32645)
+  lod(m)
+}
+\keyword{manip}
diff --git a/man/massart97ex3.Rd b/man/massart97ex3.Rd
index 1dabf90..2618709 100644
--- a/man/massart97ex3.Rd
+++ b/man/massart97ex3.Rd
@@ -10,6 +10,32 @@
   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)
+# The following concords with the book
+inverse.predict(m, 15, ws = 1.67)
+inverse.predict(m, 90, ws = 0.145)
+
+calplot(m)
+
+m0 <- lm(y ~ x)
+lod(m0)
+lod(m)
+
+# 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)
+}
 \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