R/MADII_layout_function.R

   1 ### Function for generating a randomized, MADII experimental design
   2 ## Tyler Tiede. University of Minnesota. March 7, 2014
   3 #cat (" \n** Before using MADIIdgn be sure
   4 #that there are no files and/or folders in your Working Directory named\nthe same as
   5 #the character string you plan to designate for 'designID'; they will be overwritten.\n ")
   6
   7 #List of checks starting with the primary; assumes you have a primary check
   8
   9 MADIIdgn <- function(enviro="Eretz",
  10                      entries= NULL,
  11                      num.entries= NULL,
  12                      chk.names= NULL,
  13                      num.sec.chk= NULL,
  14                      num.rows= NULL,
  15                      num.cols= NULL,
  16
  17                      plot.start = 1001,
  18                      designID=NULL,
  19                      annoy=F)
  20   {
  21
  22   ## Ensure that no folders/files in WD will be overwritten and then create folder for results
  23   #JL "is.null(designID)==F"  is equivalent to "!is.null(designID)"
  24   if(annoy==T & is.null(designID)==F){
  25   YorN <- NA
  26     while(YorN %in% c("y","n","Y","N","yes","no","YES","NO","Yes","No") == F){
  27       YorN <- readline("Confirm there is Not a folder and/or file in your Working Directory with\nthe same name as 'designID' (type 'y' to proceed, or 'n' to stop): ")
  28     }
  29     #JL It doesn't look like typing 'n' will actually stop the function...
  30     # I would put in something like if(YorN %in% c("n", "N", "no", "No", "NO")) stop("User requested stop")
  31     if(YorN %in% c("y","Y","yes","YES","Yes")==T){
  32       unlink(designID, recursive=T)
  33       dir.create(path=designID, recursive=F)
  34     }
  35   } else{
  36     #if(is.null(designID)==F){
  37       unlink(designID, recursive=T)
  38       dir.create(path=designID, recursive=F)
  39     #}
  40   }
  41
  42
  43   ## Load necessary packag#es
  44   require("grid", quietly=T) ; library(grid)
  45
  46   ## QC of function inputs
  47   #JL why are these only warnings and not errors that cause the function to stop?
  48   if(is.null(entries) & is.null(num.entries)){
  49     warning("Must provide an entry list (entries=) OR the number of entries desired (num.entries=).")
  50   }
  51
  52   if(is.null(chk.names) & is.null(num.sec.chk)){
  53     warning("Must provide a list of check names (chk.names=) with the primary check listed first\n OR the number of SECONDARY checks desired (num.sec.chk=).")
  54   }
  55
  56   if(is.null(num.rows)){
  57     warning("Provide the number of rows (sometimes called ranges or beds) (num.rows=).")
  58   }
  59
  60   if(num.rows %% 3 != 0){
  61     warning("The MADII design requires that the number of rows be a multiple of 3.")
  62   }
  63
  64
  65   ## Develop other non-input function parameters
  66   if(is.null(entries)){
  67     entries <- as.matrix(paste("entry", 1:num.entries, sep="_"))
  68   }
  69
  70   if(is.null(num.entries)){
  71     entries <- as.matrix(entries) # If user input of entries was a list, will convert it to matrix
  72     num.entries <- nrow(entries)
  73   }
  74
  75
  76   ## This warning is dependent on the number of entries
  77   if(is.null(num.cols)==F){
  78     if(((num.cols / 5) * (num.rows / 3) + num.entries) > num.cols*num.rows){
  79       warning("The minimum number of plots has not been met by the given row/column dimensions.")
  80     }
  81   }
  82
  83   if(is.null(chk.names)){
  84     sec.chks <- as.character(2:(num.sec.chk+1)) ## Do need as seperate object for later on in function
  85     chk.names <- paste("chk", c(1,sec.chks), sep="") ## All generic check names
  86   }
  87
  88   if(is.null(num.sec.chk)){
  89     sec.chks <- chk.names[-1] ## The primary check must be listed first in the function input
  90     num.sec.chk <- length(sec.chks)
  91   }
  92
  93   blk.rows <- num.rows / 3 # This is the standard for the MADII design
  94
  95   ## If the number of columns is provided then it is straight forward, otherwise a set of algorithms will develop and optimize the paramaters
  96   if(is.null(num.cols)==F){
  97
  98     if(num.cols%%5 != 0){
  99       warning("The MADII design requires that the number of columns be a multiple of 5.")
 100     }
 101
 102     blk.cols <- num.cols / 5 # This is the standard for the MADII design
 103
 104     num.blks <- blk.rows * blk.cols
 105     exp.size <- num.blks * 15
 106     num.chks <- exp.size - num.entries
 107     #JL is there a rule that if a block has secondary checks, then it must have _all_ secondary checks?
 108     # That would not seem like a good rule.  Maybe it is required by some MADII analyses but it would not be
 109     # required by the moving average analysis or by an AR1 analysis, and my guess is those are better.
 110     num.sec.chk.blks <- floor((num.chks - num.blks) / num.sec.chk) # one primary check in each block
 111     # Number of checks to total number of plots: I think per.chks gets awfully high depending on the design.  I think there
 112     # should be no more than two check plots per block of 15.  It's just too much effort otherwise.
 113     per.chks <- (num.blks + (num.sec.chk.blks*num.sec.chk)) / (num.entries + num.blks + (num.sec.chk.blks*num.sec.chk))
 114     num.fill <- exp.size - (num.blks + num.entries + num.sec.chk.blks*num.sec.chk) # Fill lines are empty plots at the end of the experiment
 115     #JL do Fill plots actually get planted to something?  I am thinking that they would be.  If they are, it would seem better to
 116     #put secondary checks in them and distribute them randomly around the experiment
 117
 118     if(is.null(designID)==F){
 119       write.table(per.chks, paste(designID, "/", "%checks_in_", designID, ".txt", sep=""))
 120     }else{
 121       write.table(per.chks, "%checks.txt")
 122     }
 123
 124
 125   }else{
 126     ## If the number of columns is not specified, below algorithms will develop the necessary design
 127     ## Calculate starting (non-optimized paramaters)
 128     per.chks <- 0.10 #JL yes, this is a good number
 129
 130     ## Number of total checks in experiment (primary + secondary) ; calculated as percent of entries
 131     num.chks <- ceiling((per.chks * num.entries) / (1-per.chks))
 132     entries.plus.chks <- num.entries + num.chks
 133
 134     num.cols <- ceiling(entries.plus.chks / num.rows)
 135     #JL num.cols <- ceiling(num.cols / 5) * 5
 136     # so the statement above could just be "num.cols <- ceiling(entries.plus.chks / num.rows / 5) * 5"
 137     while(num.cols %% 5 != 0){
 138       num.cols <- num.cols + 1
 139     }
 140
 141     blk.cols <- num.cols / 5 # This is the standard for the MADII design
 142     num.blks <- blk.rows * blk.cols
 143     exp.size <- num.blks*15 # 15 plots per block
 144     num.sec.chk.blks <- ceiling((num.chks - num.blks) / num.sec.chk) # one primary check in each block
 145
 146     ## If the ratio of blk.cols to num.sec.chk.blks does not allow each blk.col to have a sec.chk blk in it, then optimize per.chks
 147     while((blk.cols > num.sec.chk.blks) & (num.sec.chk.blks <  num.blks)){
 148       per.chks <- per.chks + 0.0001
 149       num.chks <- ceiling((per.chks * num.entries) / (1-per.chks))
 150       entries.plus.chks <- num.entries + num.chks
 151
 152       num.cols <- ceiling(entries.plus.chks / num.rows)
 153
 154       while(num.cols %% 5 != 0){
 155         num.cols <- num.cols + 1
 156       }
 157
 158       blk.cols <- num.cols / 5 # This is the standard for the MADII design
 159       num.blks <- blk.rows * blk.cols
 160       exp.size <- num.blks*15 # 15 plots per block
 161       num.sec.chk.blks <- ceiling((num.chks - num.blks) / num.sec.chk) # one primary check in each block
 162
 163     }
 164
 165     num.fill <- num.blks*15 - (num.blks + num.entries + num.sec.chk.blks*num.sec.chk) # Fill lines are empty plots at the end of the experiment
 166
 167     ## Increase number of checks to minimize number of Fill plots
 168     while(num.fill >= num.sec.chk & (num.sec.chk.blks <  num.blks)){
 169       per.chks <- per.chks + 0.0001
 170       num.chks <- ceiling((per.chks * num.entries) / (1-per.chks))
 171       entries.plus.chks <- num.entries + num.chks
 172
 173       num.sec.chk.blks <- floor((num.chks - num.blks) / num.sec.chk) # one primary check in each block
 174
 175       num.fill <-num.blks*15 - (num.blks + num.entries + num.sec.chk.blks*num.sec.chk) # Fill lines are empty plots at the end of the experiment
 176
 177     }
 178
 179     if(is.null(designID)==F){
 180       write.table(per.chks, paste(designID, "/", "%checks_in_", designID, ".txt", sep=""))
 181     }else{
 182       write.table(per.chks, "%checks.txt")
 183     }
 184
 185   }
 186
 187
 188   #################################################################
 189   ############### Build Field File ################################
 190
 191   ## Put together field design columns; plot, row, col, blk, row.blk, col.blk
 192   fld.dgn <- as.data.frame(matrix(data=c((plot.start:(plot.start-1+exp.size)),
 193                                           rep(1:num.rows, each=num.cols), rep(c(1:num.cols, num.cols:1), length.out=exp.size),
 194                                           rep(NA, times=exp.size), rep(1:blk.rows, each=(exp.size/blk.rows)),
 195                                           rep(c(1:blk.cols, blk.cols:1), each=5, length.out=exp.size),
 196                                           rep(NA, times=2*exp.size)), nrow=exp.size, byrow=F))
 197
 198   colnames(fld.dgn) <- c("Plot", "Row", "Col", "Blk", "Row.Blk", "Col.Blk", "Line.Code", "Entry")
 199
 200   if(num.fill>0){
 201     fld.dgn[(exp.size-num.fill+1):exp.size, 7] <- "F"
 202   }
 203
 204   blk.list <- 1:num.blks
 205   for(b in 1:blk.rows){
 206     if((b %% 2 == 0)==T){
 207       blk.list[(1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols)] <- rev((1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols))
 208     } else{
 209       blk.list[(1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols)] <- ((1+blk.cols*(b-1)):((blk.cols*(b-1))+blk.cols))
 210     }
 211   }
 212
 213   ## Assign plots to blocks
 214   count <- 1
 215   for(b in 1:blk.rows){
 216     for(c in 1:blk.cols){
 217       blk <- blk.list[count]
 218       r1 <- (1+seq(0,1000,by=3))[b]
 219       r2 <- (seq(0,1000,by=3))[b+1]
 220
 221       c1 <- (1+seq(0,250000,by=5))[c]
 222       c2 <- (seq(0,250000,by=5))[c+1]
 223
 224       fld.dgn[which(fld.dgn$Row %in% r1:r2 & fld.dgn$Col %in% c1:c2), 4] <- blk
 225
 226       count <- count+1
 227     }
 228   }
 229
 230
 231   ## Selecting secondary blocks randomly, but ensuring that each blk.row and each blk.col is represented at least once
 232
 233   if((num.sec.chk.blks > blk.cols)==T){
 234     satisfied <- F
 235     col.blk.list <- c(1:blk.cols)
 236     row.blk.list <- c(1:blk.rows)
 237     while(satisfied != T){
 238       blk.cols.rep <- NULL
 239       blk.rows.rep <- NULL
 240       sample.kept <- c()
 241       new.length <- c()
 242       while((length(blk.cols.rep) < blk.cols)==T){
 243         length <- new.length
 244         sample <- sample(1:(num.blks), 1, replace=F)
 245         blk.cols.rep <- c(blk.cols.rep, unique(fld.dgn[which(fld.dgn$Blk==sample), 6]))
 246         blk.cols.rep <- unique(blk.cols.rep)
 247         new.length <- length(blk.cols.rep)
 248         if(new.length == 1){
 249           sample.kept <- c(sample.kept, sample)
 250         }else{
 251           if(new.length == (length+1))
 252             sample.kept <- c(sample.kept, sample)
 253         }
 254       }
 255
 256       for(s in sample.kept){
 257         blk.rows.rep <- c(blk.rows.rep, unique(fld.dgn[which(fld.dgn$Blk==s), 5]))
 258         blk.rows.rep <- unique(blk.rows.rep)
 259       }
 260
 261       if(all(col.blk.list %in% blk.cols.rep) & all(row.blk.list %in% blk.rows.rep)){
 262         satisfied <- T
 263       }
 264
 265     }
 266
 267     ## Need to bring "sample" back up to the number of secondary check blocks with some randomly chosen blocks
 268     cat((num.sec.chk.blks - length(sample.kept)), "\n")
 269     sample2 <- sample(setdiff(1:num.blks, sample.kept), (num.sec.chk.blks - length(sample.kept)))
 270     sample <- c(sample.kept, sample2)
 271
 272   }else{
 273     ## Selecting secondary blocks randomly, but ensuring that if sec.chk.blk < blk.cols then >1 sec.chk.block does not end up in a column
 274     satisfied <- F
 275     col.blk.list <- c(1:blk.cols)
 276     row.blk.list <- c(1:blk.rows)
 277     while(satisfied != T){
 278       blk.cols.rep <- NULL
 279       blk.rows.rep <- NULL
 280       #sample <- sample(1:(num.blks-(ceiling(num.fill/5))), num.sec.chk.blks, replace=F)
 281       # Use this line instead of below line if secondary checks are not wanted in partial plots due to Filler
 282       sample <- sample(1:(num.blks), num.sec.chk.blks, replace=F) # Using this line will allow secondary blocks to contain fill plots
 283       for(s in sample){
 284         blk.cols.rep <- c(blk.cols.rep, unique(fld.dgn[which(fld.dgn$Blk==s), 6]))
 285         blk.rows.rep <- c(blk.rows.rep, unique(fld.dgn[which(fld.dgn$Blk==s), 5]))
 286       }
 287       if(length(unique(blk.cols.rep))==length(blk.cols.rep) & all(row.blk.list %in% blk.rows.rep)){
 288         satisfied <- T
 289       }
 290     }
 291   }
 292
 293   ## Assign primary checks to field design
 294   for(b in 1:num.blks){
 295     blk <- fld.dgn[which(fld.dgn$Blk==b), ]
 296     row <- which(fld.dgn$Row == mean(blk$Row))
 297     col <- which(fld.dgn$Col == mean(blk$Col))
 298     fld.dgn[row[which(row  %in% col ==T)], 7] <- 1
 299   }
 300
 301   ## Assign secondary checks to field
 302   for(s in sample){
 303     sec.blk <- fld.dgn[which((fld.dgn$Blk==s) & is.na(fld.dgn$Line.Code)),][,1]
 304     sec.plots <- sec.blk[sample(1:length(sec.blk), num.sec.chk)]
 305     for(i in 1:length(sec.plots)){
 306       plot <- sec.plots[i]
 307       fld.dgn[which(fld.dgn$Plot==plot), 7] <- (i+1)
 308     }
 309   }
 310
 311
 312   fld.dgn[which(is.na(fld.dgn$Line.Code)), 7] <- 0
 313
 314   options(warn=-1) ## If no fill plots then the lines below will throw an error
 315   ## Assign entry names and check names
 316   fld.dgn[which(fld.dgn$Line.Code == 0), 8] <- entries[order(sample(1:length(entries), length(entries)))]
 317   if(num.fill>0){
 318     fld.dgn[which(fld.dgn$Line.Code == "F"), 8] <- "Fill"
 319   }
 320   fld.dgn[which(fld.dgn$Line.Code == "F"), 7] <- 0 # Change the fill lines back to experimental entries for purposes of MADIIadj
 321   options(warn=0) # Turn warnings back on
 322
 323
 324   for(c in 1:length(chk.names)){
 325     chk <- chk.names[c]
 326     fld.dgn[which(fld.dgn$Line.Code==c), 8] <- chk
 327   }
 328
 329   # Instead of re-working code, just re-order fld.dgn df before reading out
 330
 331
 332   fld.dgn.mod <- cbind(matrix(enviro, nrow=num.rows, ncol=1), fld.dgn) ; colnames(fld.dgn.mod)[c(1,8)] <- c("Enviro", "Check")
 333   fld.dgn.mod <- fld.dgn.mod[,c(1,2,9,8, 3:7)]
 334
 335   if(is.null(designID)==F){
 336     assign(designID, fld.dgn.mod, envir=.GlobalEnv)
 337     cat("\nYour MADII design is available in your working directory and R environment; named per ID provided.")
 338     write.csv(get(designID), paste(designID,"/",designID, ".csv", sep=""), row.names=F)
 339   } else{
 340     assign("YourMADIIdgn", fld.dgn.mod, envir=.GlobalEnv)
 341     cat("\nYour MADII design is available in your working directory and R environment; named 'YourMADIIdgn'.")
 342     write.csv(YourMADIIdgn, "YourMADIIdgn.csv", row.names=F)
 343   }
 344
 345   ########## Generating Visual Field Map ##########
 346
 347   ## First will generate a series of .csv files for excel
 348   ### Plots
 349   for.fld.plot <- fld.dgn[order(fld.dgn$Row, fld.dgn$Col),1:3]
 350   fld.plot <- apply(t(matrix(for.fld.plot$Plot, nrow=num.rows, ncol=num.cols, byrow=T)), 1, rev)
 351   # transposes matrix then flips all rows horizontally
 352   fld.plot <- cbind(matrix(num.rows:1, nrow=num.rows, ncol=1), fld.plot)
 353   fld.plot <- rbind(fld.plot, matrix(c(" ", 1:num.cols), nrow=1, ncol=ncol(fld.plot)))
 354
 355   ### Line codes
 356   for.line.code <- fld.dgn[order(fld.dgn$Row, fld.dgn$Col), c(2,3,7)]
 357   line.code <- apply(t(matrix(for.line.code$Line.Code, nrow=num.rows, ncol=num.cols, byrow=T)), 1, rev)
 358   #transposes matrix then flips all rows horizontally
 359   line.code <- cbind(matrix(num.rows:1, nrow=num.rows, ncol=1), line.code)
 360   line.code <- rbind(line.code, matrix(c(" ", 1:num.cols), nrow=1, ncol=ncol(line.code)))
 361
 362   ### Write out the design .csv files
 363   if(is.null(designID)==F){
 364     write.table(fld.plot, paste(designID,"/",designID, "_plot.lay.csv", sep=""), row.names=F, col.names=F, sep=",")
 365     write.table(line.code, paste(designID, "/", designID, "_line.codes.lay.csv", sep=""), row.names=F, col.names=F, sep=",")
 366   }else{
 367     write.table(fld.plot, paste(designID, "plot.lay.csv", sep="."), row.names=F, col.names=F, sep=",")
 368     write.table(line.code, paste(designID, "line.codes.lay.csv", sep="."), row.names=F, col.names=F, sep=",")
 369   }
 370
 371
 372   ###### Now a .png image ######
 373   options(warn=-1)
 374   plot.new()
 375
 376   if(is.null(designID)==F){
 377     dev.copy(png, paste(designID, "/", designID, "_plot.png", sep=""))
 378   }else{
 379     dev.copy(png, "YourMADIIdgn.plog.png")
 380   }
 381
 382
 383   field.plot <- matrix(nrow=num.rows, ncol=num.cols)
 384   field.NA.plot <- t(matrix(NA, nrow=num.rows, ncol=num.cols))
 385   for(r in 1:num.rows){
 386     for(c in 1:num.cols){
 387       field.plot[r,c] <- fld.dgn[which(fld.dgn$Row==r & fld.dgn$Col==c), 1]
 388     }
 389   }
 390   field.plot <- t(field.plot)
 391
 392   for(i in 1:(num.sec.chk+1)){
 393     chk.plots <- as.list(fld.dgn$Plot[which(fld.dgn$Line.Code==i)])
 394     field.NA.plot[which(field.plot %in% chk.plots)] <- i
 395   }
 396
 397   fill.plots <- as.list(fld.dgn$Plot[which(fld.dgn$Entry=="Fill")])
 398   field.NA.plot[which(field.plot %in% fill.plots)] <- (i+1)
 399
 400   colors <- c("red", "blue", "green", "orange", "gray", "brown", "yellow", "purple", "cyan")
 401   par(new=T,mar=c(4.5,2,.5,.5), srt=0, xpd=T, font=2, bg="white")
 402   s
 403   image(field.plot, add=F, col=0, xaxt='n', yaxt='n' ,xlab="Columns", ylab="Rows", line=0.5)
 404   image(field.NA.plot, add=T, col=colors[1:sum(length(chk.names),(num.fill/num.fill), na.rm=T)])
 405
 406   box(lwd=3)
 407
 408   if(num.fill>0){
 409     legend("bottom",legend=c(chk.names, "Fill"), bty="n", horiz=T, cex=.75, inset=c(0,-.125),
 410            xjust=0.5, text.col=colors[1:sum(length(chk.names),(num.fill/num.fill), na.rm=T)])}
 411   else{
 412     legend("bottom",legend=c(chk.names), bty="n", horiz=T, cex=.75, inset=c(0,-.125),
 413            xjust=0.5, text.col=colors[1:sum(length(chk.names),(num.fill/num.fill), na.rm=T)])
 414   }
 415
 416
 417   par(xpd=F)
 418   grid(ny=num.rows, nx=num.cols, col=1, lty=1 )
 419   grid(ny=blk.rows, nx=blk.cols, col=1, lty=1, lwd=3)
 420
 421   dev.off()
 422
 423   cat("\n\nA graphic of the design can also be found in the project's folder.\n ")
 424
 425   options(warn=0)
 426
 427 } # End of MADIIdgn loop
 428
 429 test <- MADIIdgn(num.entries=330, num.rows=9, num.cols=NULL, num.sec.chk=3, designID="tester1", annoy=T)
 430
 431 #test2 <- MADIIdgn(num.entries=100, num.rows=6, num.cols=5, num.sec.chk=3, designID="tester1", annoy=T)
 432
 433 #layoutParameters <- function(num.entries=330, num.rows=9, num.cols=NULL, num.sec.chk=3, designID="tester1", annoy=T){
 434 #  minBlks <- ceiling(num.entries / 14)
 435 #
 436 #}