[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.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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.
Note 2: The 9 treatment design? It was found via group theory. Apn gave a link Some New Row-Complete Latin Squares, Archdeacon, Dinitz, Stinson and Tillson
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)
}
To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.