1blockpick <- function(k, gen, k.block, design=NULL, show=10, alias.block.2fis=FALSE, select.catlg=catlg){ 2 ## function that picks design with right number of blocks (power of 2, 2^k.block) 3 4 ## k.block is the number of independent block factors needed for block construction 5 ## (log2 of number of blocks) 6 ## usual assumption: their interactions with treatment factors can be neglected 7 8 ## 2^k.block-1 columns are used exclusively for blocks 9 10 ## if requested (default), alias.block.2fis=FALSE stops confounding of treatment 2fis with blocks 11 ## if this is not requested, no attempt is made to reduce the number of aliased 2fis, 12 ## but they are reported 13 14 if (!is.logical(alias.block.2fis)) stop("alias.block.2fis must be a logical.") 15 if (!is.null(design)){ 16 if (!is.character(design)) stop("If specified, design must be a character string.") 17 if (!design %in% names(select.catlg)) stop("If specified, design must be the name of a design in select.catlg.") 18 if (!missing(gen)) warning("If design is specified, gen is ignored.") 19 } 20 if (is.null(design) & missing(gen)) stop("One of gen or design must be given.") 21 if (!"catlg" %in% class(select.catlg)) stop("select.catlg must be a catalogue of designs of class catlg.") 22 23 if (is.character(design)) { 24 if (!design %in% names(select.catlg)) stop ("Design must be a valid design name.") 25 if (select.catlg[[design]]$res == 3 & !alias.block.2fis) 26 stop("Resolution III design. Package FrF2 does not allow to request all 2fis to be clear of blocks.") 27 if (!select.catlg[[design]]$nruns == 2^k) stop("mismatch between k and design") 28 gen <- select.catlg[[design]]$gen 29 res <- select.catlg[[design]]$res 30 } 31 if (identical(gen,0) | length(gen)==0) { 32 ## treat full factorials by catalogue (block.catlg) 33 gen <- numeric(0) 34 if (k <= max(block.catlg$k) & k.block <= ncol(block.catlg)-2){ 35 if (length(which(block.catlg$k==k & block.catlg$k.block==k.block))>0) 36 bcols <- block.catlg[which(block.catlg$k==k & block.catlg$k.block==k.block)[1], 37 paste("b",1:k.block,sep="")] 38 else stop("This full factorial is not in the catalogue.") 39 } 40 else bcols <- catlg[[paste(k+k.block,"-",k.block,".",1,sep="")]]$gen 41 } 42 g <- length(gen) 43 minus <- 1 44 if (g>0){ 45 minus <- which(gen<0) 46 gen <- abs(gen) 47 } 48 49 hilf <- c(k,gen,k.block,show) 50 if (!is.numeric(hilf)) 51 stop ("k, gen, k.block and show must be numeric.") 52 if (!all(hilf == floor(hilf) & hilf > 0)) 53 stop ("k, gen, k.block and show must contain integer numbers.") 54 if (!k >= 3) stop ("blockpick requires k>=3.") 55 if (!k.block < k) stop ("blockpick requires k.block < k.") 56 if (any(gen %in% 2^(0:(k-1)))) 57 stop ("gen must not contain column numbers of base factors.") 58 if (any(!gen %in% 3:(2^k-1))) 59 stop ("Column numbers in gen must be smaller than ", 2^k,".") 60 61 ## assignment of k.block factors to the remaining columns 62 ## 63 hilf <- cols(k, gen) 64 fi2s <- hilf$fi2s 65 names(fi2s) <- apply(combn(k+g,2),2,function(obj) paste(Letters[obj],collapse="")) 66 minus2fis <- which(apply(combn(k+g,2),2,function(obj) length(obj %in% minus)==1)) 67 names(fi2s)[minus2fis] <- paste("-",names(fi2s)[minus2fis],sep="") 68 nfi2s <- table(fi2s) 69 fi2cols <- as.numeric(names(nfi2s)) 70 if (is.null(design)) { 71 if (any(2^(0:(k-1)) %in% fi2cols)) res <- 3 72 else if (any(nfi2s>1)) res <- 4 73 else res="5+" 74 } 75 if (alias.block.2fis) {banned <- hilf$main; eligible <- c(fi2cols, hilf$freecols)} 76 else {banned <- c(hilf$main,fi2cols); eligible <- hilf$freecols} 77 if (length(eligible)<2^k.block-1 & !alias.block.2fis) 78 stop("no adequate block design found with 2fis unconfounded with blocks") 79 if (length(eligible)<2^k.block-1) stop("no adequate block design found") 80 if (g>0) 81 perm <- t(combn(length(eligible),k.block)) ## rows are possible selections of k.block 82 else { 83 if (g==0 & !alias.block.2fis) if (!all(bcols %in% eligible)) 84 stop("no adequate block design found with 2fis unconfounded with blocks") 85 perm <- matrix(which(eligible %in% bcols),nrow=1) 86 } 87 88 hilf <- digitsBase(eligible,ndigits=k) ## always k rows 89 ## potential block columns as binary numbers 90 ergeb <- matrix(0,nrow=nrow(perm),ncol=2^k.block-1) 91 ## will contain generating and main effects columns for blocks 92 banned.block <- 99 93 ## will contain number of block main effects aliased with banned columns 94 dependent.block <- TRUE 95 ## will contain number of block main effects aliased with banned columns 96 n2fis.block <- NA 97 ## number of 2fis aliased with block main effects 98 n2fis.clear <- NA 99 ## number of 2fis aliased with block main effects 100 hilfc <- Yates[1:(2^k.block-1)] 101 ## combination patterns of block factors into all block main effects 102 count <- 0 ## number of possibilities found 103 i <- 0 ## current position in vectors 104 stopp <- FALSE 105 last <- numeric(0) 106 107 ## loop through selections of block columns 108 while (!stopp){ 109 perm <- perm[setdiff(1:nrow(perm),last),,drop=FALSE] 110 i <- i+1 111 bg <- eligible[perm[1,]] ## block generators 112 bcols <- sapply(lapply(hilfc, function(obj) mult.gen(Yates[bg[obj]])), 113 function(obj) {h2 <- which(sapply(Yates[1:(2^k-1)],function(obj2){ 114 hh <- FALSE 115 if (length(obj2)==length(obj)) 116 if (all(obj2==obj)) hh <- TRUE 117 hh})) 118 if (length(h2)==0) 0 else h2 119 } 120 ) 121 122 last <- which(apply(perm,1,function(obj) all(eligible[obj] %in% bcols))) 123 ### remove redundant rows from perm 124 ergeb[i,] <- bcols ## 2^k.block-1 potential block main effect columns 125 126 if (any(bcols==0)) next 127 ## dependent.block is TRUE per default and remains so for bcols with 0 128 banned.block[i] <- sum(bcols %in% banned) ## generated block main effects confounded with banned columns 129 dependent.block[i] <- length(table(bcols)) < length(bcols) ## k.block columns not independent 130 n2fis.block[i] <- sum(nfi2s[which(as.numeric(names(nfi2s)) %in% bcols)]) 131 n2fis.clear[i] <- length(nfi2s[which(nfi2s==1 & !as.numeric(names(nfi2s)) %in% c(2^(0:(k-1)),gen,bcols))]) 132 if (banned.block[i] == 0 & !dependent.block[i]){ 133 count <- count + 1 134 if (count>=show) break 135 } 136 if (nrow(perm) <= length(last)) stopp <- TRUE 137 } ## end of loop over permutations 138 139 pick <- which(banned.block==0 & !dependent.block) 140 141 if (length(pick)<show) show <- length(pick) 142 if (show==0 & !alias.block.2fis) stop("no adequate block design found with 2fis unconfounded with blocks") 143 if (show==0) stop("no adequate block design found") else { 144 ntreat <- k+length(gen) 145 blocks <- 2^k.block 146 blockcols <- ergeb[pick[1:show],,drop=FALSE] 147 if (alias.block.2fis){ 148 alias.2fis.block <- vector("list",show) 149 for (i in 1:show) 150 alias.2fis.block[[i]] <- names(fi2s)[which(abs(fi2s) %in% blockcols[i,])] 151 } 152 else alias.2fis.block <- "none" 153 nclear.2fis <- n2fis.clear[pick[1:show]] 154 nblock.2fis <- n2fis.block[pick[1:show]] 155 if (sum(nclear.2fis)>0) { 156 clear.2fis <- vector("list",show) 157 for (i in 1:show) 158 clear.2fis[[i]] <- names(fi2s)[which(fi2s %in% 159 setdiff(as.numeric(names(nfi2s))[which(nfi2s==1)], c(2^(0:(k-1)),gen,blockcols[i,])))] 160 } 161 else clear.2fis <- character(0) 162 gen[minus] <- -gen[minus] 163 list(gen=gen, 164 basics = c(nruns=2^k,nblocks=blocks, ntreat=ntreat, res.base=res), 165 blockcols=blockcols[,2^(0:(k.block-1))], 166 alias.2fis.block=alias.2fis.block, 167 nblock.2fis=nblock.2fis, 168 nclear.2fis=nclear.2fis, 169 clear.2fis=clear.2fis) 170 } 171} 172 173