1add.center <- function(design, ncenter, distribute=NULL, ...){ 2 if (!"design" %in% class(design)) stop("design must be of class design") 3 di <- design.info(design) 4 if (missing(ncenter)) stop("ncenter must be specified") 5 if (!(substr(di$type,1,4)=="FrF2" | substr(di$type,1,2)=="pb" | (substr(di$type,1,14)=="full factorial" & all(di$nlevels==2)))) 6 stop("center points only available for FrF2 and pb type designs") 7 if (di$type=="FrF2.splitplot") 8 stop("currently, center points for splitplot designs are not supported") 9 if (di$type=="FrF2.param") 10 stop("currently, center points for long version parameter designs are not supported") 11 if (length(grep("center",di$type))>0) 12 stop("design has center points already") 13 if (!is.numeric(ncenter)) stop("ncenter must be a number") 14 if (!length(ncenter)==1) stop("ncenter must be a number") 15 if (!ncenter==floor(ncenter)) stop("ncenter must be an integer number") 16 if (is.null(distribute)){ 17 if (ncenter==0 | !di$randomize) distribute <- 1 18 else distribute <- min(ncenter, 3)} 19 if (!is.numeric(distribute)) stop("distribute must be a number") 20 if (!distribute==floor(distribute)) stop("distribute must be an integer number") 21 if (distribute < 1 | distribute > max(1,min(di$nruns+1, ncenter))) 22 stop("distribute must be at least 1 and at most min(ncenter, nruns+1 of the design)") 23 if (di$randomize & distribute==1 & ncenter > 1) warning("running all center point runs together is usually not a good idea.") 24 25 if (any(is.na(sapply(factor.names(design),"is.numeric")))) 26 stop("Center points are implemented for experiments with all factors quantitative only.") 27 28 29 design <- qua.design(design, quantitative="all") 30 di <- design.info(design) 31 clevels <- sapply(factor.names(design), "mean") 32 di$ncube <- di$nruns 33 di$ncenter <- ncenter 34 if (!is.null(di$blocksize)) di$nruns <- di$nruns + di$nblocks*ncenter else 35 di$nruns <- di$nruns+ncenter 36 37 desnum <- desnum(design) 38 ro <- run.order(design) 39 if (!all(ro$run.no == sort(ro$run.no)) & (di$replications>0 | !is.null(di$blocksize))) 40 stop("center points cannot be added to designs that are not in original order") 41 hilf <- undesign(design) 42 43 ### 15 Feb 2013: removed restriction to bbreps==1 (as now always factor) 44 if (!is.null(di$blocksize)) block.contrasts <- contrasts(hilf[,di$block.name]) 45 46 ## positions for adding center points (added after these positions) 47 ## !!! blocks, replications and repeated measurements to be treated correctly 48 if (!is.null(di$blocksize)){ 49 ## distribute within blocks 50 if (distribute==1) 51 addpos <- di$blocksize 52 else 53 addpos <- c(0,round((1:(distribute-1))*di$blocksize/(distribute-1))) 54 } 55 else{ 56 if (distribute==1) 57 addpos <- di$ncube 58 else 59 addpos <- c(0,round((1:(distribute-1))*di$ncube/(distribute-1))) 60 } 61 ## determine numbers of center points to be added at each position 62 nrest <- ncenter%%distribute 63 if (nrest==0) nadd <- rep(round(ncenter/distribute),distribute) 64 if (nrest==1) nadd <- c(1,rep(0,distribute-1))+rep(floor(ncenter/distribute),distribute) 65 if (nrest==2) nadd <- c(1,rep(0,distribute-2),1)+rep(floor(ncenter/distribute),distribute) 66 if (nrest>2) nadd <- rep(floor(ncenter/distribute),distribute) + sample(c(rep(1,nrest),rep(0,distribute-nrest))) 67 ## repeated measurements also for center points (assuming they are done for measurement accuracy) 68 if (di$repeat.only) if (!is.null(di$blocksize)) nadd <- nadd*di$wbreps else nadd <- nadd*di$replications 69 70 ## split vector for design and run order data frame 71 if (!is.null(di$blocksize)){ 72 if (di$wbreps==1 | di$repeat.only) 73 together <- paste(rep(1:(di$nblocks*di$bbreps),each=di$blocksize*di$wbreps), 74 rep(cut(1:di$blocksize,unique(c(0,addpos))),each=di$wbreps,di$nblocks*di$bbreps),sep=".") 75 else 76 together <- paste(rep(1:(di$nblocks*di$bbreps*di$wbreps),each=di$blocksize), 77 rep(cut(1:di$blocksize,unique(c(0,addpos))),di$nblocks*di$bbreps*di$wbreps),sep=".") 78 } 79 else {if (di$replications==1 | (di$replications > 1 & di$repeat.only)) 80 together <- cut(rep(1:di$ncube,each=di$replications), unique(c(0,addpos))) 81 else together <- paste(rep(1:di$replications,each=di$ncube),rep(cut(1:di$ncube,unique(c(0,addpos))),di$replications),sep=".") 82 } 83 ## the following command prevents split from reordering the data 84 together <- factor(together, levels=unique(together)) 85 ## split the data frame 86 getrennt <- split(hilf,together) 87 ## added 15 Feb 2013 88 ## obtaine levels for the new factors including the center points 89 rostdlev <- levels(ro$run.no.in.std.order) 90 rostdrplev <- levels(ro$run.no.std.rp) 91 levgleich <- identical(ro$run.no.in.std.order,ro$run.no.std.rp) 92 93 ## extend levels to include all center point related ones 94 ## for run.no.in.std.order 95 if (all(sapply(rostdlev, function(obj) length(grep(".", obj, fixed=TRUE)) > 0 ) )){ 96 punkt1 <- sapply(rostdlev, function(obj) regexpr(".",obj,fixed=TRUE)) 97 if (all(sapply(substr(rostdlev,punkt1+1,999), function(obj) length(grep(".", obj, fixed=TRUE))>0))) 98 punkt2 <- sapply(substr(rostdlev, punkt1+1, 999), function(obj) regexpr(".", obj, fixed=TRUE))+punkt1 99 else punkt2 <- 1000 ## repeated measurements or replications only 100 if (punkt2[1] < 1000) rostdlev <- c(unique(paste("0", substr(rostdlev, punkt1+1, punkt2-1), "0", sep=".")), 101 rostdlev) else rostdlev <- c(unique(paste("0", substr(rostdlev, punkt1+1, punkt2-1), sep=".")), 102 rostdlev) 103 } 104 else rostdlev <- c("0", rostdlev) 105 ## for run.no.std.rp 106 if (levgleich) rostdrplev <- rostdlev else{ 107 if (all(sapply(rostdrplev, function(obj) length(grep(".", obj, fixed=TRUE))>0))){ 108 punkt1 <- sapply(rostdrplev, function(obj) regexpr(".",obj,fixed=TRUE)) 109 if (all(sapply(substr(rostdrplev,punkt1+1,999), function(obj) length(grep(".", obj, fixed=TRUE))>0))){ 110 punkt2 <- sapply(substr(rostdrplev, punkt1+1, 999), function(obj) regexpr(".", obj, fixed=TRUE))+punkt1 111 punkt3 <- sapply(substr(rostdrplev, punkt2+1, 999), function(obj) regexpr(".", obj, fixed=TRUE))+punkt2 112 if (all(sapply(substr(rostdrplev,punkt3+1,999), function(obj) length(grep(".", obj, fixed=TRUE))>0))) 113 punkt4 <- sapply(substr(rostdrplev, punkt3+1, 999), function(obj) regexpr(".", obj, fixed=TRUE))+punkt3 else 114 punkt4 <- 1000 115 if (punkt4[1] == 1000) rostdrplev <- c(unique(paste("0", 116 substr(rostdrplev, punkt1+1, punkt2-1), "0", 117 substr(rostdrplev, punkt3+1, punkt4-1), sep=".")), rostdrplev) 118 else {if (di$repeat.only) rostdrplev <- c(unique(paste("0", 119 substr(rostdrplev, punkt1+1, punkt2-1), "0", 120 substr(rostdrplev, punkt3+1, punkt4-1), 1:di$wbreps, sep=".")), rostdrplev) 121 else rostdrplev <- c(unique(paste("0", 122 substr(rostdrplev, punkt1+1, punkt2-1), "0", 123 substr(rostdrplev, punkt3+1, punkt4-1), substr(rostdrplev, punkt4+1, 999), sep=".")), rostdrplev) 124 } 125 } 126 else {punkt2 <- 1000 ## repeated measurements or replications only 127 rostdrplev <- c(unique(paste("0", substr(rostdrplev, punkt1+1, punkt2-1), 128 sep=".")), rostdrplev) 129 } 130 } 131 else rostdrplev <- c("0", rostdrplev) 132 } 133 ## end of added 15 Feb 2013 134 ro.getrennt <- split(ro,together) 135 ## wird das �berhaupt gebraucht??? 136 ro$run.no.in.std.order <- as.character(ro$run.no.in.std.order) 137 ro$run.no.std.rp <- as.character(ro$run.no.std.rp) 138 139 blockid <- rep("1", length(getrennt)) 140 141 von <- cumsum(c(1,lapply(getrennt, "nrow")))[-(length(getrennt)+1)] 142 bis <- cumsum(lapply(getrennt, "nrow")) 143 144 ## added 15 Feb 2013 145 std <- lapply(ro.getrennt, function(obj) { 146 hilf <- t(as.matrix(sapply(obj$run.no.in.std.order, function(obj2) unlist(strsplit(as.character(obj2),".",fixed=TRUE))))) 147 if (nrow(hilf)==1) hilf <- t(hilf) 148 hilf 149 }) 150 stdrp <- lapply(ro.getrennt, function(obj) { 151 hilf <- t(as.matrix(sapply(obj$run.no.std.rp, function(obj2) unlist(strsplit(as.character(obj2),".",fixed=TRUE))))) 152 if (nrow(hilf)==1) hilf <- t(hilf) 153 hilf 154 }) 155 ## get rid of rp-entry 156 if (is.null(di$blocksize) & di$replications > 1 & di$repeat.only) 157 stdrp <- lapply(stdrp, function(obj) obj[,1:(ncol(obj)-1)]) 158 if (!is.null(di$blocksize)){ 159 blockid <- sapply(getrennt, function(obj) as.character(obj[1,di$block.name])) 160 ## get rid of rp-entry in stdrp for blocked designs 161 if (di$wbreps > 1 & di$repeat.only) 162 stdrp <- lapply(stdrp, function(obj) obj[,1:(ncol(obj)-1)]) 163 } 164 ## end of added 15 Feb 2013 165 if (!is.null(di$blocksize)) 166 blockid <- sapply(getrennt, function(obj) as.character(obj[1, 167 di$block.name])) 168 if (di$replications > 1 & is.null(di$blocksize) & !di$repeat.only){ 169 blockid <- rep(1:di$replications, each=max(1,distribute-1)) 170 } 171 if (!is.null(di$blocksize)) if (di$wbreps > 1 & !di$repeat.only){ 172 blockid <- paste(blockid, rep(1:di$wbreps, 173 each=max(1,distribute-1), 174 times=di$nblocks*di$bbreps),sep=".") 175 } 176 ## create center point data frames to be interleaved 177 # if (!is.null(di$nblocks)) no.center.groups <- distribute*di$nblocks*di$bbreps else 178 # no.center.groups <- distribute*di$replications 179 centers <- vector("list") 180 ros <- lapply(1:length(std), function(obj) vector("list")) 181 182 ## centers and ros for WITHIN each block 183 more <- setdiff(colnames(hilf), c(di$block.name, names(di$factor.names))) 184 ## columns other than design factors, e.g. responses 185 if (ncenter>0){ 186 for (i in 1:distribute){ 187 cnext <- data.frame(matrix(clevels, ncol=di$nfactors, nrow=nadd[i], 188 byrow=TRUE, dimnames=list(NULL, names(di$factor.names)))) 189 if (length(more)>0) cnext <- cbind(cnext, matrix(NA,nrow=nadd[i],ncol=length(more),dimnames=list(NULL,more))) 190 centers <- c(centers, list(cnext)) 191 192 ## simultaneous creation of inserts for all blocks via lapply 193 if (!is.null(di$blocksize)){ 194 stdnext <- lapply(1:length(std), function(obj) rep(paste(c(0,std[[obj]][1,2],0),collapse="."), nadd[i])) 195 stdrpnext <- stdnext 196 if (di$bbreps>1 | (di$wbreps>1 & !di$repeat.only)) stdrpnext <- 197 mapply(paste, stdrpnext, lapply(stdrp, function(obj) paste(obj[1,4], collapse=".")), 198 sep=".", SIMPLIFY=FALSE) 199 if (di$bbreps>1 & di$wbreps>1 & !di$repeat.only) stdrpnext <- 200 mapply(paste, stdrpnext, lapply(stdrp, function(obj) paste(obj[1,5], collapse=".")), 201 sep=".", SIMPLIFY=FALSE) 202 if (di$wbreps>1 & di$repeat.only) 203 stdrpnext <- lapply(stdrpnext, function(obj) paste(obj, rep(1:di$wbreps, round(nadd[i]/di$wbreps)),sep=".")) 204 } 205 else { 206 stdnext <- stdrpnext <- lapply(1:length(std), function(obj) rep(0, nadd[i]) ) 207 if (di$replications > 1){ 208 if (!di$repeat.only) 209 stdrpnext <- mapply(paste, stdrpnext, lapply(stdrp, function(obj) obj[1,2]),sep=".", SIMPLIFY=FALSE) 210 else stdrpnext <- lapply(stdrpnext, function(obj) paste(obj, rep(1:di$replications,round(nadd[i]/di$replications)),sep=".")) 211 } 212 } 213 stdnext <- lapply(stdnext, function(obj) factor(obj, levels=rostdlev)) 214 stdrpnext <- lapply(stdrpnext, function(obj) factor(obj, levels=rostdrplev)) 215 bothnext <- mapply(function(...) list(data.frame(...)), stdnext, lapply(1:length(std), function(obj) rep(0,nadd[i])), 216 stdrpnext, SIMPLIFY=FALSE) 217 bothnext <- lapply(bothnext, function(obj){colnames(obj[[1]]) <- c("run.no.in.std.order","run.no","run.no.std.rp"); obj}) 218 ros <- mapply(c,ros, bothnext, SIMPLIFY=FALSE) 219 } 220 ## interleave 221 ## first center point entry of first (replication) block 222 if (!is.null(di$blocksize)) 223 new <- cbind(matrix(hilf[1,di$block.name],nrow=nadd[1],ncol=1, dimnames=list(NULL, di$block.name)),centers[[1]]) 224 else new <- centers[[1]] 225 ronew <- ros[[1]][[1]] 226 blockids <- unique(blockid) 227 if (distribute==1){ 228 ## center points appended at end of each (replication) block 229 if (length(blockids)==1){ 230 ## one (replication) block only 231 new <- rbind(hilf, new) 232 ronew <- rbind(ro, ronew) 233 } 234 else{ 235 ## more than one (replication) block 236 ## as distribute equals 1, blockids == blockid 237 new <- rbind(getrennt[[1]], new) 238 ronew <- rbind(ro.getrennt[[1]], ronew) 239 for (i in 2:length(blockids)){ 240 if (!is.null(di$blocksize)) 241 new <- rbind(new, getrennt[[i]], 242 cbind(matrix(getrennt[[i]][1,di$block.name],nrow=nadd[1],ncol=1, dimnames=list(NULL, di$block.name)),centers[[1]])) 243 else new <- rbind(new, getrennt[[i]], centers[[1]]) 244 ronew <- rbind(ronew, ro.getrennt[[i]], ros[[i]][[1]]) 245 } 246 } 247 } 248 else{ 249 ### distribute > 1 250 for (i in 1:length(blockids)){ 251 iblock <- which(blockid==blockids[i]) 252 if (i>1){ 253 ## append first center points series of new block 254 ## for first block done already 255 if (!is.null(di$blocksize)) 256 new <- rbind(new, 257 cbind(matrix(getrennt[[iblock[1]]][1,di$block.name],nrow=nadd[1],ncol=1, dimnames=list(NULL, di$block.name)),centers[[1]])) 258 else new <- rbind(new, centers[[1]]) 259 ronew <- rbind(ronew, ros[[iblock[1]]][[1]]) 260 } 261 for (j in 1:length(iblock)){ 262 ## append next data with subsequent center points for the current (replication) block 263 if (!is.null(di$blocksize)) 264 new <- rbind(new, getrennt[[iblock[j]]], 265 cbind(matrix(getrennt[[iblock[j]]][1,di$block.name],nrow=nadd[j+1],ncol=1, dimnames=list(NULL, di$block.name)),centers[[j+1]])) 266 else new <- rbind(new, getrennt[[iblock[j]]], centers[[j+1]]) 267 ronew <- rbind(ronew, ro.getrennt[[iblock[j]]], ros[[iblock[1]]][[j+1]]) 268 } 269 } 270 } 271 if (!is.null(di$blocksize)){ 272 new[[di$block.name]] <- factor(new[[di$block.name]], levels(design[[di$block.name]])) 273 if (exists("block.contrasts")) contrasts(new[,di$block.name]) <- block.contrasts 274 desnum <- model.matrix(formula(paste("~",paste(c(di$block.name,names(di$factor.names)),collapse="+"))),data=new)[,-1] 275 } 276 else 277 desnum <- model.matrix(formula(paste("~",paste(names(di$factor.names),collapse="+"))),data=new)[,-1] 278 if (length(setdiff(colnames(new),c(di$block.name, names(di$factor.names))))>0){ 279 anhaeng <- as.matrix(new[,setdiff(colnames(new),c(di$block.name, names(di$factor.names))),drop=FALSE]) 280 storage.mode(anhaeng) <- "numeric" 281 desnum <- cbind(desnum, anhaeng) 282 } 283 284 rownames(new) <- rownames(desnum) <- ronew$run.no <- 1:nrow(new) 285 class(new) <- c("design","data.frame") 286 desnum(new) <- desnum 287 288 run.order(new) <- ronew 289 di$type <- paste(di$type,"center",sep=".") 290 ## added in order to support usage of function code.design and steepest ascent analysis 291 di$coding <- make.formulas(paste("x", 1:length(di$factor.names), sep=""), di$factor.names) 292 design.info(new) <- di 293 } 294 else { 295 ### ??? was wird noch ben�tigt f�r ccd.augment??? 296 new <- design 297 ronew <- run.order(design) 298 if (!any(is.na(as.numeric(as.character(ronew$run.no.in.std.order))))) 299 ronew$run.no.in.std.order <- as.numeric(as.character(ronew$run.no.in.std.order)) 300 design.info(new) <- c(design.info(new), ncenter=0, ncube=nrow(new)) 301 run.order(new) <- ronew 302 } 303 new 304}