Substantial performance improvements from comments by Eytan Bakshy. U…

…sing foreach's combine option instead of apply and rbind.

multiway_boot.R

@@ -1,30 +1,31 @@
 library(plyr)
+library(foreach)
 
 r.double.or.nothing <- function(n) {
   2 * rbinom(n, 1, .5)
 }
 
 multiway.boot <- function(
-  statistic, R, N,
+  statistic, R,
   groups = as.matrix(1:N),
   verbose = FALSE,
   RNG = r.double.or.nothing,
   .parallel = FALSE,
   .progress = 'none',
   ...
   ) {
-  groups.num <- apply(groups, 2, function(x) as.numeric(as.factor(x)))
-  N.groups <- apply(groups, 2, function(x) length(unique(x)))
+  groups <- apply(groups, 2, function(x) as.numeric(as.factor(x)))
+  N.groups <- apply(groups, 2, function(x) max(x))
   N.groupingFactors <- ncol(groups)
+
   llply(
     1:R,
-    function(i) {
-      W <- matrix(nrow = 0, ncol = N)
-      for (j in 1:N.groupingFactors) {
-        W <- rbind(W, RNG(N.groups[j])[groups.num[,j]])
-      }
+    function(r) {
       # Observation weights are products of weights for each factor
-      w <- apply(W, 2, prod)
+      w <- foreach(i = 1:N.groupingFactors, .combine = `*`) %do% {
+        RNG(N.groups[i])[groups[, i]]
+      }
+
       if (verbose) cat(i, " ")
       statistic(..., weights = w)
     },

multiway_boot_ex.R

@@ -24,13 +24,12 @@ q <- c(alpha / 2, 1 - alpha / 2)
 mb.2 <- multiway.boot(
   statistic = wtd.mean,
   R = 500,
-  N = N,
   groups = cbind(subject, stimulus),
   .progress = 'text',
   x = x
   )
 mb.2 <- unlist(mb.2)
-var(mb.2) # bootstrap estimate of variance of the mean
+sd(mb.2) # bootstrap estimate of standard error of the mean
 qnorm(q, mean(x), sd(mb.2)) # normal approximation
 quantile(mb.2, q) # percentile bootstrap CI
 
@@ -44,13 +43,12 @@ registerDoMC()
 mb.1 <- multiway.boot(
   statistic = wtd.mean,
   R = 500,
-  N = N,
   groups = cbind(subject), # only cluster/block on subject
   .parallel = TRUE, # requires foreach and doMC libraries
   x = x
   )
 mb.1 <- unlist(mb.1)
-var(mb.1)
+sd(mb.1)
 qnorm(q, mean(x), sd(mb.1))
 quantile(mb.1, q)