1 files changed, 72 insertions, 23 deletions
diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R
index f438d99..f1921e4 100644
--- a/R/inverse.predict.lm.R
+++ b/R/inverse.predict.lm.R
@@ -1,40 +1,89 @@
 # This is an implementation of Equation (8.28) in the Handbook of Chemometrics
-# and Qualimetrics, Part A, Massart et al, page 200, validated with Example 8
-# on the same page
+# and Qualimetrics, Part A, Massart et al (1997), page 200, validated with
+# Example 8 on the same page
 
-inverse.predict <- function(object, newdata, alpha=0.05)
+inverse.predict <- function(object, newdata, 
+  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
+  alpha=0.05, ss = "auto")
 {
   UseMethod("inverse.predict")
 }
 
-inverse.predict.default <- function(object, newdata, alpha=0.05)
+inverse.predict.default <- function(object, newdata, 
+  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
+  alpha=0.05, ss = "auto")
 {
   stop("Inverse prediction only implemented for univariate lm objects.")
 }
 
-inverse.predict.lm <- function(object, newdata, alpha=0.05)
+inverse.predict.lm <- function(object, newdata, 
+  ws = ifelse(length(object$weights) > 0, mean(object$weights), 1),
+  alpha=0.05, ss = "auto")
 {
-  if (is.list(newdata)) {
-    if (!is.null(newdata$y))
-      newdata <- newdata$y
-    else
-      stop("Newdata list should contain element newdata$y")
-  }
+  if (length(object$coef) > 2)
+    stop("More than one independent variable in your model - not implemented")
+
+  if (alpha <= 0 | alpha >= 1)
+    stop("Alpha should be between 0 and 1 (exclusive)")
+
+  ybars <- mean(newdata)
+  m <- length(newdata)
 
-  if (is.matrix(newdata)) {
-    Y <- newdata
-    Ybar <- apply(Y,1,mean)
-    nrepl <- ncol(Y)
+  yx <- split(object$model$y,object$model$x)
+  n <- length(yx)
+  x <- as.numeric(names(yx))
+  ybar <- sapply(yx,mean)
+  if (length(object$weights) > 0) {
+    wx <- split(object$weights,object$model$x)
+    w <- sapply(wx,mean)
+  } else {
+    w <- rep(1,n)
   }
-  else {
-    Y <- as.vector(newdata)
-    Ybar <- Y
-    nrepl <- 1 
+  yhatx <- split(object$fitted.values,object$model$x)
+  yhat <- sapply(yhatx,mean)
+  se <- sqrt(sum(w*(ybar - yhat)^2)/(n-2))
+  if (ss == "auto") {
+    ss <- se
+  } else {
+    ss <- ss
   }
 
-  if (length(object$coef) > 2)
-    stop("Inverse prediction not yet implemented for more than one independent variable")
+  b1 <- object$coef[["x"]]
 
-  if (alpha <= 0 | alpha >= 1)
-    stop("Alpha should be between 0 and 1 (exclusive)")
+  ybarw <- sum(w * ybar)/sum(w)
+
+# The commented out code for sxhats is equation 8.28 without changes.  It has
+# been replaced by the code below, in order to be able to take into account a
+# precision in the sample measurements that differs from the precision in the
+# calibration samples.
+
+#  sxhats <- se/b1 * sqrt(
+#    1/(ws * m) + 
+#    1/sum(w) + 
+#    (ybars - ybarw)^2 * sum(w) /
+#      (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2))
+#  )
+
+# This is equation 8.28, but with the possibility to take into account a 
+# different precision measurement of the sample and standard solutions
+# in analogy to equation 8.26
+  sxhats <- 1/b1 * sqrt(
+    ss^2/(ws * m) + 
+    se^2 * (1/sum(w) + 
+      (ybars - ybarw)^2 * sum(w) /
+        (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2)))
+  )
+
+  if (names(object$coef)[1] == "(Intercept)") {
+    b0 <- object$coef[["(Intercept)"]]
+  } else {
+    b0 <- 0
+  }
+
+  xs <- (ybars - b0) / b1
+  t <- qt(1-0.5*alpha, n - 2)
+  conf <- t * sxhats
+  result <- list("Prediction"=xs,"Standard Error"=sxhats,
+    "Confidence"=conf, "Confidence Limits"=c(xs - conf, xs + conf))
+  return(result)
 }