From e636c17f0d354a8e74546fc1469431dbe502dc76 Mon Sep 17 00:00:00 2001
From: Johannes Ranke <jranke@uni-bremen.de>
Date: Mon, 16 Jul 2018 17:17:26 +0200
Subject: Attempt to fix the problem discovered by Anna Burniol Figols

but then the tests fail...
---
 R/inverse.predict.lm.R | 25 +++++++++++--------------
 man/inverse.predict.Rd |  2 +-
 vignettes/chemCal.Rmd  |  6 +++---
 3 files changed, 15 insertions(+), 18 deletions(-)

diff --git a/R/inverse.predict.lm.R b/R/inverse.predict.lm.R
index 927e672..77d548f 100644
--- a/R/inverse.predict.lm.R
+++ b/R/inverse.predict.lm.R
@@ -23,10 +23,9 @@ inverse.predict.lm <- function(object, newdata, ...,
     ws <- ifelse(length(object$weights) > 0, mean(object$weights), 1)
   }
   if (length(object$weights) > 0) {
-    wx <- split(object$weights,object$model[[xname]])
-    w <- sapply(wx,mean)
+    w <- object$weights
   } else {
-    w <- rep(1,length(split(object$model[[yname]],object$model[[xname]])))
+    w <- rep(1, nrow(object$model))
   }
   .inverse.predict(object = object, newdata = newdata, 
     ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname)
@@ -40,8 +39,7 @@ inverse.predict.rlm <- function(object, newdata, ...,
   if (ws == "auto") {
     ws <- mean(object$w)
   }
-  wx <- split(object$weights,object$model[[xname]])
-  w <- sapply(wx,mean)
+  w <- object$w
   .inverse.predict(object = object, newdata = newdata, 
     ws = ws, alpha = alpha, var.s = var.s, w = w, xname = xname, yname = yname)
 }
@@ -57,14 +55,13 @@ inverse.predict.rlm <- function(object, newdata, ...,
   ybars <- mean(newdata)
   m <- length(newdata)
 
-  yx <- split(object$model[[yname]], object$model[[xname]])
-  n <- length(yx)
+  n <- nrow(object$model)
   df <- n - length(object$coef)
-  x <- as.numeric(names(yx))
-  ybar <- sapply(yx,mean)
-  yhatx <- split(object$fitted.values, object$model[[xname]])
-  yhat <- sapply(yhatx, mean)
-  se <- sqrt(sum(w * (ybar - yhat)^2)/df)
+  xi <- object$model[[2]]
+  yi <- object$model[[1]]
+  yihat <- object$fitted.values
+
+  se <- sqrt(sum(w * (yi - yihat)^2)/df)
 
   if (var.s == "auto") {
     var.s <- se^2/ws
@@ -72,14 +69,14 @@ inverse.predict.rlm <- function(object, newdata, ...,
 
   b1 <- object$coef[[xname]]
 
-  ybarw <- sum(w * ybar)/sum(w)
+  ybarw <- sum(w * yi)/sum(w)
 
 # This is the adapted form of equation 8.28 (see package vignette)
   sxhats <- 1/b1 * sqrt(
     (var.s / m) + 
     se^2 * (1/sum(w) + 
       (ybars - ybarw)^2 * sum(w) /
-        (b1^2 * (sum(w) * sum(w * x^2) - sum(w * x)^2)))
+        (b1^2 * (sum(w) * sum(w * xi^2) - sum(w * xi)^2)))
   )
 
   if (names(object$coef)[1] == "(Intercept)") {
diff --git a/man/inverse.predict.Rd b/man/inverse.predict.Rd
index 347d670..26ba6b8 100644
--- a/man/inverse.predict.Rd
+++ b/man/inverse.predict.Rd
@@ -51,7 +51,7 @@
   This is elaborated in the package vignette.
 }
 \note{
-  The function was validated with examples 7 and 8 from Massart et al. (1997).
+  The function was validated with examples 7 and 8 from Massart et al. (1997). 
 }
 \references{
   Massart, L.M, Vandenginste, B.G.M., Buydens, L.M.C., De Jong, S., Lewi, P.J.,
diff --git a/vignettes/chemCal.Rmd b/vignettes/chemCal.Rmd
index cea1678..2515abb 100644
--- a/vignettes/chemCal.Rmd
+++ b/vignettes/chemCal.Rmd
@@ -107,9 +107,9 @@ with
 s_e = \sqrt{ \frac{\sum w_i (y_i - \hat{y_i})^2}{n - 2}}
 \end{equation}
 
-where $w_i$ is the weight for calibration standard $i$, $y_i$ is the mean $y$
-value (!) observed for standard $i$, $\hat{y_i}$ is the estimated value for
-standard $i$, $n$ is the number calibration standards, $w_s$ is the weight
+where $w_i$ is the weight for calibration standard $i$, $y_i$ is the $y$
+value observed for standard $i$, $\hat{y_i}$ is the estimated value for
+standard $i$, $n$ is the number of calibration samples, $w_s$ is the weight
 attributed to the sample $s$, $m$ is the number of replicate measurements of
 sample $s$, $\bar{y_s}$ is the mean response for the sample,
 $\bar{y_w} = \frac{\sum{w_i y_i}}{\sum{w_i}}$ is the weighted mean of responses
-- 
cgit v1.2.1