Fix inverse predictions for replicate measurements

For details, see NEWS.md

4 files changed, 60 insertions, 21 deletions

diff --git a/tests/testthat/test_compare_investr.R b/tests/testthat/test_compare_investr.R
new file mode 100644
index 0000000..6191c89
--- /dev/null
+++ b/tests/testthat/test_compare_investr.R
@@ -0,0 +1,30 @@
+context("Compare with investr::calibrate")
+
+library(chemCal)
+library(investr)
+
+test_that("Unweighted regressions give same results as investr::calibrate using the Wald method", {
+  compare_investr <- function(object, y_sample) {
+    pred_chemCal <- inverse.predict(object, y_sample)
+    pred_investr <- calibrate(object, y_sample, interval = "Wald")
+    expect_equivalent(pred_chemCal[["Prediction"]], 
+                      pred_investr$estimate)
+    expect_equivalent(pred_chemCal[["Standard Error"]], 
+                      pred_investr$se)
+    expect_equivalent(pred_chemCal[["Confidence Limits"]][1], 
+                      pred_investr$lower)  
+    expect_equivalent(pred_chemCal[["Confidence Limits"]][2],
+                      pred_investr$upper)  
+  }
+  m_tol <- lm(peak_area ~ amount, data = rl95_toluene)
+  compare_investr(m_tol, 1000)
+
+  m_din <- lm(y ~ x, din32645)
+  compare_investr(m_din, 5000)
+
+  m_m1 <- lm(y ~ x, massart97ex1)
+  compare_investr(m_m1, 15)
+
+  m_m3 <- lm(y ~ x, massart97ex3)
+  compare_investr(m_m3, 15)
+})
diff --git a/tests/testthat/test_inverse.predict.R b/tests/testthat/test_inverse.predict.R
index 61484dc..a4e1bfd 100644
--- a/tests/testthat/test_inverse.predict.R
+++ b/tests/testthat/test_inverse.predict.R
@@ -23,20 +23,23 @@ test_that("Inverse predictions for unweighted regressions are stable", {
 })
 
 test_that("Inverse predictions for weighted regressions are stable", {
-  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)]
-  m3 <- lm(y ~ x, w = weights)
-
-  p3.1 <- inverse.predict(m3, 15, ws = 1.67)
+  weights <- with(massart97ex3, {
+    yx <- split(y, x)
+    ybar <- sapply(yx, mean)
+    s <- round(sapply(yx, sd), digits = 2)
+    w <- round(1 / (s^2), digits = 3)
+  })
+
+  massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean)
+
+  m3.means <- lm(y ~ x, w = weights, data = massart97ex3.means)
+
+  p3.1 <- inverse.predict(m3.means, 15, ws = 1.67)
   expect_equal(signif(p3.1$Prediction, 7), 5.865367)
   expect_equal(signif(p3.1$`Standard Error`, 7), 0.8926109)
   expect_equal(signif(p3.1$Confidence, 7), 2.478285)
 
-  p3.2 <- inverse.predict(m3, 90, ws = 0.145)
+  p3.2 <- inverse.predict(m3.means, 90, ws = 0.145)
   expect_equal(signif(p3.2$Prediction, 7), 44.06025)
   expect_equal(signif(p3.2$`Standard Error`, 7), 2.829162)
   expect_equal(signif(p3.2$Confidence, 7), 7.855012)
diff --git a/tests/testthat/test_lod_loq.R b/tests/testthat/test_lod_loq.R
index 6ba0ad0..da0f45e 100644
--- a/tests/testthat/test_lod_loq.R
+++ b/tests/testthat/test_lod_loq.R
@@ -15,7 +15,11 @@ test_that("lod is stable across chemCal versions", {
 })
 
 test_that("loq is stable across chemCal versions", {
-  m2 <- lm(y ~ x, data = massart97ex3)
+  # Actually it was not stable between chemCal <0.2 and 
+  # chemCal > 0.2, so we needed to adapt the test
+  # to work on a model on the means
+  massart97ex3_means <- aggregate(y ~ x, massart97ex3, mean)
+  m2 <- lm(y ~ x, data = massart97ex3_means)
   loq_1 <- loq(m2)
   expect_equal(signif(loq_1$x, 7), 13.97764)
   expect_equal(signif(loq_1$y, 7), 30.6235)
diff --git a/tests/testthat/test_massart.R b/tests/testthat/test_massart.R
index 791c2e7..6e4ce75 100644
--- a/tests/testthat/test_massart.R
+++ b/tests/testthat/test_massart.R
@@ -19,22 +19,24 @@ test_that("Inverse predictions for example 1 are correct",{
   expect_equal(round(p1.3$Confidence, 1), 3.2)
 })
 
+test_that("Inverse predictions for example data 3 are correct when regressing on means",{
+  weights <- with(massart97ex3, {
+    yx <- split(y, x)
+    ybar <- sapply(yx, mean)
+    s <- round(sapply(yx, sd), digits = 2)
+    w <- round(1 / (s^2), digits = 3)
+  })
 
-test_that("Inverse predictions for example 3 are correct",{
-  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)]
-  m3 <- lm(y ~ x, w = weights)
+  massart97ex3.means <- aggregate(y ~ x, massart97ex3, mean)
+
+  m3.means <- lm(y ~ x, w = weights, data = massart97ex3.means)
 
   # Known values are from the book
-  p3.1 <- inverse.predict(m3, 15, ws = 1.67)
+  p3.1 <- inverse.predict(m3.means, 15, ws = 1.67)
   expect_equal(round(p3.1$Prediction, 1), 5.9)
   expect_equal(round(p3.1$Confidence, 1), 2.5)
 
-  p3.2 <- inverse.predict(m3, 90, ws = 0.145)
+  p3.2 <- inverse.predict(m3.means, 90, ws = 0.145)
   expect_equal(round(p3.2$Prediction, 1), 44.1)
   expect_equal(round(p3.2$Confidence, 1), 7.9)
 })