Carry-over balanced designs for 8 treatments
[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.
Those are Williams designs you might say, but it has become clear to me that Williams designs are just a subset of all carry-over balanced designs. Not through hard work of mine, comments by Apn on my previous post Creating Williams designs with even number of products lead to this. Among other things, Apn claimed there were designs which did not confirm to the well known symmetry that each row is the reflection of another row. Related to that was a claim that there is a nine treatment square carry-over balanced design. I made my algorithm faster and in this post I can confirm the former, but the latter is too long a calculation for the algorithm I used.Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Designs
There are 12 designs, as given below. Numbers 4 and 5 can be obtained from each other by certain recoding of the treatments. The last four are the non-symmetrical ones.[[1]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 4 1 6 3 8 5 7
[3,] 3 1 5 2 7 4 8 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 8 4 7 2 5 1 3
[7,] 7 5 8 3 6 1 4 2
[8,] 8 7 6 5 4 3 2 1
[[2]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 4 8 3 6 1 5 7
[3,] 3 1 4 7 2 5 8 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 8 5 2 7 4 1 3
[7,] 7 5 1 6 3 8 4 2
[8,] 8 7 6 5 4 3 2 1
[[3]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 4 1 3 6 8 5 7
[3,] 3 8 4 7 2 5 1 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 1 5 2 7 4 8 3
[7,] 7 5 8 6 3 1 4 2
[8,] 8 7 6 5 4 3 2 1
[[4]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 4 8 6 3 1 5 7
[3,] 3 8 5 2 7 4 1 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 1 4 7 2 5 8 3
[7,] 7 5 1 3 6 8 4 2
[8,] 8 7 6 5 4 3 2 1
[[5]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 5 1 6 3 8 4 7
[3,] 3 1 4 2 7 5 8 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 8 5 7 2 4 1 3
[7,] 7 4 8 3 6 1 5 2
[8,] 8 7 6 5 4 3 2 1
[[6]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 5 8 3 6 1 4 7
[3,] 3 1 5 7 2 4 8 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 8 4 2 7 5 1 3
[7,] 7 4 1 6 3 8 5 2
[8,] 8 7 6 5 4 3 2 1
[[7]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 5 1 3 6 8 4 7
[3,] 3 8 5 7 2 4 1 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 1 4 2 7 5 8 3
[7,] 7 4 8 6 3 1 5 2
[8,] 8 7 6 5 4 3 2 1
[[8]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 5 8 6 3 1 4 7
[3,] 3 8 4 2 7 5 1 6
[4,] 4 6 2 8 1 7 3 5
[5,] 5 3 7 1 8 2 6 4
[6,] 6 1 5 7 2 4 8 3
[7,] 7 4 1 3 6 8 5 2
[8,] 8 7 6 5 4 3 2 1
[[9]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 7 1 8 3 5 4 6
[3,] 3 1 5 7 6 8 2 4
[4,] 4 8 7 5 2 1 6 3
[5,] 5 3 6 2 8 4 1 7
[6,] 6 5 8 1 4 7 3 2
[7,] 7 4 2 6 1 3 8 5
[8,] 8 6 4 3 7 2 5 1
[[10]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 7 6 3 8 1 5 4
[3,] 3 6 8 5 2 4 1 7
[4,] 4 3 5 7 1 8 6 2
[5,] 5 8 2 1 3 7 4 6
[6,] 6 1 4 8 7 3 2 5
[7,] 7 5 1 6 4 2 8 3
[8,] 8 4 7 2 6 5 3 1
[[11]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 8 5 1 7 3 6 4
[3,] 3 5 4 6 1 8 2 7
[4,] 4 1 6 8 3 7 5 2
[5,] 5 7 1 3 2 4 8 6
[6,] 6 3 8 7 4 2 1 5
[7,] 7 6 2 5 8 1 4 3
[8,] 8 4 7 2 6 5 3 1
[[12]]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 2 3 4 5 6 7 8
[2,] 2 8 5 7 4 1 3 6
[3,] 3 5 2 6 8 7 1 4
[4,] 4 7 6 2 1 5 8 3
[5,] 5 4 8 1 6 3 2 7
[6,] 6 1 7 5 3 8 4 2
[7,] 7 3 1 8 2 4 6 5
[8,] 8 6 4 3 7 2 5 1
Calculation speed
Up to 5 treatments 100 repeats are done in the blink of an eye. For 6 it is doable (0.025 sec per run). For 7 treatments I did 10 repeats (1.6 sec per run), while for 8 one run took almost 800 sec. Clearly 9 treatments is too much.> microbenchmark(gendesign(3),gendesign(4),gendesign(5),gendesign(6))
Note: no visible global function definition for ‘error’
Unit: microseconds
expr min lq median uq max neval
gendesign(3) 107.029 113.2605 119.492 125.3565 240.449 100
gendesign(4) 343.081 355.5430 369.472 405.7585 544.677 100
gendesign(5) 1141.403 1163.0290 1202.615 1258.6955 1472.755 100
gendesign(6) 24753.851 25352.7755 25595.058 26917.8980 35721.433 100
> system.time(gendesign(7))
user system elapsed
1.64 0.00 1.67
> microbenchmark(gendesign(7),times=10L)
Unit: seconds
expr min lq median uq max neval
gendesign(7) 1.624484 1.647051 1.657468 1.669608 1.768613 10
> system.time(gendesign(8))
user system elapsed
793.72 0.17 805.82
Code
Code is below. Running the JIT saves half of the time. It is basically a simplified version of my old algorithm, but with a lot of small modifications, basically avoiding loops, ifs and a load of trial and error for speed.
library(microbenchmark)
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
carover <- matrix(TRUE,nrow=nr,ncol=nc)
for (i in 1L:(nc-1L)) carover[i,i+1] <- FALSE
todo <- nextpos(desmat)
desobject <- list(desmat=desmat,carover = carover,nc=1L:n,n=n,
index =1L,npos=nrow(todo),
row=todo[,1L],
col=todo[,2L])
desresult <- list()
addpoint(desobject,desresult)
}
modify <- function(desobject,row,col,i,previous) {
desobject$desmat[row,col] <- i
desobject$carover[previous,i] <- FALSE
desobject$index <- desobject$index + 1L
desobject}
addpoint <- function(desobject,desresult) {
if (desobject$index>desobject$npos) {
l <- length(desresult)
desresult[[l+1]] <- desobject$desmat
# cat(‘#’)
return(desresult)
}
row <- desobject$row[desobject$index]
col <- desobject$col[desobject$index]
previous <- desobject$desmat[row,col-1L]
avoid <- c(desobject$desmat[row,],
desobject$desmat[,col])
nc <- desobject$nc[!is.element(desobject$nc,avoid) ]
nc <- nc[desobject$carover[previous,nc]]
for (i in nc) {
desresult <- addpoint(modify(desobject,row,col,i,previous)
,desresult)
}
desresult
}
library(compiler)
enableJIT(3)
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.