Use of permute package to calculate all permutations in a blocked design

Question

I want to calculate all permutations of a blocked design suitable for a Friedman test. Consider the following example:

thedata <- data.frame(
  score = c(replicate(4,sample(1:3))),
  judge = rep(1:4,each=3),
  wine  = rep.int(1:3,4)  
  )

Four judges ranked 3 wines and now I want to calculate every possible permutation within the data for every judge. I expect to see 1,296 permutations, as also given by:

require(permute)
CTRL <- how(within=Within("free"),
            plots=Plots(strata=factor(thedata$judge)),
            complete=TRUE,maxperm=1e9)
numPerms(12,CTRL)

However, allPerms(12,control=CTRL) produces the following error:

Error in (function (..., deparse.level = 1)  : 
  number of rows of matrices must match (see arg 2)

I tried using the block argument, but it simply returns a matrix that repeats 4 times a matrix with the 6 possible permutations of 3 values:

CTRL <- how(within=Within("free"),
            blocks=factor(thedata$judge),
            complete=TRUE,maxperm=1e9)
allPerms(12,control=CTRL)

IMPORTANT NOTE: I do have a custom function to obtain the result, using an adaptation of expand.grid() with permn() from the combinat package. I'm interested in where I misunderstand the permute package, not how I can calculate all these permutations myself.

Answer 1

The examples provided by @Joris identify two bugs in allPerms() that were not picked up by the current set of examples or unit tests (that'll be fixed soon too!).

The first issue is an obscure bug that I'll need a bit of time to think through a fix for. I have now implemented fixes for this bug too. Version 0.8-3 of permute now happily handles the Plots version of @joris' question:

R> p <- allPerms(12,control=CTRL)
R> dim(p)
[1] 1295   12
R> head(p)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    1    2    3    4    5    6    7    8    9    10    12    11
[2,]    1    2    3    4    5    6    7    8    9    11    10    12
[3,]    1    2    3    4    5    6    7    8    9    11    12    10
[4,]    1    2    3    4    5    6    7    8    9    12    10    11
[5,]    1    2    3    4    5    6    7    8    9    12    11    10
[6,]    1    2    3    4    5    6    7    9    8    10    11    12
R> packageVersion("permute")
[1] ‘0.8.3’

The second is an oversight. allPerms() generates permutation indices, but internally it works block by block. In the case @Joris reported, each block has 3 observations and hence 6 permutations of the indices 1:3. Once these permutation indices have been created, the code should have used them to index the row indices of the original data for each block. allPerms() was doing this for every conceivable combination of permutation types except the simple random permutation within blocks case. r2838 fixes this issue.

allPerms() was also not replicating each within-block permutation matrix to match each combination of rows in the other within-block permutation matrices. This requires an operation like expand.grid() but on the within-block permutation matrices. r2839 fixes this particular issue.

allPerms() works this way because it does not expect the within-block samples to be located contiguously within the original data series.

This second bug was fixed via r2838 and r2839 in the SVN sources on R-Forge.

R> require(permute)
Loading required package: permute
R> CTRL <- how(within=Within("free"),
+             blocks=factor(thedata$judge),
+             complete=TRUE,maxperm=1e9,
+             observed = TRUE)
R> numPerms(12,CTRL)
[1] 1296
R> tmp <- allPerms(12,control=CTRL)
R> dim(tmp)
[1] 1296   12
R> head(tmp)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    1    2    3    4    5    6    7    8    9    10    11    12
[2,]    1    2    3    4    5    6    7    8    9    10    12    11
[3,]    1    2    3    4    5    6    7    8    9    11    10    12
[4,]    1    2    3    4    5    6    7    8    9    11    12    10
[5,]    1    2    3    4    5    6    7    8    9    12    10    11
[6,]    1    2    3    4    5    6    7    8    9    12    11    10
R> tail(tmp)
        [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1291,]    3    2    1    6    5    4    9    8    7    10    11    12
[1292,]    3    2    1    6    5    4    9    8    7    10    12    11
[1293,]    3    2    1    6    5    4    9    8    7    11    10    12
[1294,]    3    2    1    6    5    4    9    8    7    11    12    10
[1295,]    3    2    1    6    5    4    9    8    7    12    10    11
[1296,]    3    2    1    6    5    4    9    8    7    12    11    10