address interaction bug, fixes #35

NEWS.md

@@ -1,3 +1,9 @@
+# fishpond 2.7.1
+
+* Corrected a bug in the two-group interaction (without pairing)
+  functionality, when the groups were imbalanced, as identified by
+  Samuel Chen. Fixes GitHub issue #35.
+
 # fishpond 2.5.4
 
 * As CellRanger 7 includes both spliced and unspliced counts in their

R/interaction.R

@@ -52,7 +52,7 @@ swishInterx <- function(infRepsArray, condition, covariate,
   if (!all(table(condition, covariate) > 0))
     stop("swish with interaction across two variables requires samples for each combination")
   dims <- dim(infRepsArray)
-
+  
   tab <- table(condition, covariate)
   # if sizes are equal, don't need to double or splice out columns
   all_equal <- all(tab[,2] == tab[,1])
@@ -158,15 +158,16 @@ randomSamplesToRemove <- function(tab, condition, covariate) {
   cond2 <- condition == levels(condition)[2]
   cov_lvls <- levels(covariate)
   idx <- numeric()
-  for (i in which(tab[,1] != tab[,2])) {
-    cond1small <- tab[1,i] < tab[2,i]
+  # find covariate levels which are imbalanced
+  for (j in which(tab[1,] != tab[2,])) {
+    cond1small <- tab[1,j] < tab[2,j]
     if (cond1small) {
-      idx <- c(idx, sample(which(cond2 & covariate == cov_lvls[i]),
-                           tab[2,i] - tab[1,i],
+      idx <- c(idx, sample(which(cond2 & covariate == cov_lvls[j]),
+                           tab[2,j] - tab[1,j],
                            replace=FALSE))
     } else {
-      idx <- c(idx, sample(which(cond1 & covariate == cov_lvls[i]),
-                           tab[1,i] - tab[2,i],
+      idx <- c(idx, sample(which(cond1 & covariate == cov_lvls[j]),
+                           tab[1,j] - tab[2,j],
                            replace=FALSE))
     }
   }

tests/testthat/test_interaction.R

@@ -49,7 +49,7 @@ test_that("two group interactions work", {
 
   expect_true(mcols(y)$pvalue[2] < .01)
 
-  # two groups with unbalanced sample sizes
+  # two groups with imbalanced sample sizes
   y <- makeSimSwishData(m=200, n=20, null=TRUE)
   nms <- c("counts",paste0("infRep",1:20))
   lambda1 <- rep(c(40,80,40,80),c(4,6,6,4))
@@ -60,6 +60,7 @@ test_that("two group interactions work", {
   }
   y$condition <- factor(rep(c(1,2,1,2),c(4,6,6,4)))
   y$group <- factor(rep(1:2,each=10))
+  table(y$condition, y$group)
 
   y <- scaleInfReps(y, quiet=TRUE)
   y <- labelKeep(y)