Fix for Issue #196

mixOmicsTeam · Aug 9, 2022 · b4791fd · b4791fd
1 parent 0ff2d16
commit b4791fd
Show file tree

Hide file tree

Showing 2 changed files with 59 additions and 7 deletions.
diff --git a/R/perf.R b/R/perf.R
@@ -553,8 +553,8 @@ perf.mixo_spls  <- perf.mixo_pls
             if(sum(is.na(Y.hat))>0) stop('Predicted Y values include NA')  
 
             # replaced h by 1; Y.hat is the prediction of the test samples for all q variable in comp h = 1
-            Ypred[omit, , h] = Y.hat[, , 1]
-            MSEP.mat[omit, , h] = (Y.test - Y.hat[, , 1])^2
+            Ypred[omit, nzv.Y, h] = Y.hat[, , 1]
+            MSEP.mat[omit, nzv.Y, h] = (Y.test[, nzv.Y] - Y.hat[, , 1])^2
 
 
             # Q2 criterion: buidling directly from spls object
@@ -573,7 +573,7 @@ perf.mixo_spls  <- perf.mixo_pls
             t.pred.cv[omit,h] = t.pred    # needed for tuning
             b.pred = crossprod(Y.test, t.pred)
             b.pred.cv = b.pred/ drop(sqrt(crossprod(b.pred)))
-            u.pred.cv[omit,h] = Y.test %*% b.cv  # needed for tuning, changed instead of b.pred.cv
+            u.pred.cv[omit,h] = Y.test[, nzv.Y] %*% b.cv  # needed for tuning, changed instead of b.pred.cv
 
             # predicted reg coeff, could be removed
             e.pred.cv = crossprod(as.matrix(Y.test), Y.test %*% b.pred.cv) / drop(crossprod(Y.test %*% b.pred))
@@ -585,19 +585,19 @@ perf.mixo_spls  <- perf.mixo_pls
             # deflate matrices X      
             #-- mode classic
             if (mode == "classic"){
-                Y.train = Y.train - t.cv %*% t(b.cv)  # could be pred on b
-                Y.test = Y.test - t.pred %*% t(b.cv)
+                Y.train[, nzv.Y] = Y.train[, nzv.Y] - t.cv %*% t(b.cv)  # could be pred on b
+                Y.test[, nzv.Y] = Y.test[, nzv.Y] - t.pred %*% t(b.cv)
             }
             #-- mode regression
             if (mode == "regression"){
                 Y.train = Y.train - t.cv %*% t(d.cv) # could be pred d.pred.cv? does not decrease enough
-                Y.test = Y.test - Y.hat[, , 1]   # == Y.test - t.pred %*% t(d.cv) 
+                Y.test[, nzv.Y] = Y.test[, nzv.Y] - Y.hat[, , 1]   # == Y.test - t.pred %*% t(d.cv) 
             }
 
             #-- mode canonical  ## KA added
             if (mode == "canonical"){
                 Y.train = Y.train - u.cv %*% t(e.cv)  # could be pred on e
-                Y.test = Y.test - (Y.test %*% b.cv) %*% t(e.cv)  # here u.pred = Y.test %*% b.cv (b.pred.cv gives similar results)
+                Y.test = Y.test - (Y.test[, nzv.Y] %*% b.cv) %*% t(e.cv)  # here u.pred = Y.test %*% b.cv (b.pred.cv gives similar results)
             }
             #-- mode invariant: Y is unchanged
 

diff --git a/tests/testthat/test-perf.pls.R b/tests/testthat/test-perf.pls.R
@@ -0,0 +1,52 @@
+context("perf.pls")
+
+test_that("perf() works on pls object", code = {
+  library(mixOmics)
+
+  data("liver.toxicity")
+
+  # reducing number of features to reduce run time
+  X <- liver.toxicity$gene[1:500]
+  Y <- liver.toxicity$clinic
+
+  set.seed(12)
+  pls.obg <- pls(Y, X, ncomp = 4)
+  pls.perf.obj <- perf(pls.obg, validation = "Mfold", folds = 4, 
+                       progressBar = F, 
+                       nrepeat = 3)
+
+  trueVals <- c(-0.017, -0.294, -0.431, -0.622)
+  testVals <- round(pls.perf.obj$measures$Q2.total$summary[, "mean"], 3)
+
+  expect_equal(trueVals, testVals)
+})
+
+test_that("perf() works on pls with nzv features (all modes)", code = {
+  library(mixOmics)
+
+  data("liver.toxicity")
+
+  # reducing number of features to reduce run time
+  X <- liver.toxicity$gene[, 1:1000]
+  Y <- liver.toxicity$clinic
+
+  # to reproduce error, we need to induce some features to have near zero variance
+  X[, c(1, 23, 62, 234, 789)] <- 0
+
+  modes <- c("classic", "regression", "canonical")
+  trueVals <- list(c(0.031,  0.007,  0.003, -0.006),
+                   c(0.006, -0.222, -0.379, -0.553),
+                   c(0.089, -0.473, -1.238, -2.228))
+
+  for (m in 1:3) {
+    set.seed(12)
+    suppressWarnings(pls.obg <- pls(Y, X, ncomp = 4, mode = modes[m]))
+    suppressWarnings(pls.perf.obj <- perf(pls.obg, validation = "Mfold", folds = 4, 
+                         progressBar = F, 
+                         nrepeat = 3))
+
+    testVals <- round(pls.perf.obj$measures$Q2.total$summary[, "mean"], 3)
+    expect_equal(trueVals[[m]], testVals)
+  }
+
+})