[This article was first published on

** Wiekvoet**, and kindly contributed to

R-bloggers]. (You can report issue about the content on this page

here)

Want to share your content on R-bloggers?

click here if you have a blog, or

here if you don't.

Based on some comments I am looking at using symmetry to obtain Williams style designs. Symmetry allows reduction of the number of combinations examined hence faster calculation times. Two avenues are examined. Both work for a low number of treatments. For eight treatments they give different sets, which together define all 12 possible designs.

### Horizontal reflection

It is well known that each row in an even number of treatments Williams design is also present reversed. But apn noticed something else, within a row there is a horizontal symmetry; the values of the cell i and n+1-i add up to n+1. That is a bit abstract. Lets look at a six treatment design. Cell [2,2] is 4 and cell [2,5] is 3, they add to 7. And so it continues. If this holds, then the permutations to check is decreased enormously. Especially, since cell [2,2] also by the reversal of rows is equal to cell [5,5] and hence cell [5,2] is 3.

[,1] [,2] [,3] [,4] [,5] [,6]

[1,] 1 2 3 4 5 6

[2,] 2 4 1 6 3 5

[3,] 3 1 5 2 6 4

[4,] 4 6 2 5 1 3

[5,] 5 3 6 1 4 2

[6,] 6 5 4 3 2 1

After putting this in an algorithm it appeared it worked. Run time for 8 treatments only 100 miliseconds.

microbenchmark(gendesign(4),gendesign(6),gendesign(8))

Unit: milliseconds

expr min lq median uq max neval

gendesign(4) 1.576118 1.609474 1.638796 1.717969 2.896391 100

gendesign(6) 8.695041 8.876843 9.054615 9.229821 14.714345 100

gendesign(8) 100.170001 102.489826 106.144955 108.703397 121.015079 100

The algorithm gave 8 designs for 8 treatments, 192 for 10 treatments, same as

previously.

### Symmetrical

I am also looking for a smart way to tackle designs for odd products. There is a nine treatment solution, I want to be able to find it. On a whim, I decided to check if symmetrical designs give is the solution. It is not. However, this approach can also be used for the even number of treatments. Again there are 8 designs for 8 treatments. Four of these have the horizontal semi-reflection, four do not have it. Hence all 12 designs which I found previously by brute force are covered. With 13 seconds run time for 8 treatments, I did not check for the number of solutions for 10 treatments..

Note 1: I do know symmetrical has no meaning in a Williams design. I have only used the designs with rows as subjects, columns as periods/rounds and numbers as treatments. And I do randomize subjects over rows, treatments over numbers. Symmetry is thus just a means.

### Code

The two sets of code are quite similar. This is only logical since one is derived from the other. Each section should be able to be run separately as displayed here.

#### Horizontal semi reflection

nextpos <- function(desmat) which(desmat==0,arr.ind=TRUE)

gendesign <- function(n=6) {

nr <- as.integer(n)

nc <- nr

desmat <- matrix(0L,nrow=nr,ncol=nc)

desmat[1,] <- 1L:nc

desmat[,1] <- 1L:nr

desmat[n,] <- nc:1L

desmat[,n] <- nr:1L

desobject <- list(desmat=desmat,nc=1L:n,n=n,np=n+1L)

desresult <- list()

addpoint(desobject,desresult)

}

modify <- function(desobject,row,col,i) {

desobject$desmat[desobject$np-row,desobject$np-col] <-

desobject$desmat[row,col] <- i

desobject$desmat[desobject$np-row,col] <-

desobject$desmat[row,desobject$np-col] <- desobject$np-i

desobject}

addpoint <- function(desobject,desresult) {

if (!cocheck(desobject)) return(desresult)

rc <- nextpos(desobject$desmat)

if (length(rc)==0) {

if (var(colSums(desobject$desmat))==0) {

l <- length(desresult)

desresult[[l+1]] <- desobject$desmat

}

return(desresult)

}

row <- rc[1,1]

col <- rc[1,2]

avoid <- c(desobject$desmat[row,],

desobject$desmat[,col])

nc <- desobject$nc[!is.element(desobject$nc,avoid) ]

for (i in nc) {

desresult <- addpoint(modify(desobject,row,col,i)

,desresult)

}

desresult

}

cocheck <- function(desobject) {

dm2 <- desobject$desmat*desobject$n

dm2 <- dm2[,-1]

dm3 <- desobject$desmat[,-desobject$n]

dm4 <- dm2+dm3

dm4 <- c(dm4[dm2!=0 & dm3!=0])

all(table(dm4)==1)

}

#### symmetrical

nextpos <- function(desmat) which(desmat==0,arr.ind=TRUE)

gendesign <- function(n=5) {

nr <- as.integer(n)

nc <- nr

desmat <- matrix(0L,nrow=nr,ncol=nc)

desmat[1,] <- 1L:nc

desmat[,1] <- 1L:nr

desobject <- list(desmat=desmat,nc=1L:n,n=n,np=n+1L)

desresult <- list()

addpoint(desobject,desresult)

}

modify <- function(desobject,row,col,i) {

desobject$desmat[col,row] <-

desobject$desmat[row,col] <- i

desobject}

addpoint <- function(desobject,desresult) {

if (!cocheck(desobject)) return(desresult)

rc <- nextpos(desobject$desmat)

if (length(rc)==0) {

if (var(colSums(desobject$desmat))==0) {

l <- length(desresult)

desresult[[l+1]] <- desobject$desmat

}

return(desresult)

}

row <- rc[1,1]

col <- rc[1,2]

avoid <- c(desobject$desmat[row,],

desobject$desmat[,col])

nc <- desobject$nc[!is.element(desobject$nc,avoid) ]

for (i in nc) {

desresult <- addpoint(modify(desobject,row,col,i)

,desresult)

}

desresult

}

cocheck <- function(desobject) {

dm2 <- desobject$desmat[,-1]*desobject$n

dm3 <- desobject$desmat[,-desobject$n]

dm4 <- dm2+dm3

dm4 <- c(dm4[dm2!=0 & dm3!=0])

all(table(dm4)==1)

}

*Related*

If you got this far, why not

__subscribe for updates__ from the site? Choose your flavor:

e-mail,

twitter,

RSS, or

facebook...