3

My question is whether any R packages out there provide function(s) that make it reasonably easy to randomize a standard experimental design that may involve crossed factors, nesting, and/or blocking.

For concreteness, show me specifically how to produce a new randomization of the Oats experiment provided as a dataset in the nlme package.

> data(Oats, package="nlme")
> summary(Oats)
 Block           Variety       nitro          yield      
 VI :12   Golden Rain:24   Min.   :0.00   Min.   : 53.0  
 V  :12   Marvellous :24   1st Qu.:0.15   1st Qu.: 86.0  
 III:12   Victory    :24   Median :0.30   Median :102.5  
 IV :12                    Mean   :0.30   Mean   :104.0  
 II :12                    3rd Qu.:0.45   3rd Qu.:121.2  
 I  :12                    Max.   :0.60   Max.   :174.0 

In that experiment, there are six blocks. Each block is divided into three plots which are randomly assigned to varieties (one plot per variety in each block, each block separately randomized). Each plot is subdivided into 4 subplots, and randomly assigned to four amounts of nitrogen (0, 0.2, 0.4, and 0.6), one subplot per nitro level randomized separately in each plot. In the dataset, the plots are identifiable as the combinations of Block and Variety. (The response variable is yield, so that is not really part of the treatment design.)

Second question: Given you can randomize Oats, can you use the same package to easily randomize other classical designs, e.g., a three-factor CRD, a nested design, a three-period crossover design, or a 5x5 Graeco-Latin square?

I actually already know how to do this using base functions in the R language; so I'm not particularly interested in seeing a programmatic answer. I want to know if any existing packages make this easy. I can identify a few packages that may help, e.g., randomizeR and randomizr, but a quick reading of the documentation of these still does not make it very obvious (to me) how to do this.

I have the makings of a general-purpose randomization package I developed a few years ago for my students, and I am trying to decide whether to develop it further for release on CRAN.

Russ Lenth
  • 5,922
  • 2
  • 13
  • 21
  • I see that there are "close" votes as this question is viewed as off-topic by some people. I guess I had thought it was about programming, but am I correct that I should put it on CrossValidated instead? I can take this down and move it. – Russ Lenth Jan 04 '17 at 22:57
  • Questions asking for code / packages are off topic on [stats.SE]. This question will be off topic there. It does seem to be a coding question to me, so I would think it would be on topic here, but it certainly would not be there. – gung - Reinstate Monica Jan 05 '17 at 03:42

1 Answers1

2

Here's how I'd do it with randomizr:

data(Oats, package="nlme")

# get the latest version from github      
install.packages("devtools")
devtools::install_github("DeclareDesign/randomizr")

library(randomizr)
Oats <- within(Oats,{
  Variety_new <- block_ra(block_var = Block, 
                          condition_names = c("Golden Rain", "Marvellous", "Victory"))
  nitro_new <- block_ra(block_var = paste0(Block, Variety_new), 
                        condition_names = c(0, 0.2, 0.4, 0.6))
  })

# Original Random Assignment
with(Oats, table(Block, Variety))
with(Oats, table(Block, nitro))
with(Oats, table(Block, nitro, Variety))

# New Random Assignment
with(Oats, table(Block, Variety_new))
with(Oats, table(Block, nitro_new))
with(Oats, table(Block, nitro_new, Variety_new))

The key is to realize that you need to block on both variety and Block to randomize the subplots into nitro conditions. (that's why we need the call to paste0).

EDIT

here is another way to do this (see comments)

library(randomizr) 
des <- rev(expand.grid(subplot=1:4, wholeplot=1:3, block=1:6)) 
des <- within(des,{ 
   plot_id <- paste0(block, "_", wholeplot) 
   Variety <- block_and_cluster_ra(
                block_var = block, 
                clust_var = plot_id, 
                condition_names = c("Golden Rain", "Marvellous", "Victory")) 
   nitro <- block_ra(block_var = plot_id, 
                     condition_names = c(0, 0.2, 0.4, 0.6)) 
}) 

with(des, table(Variety, block)) 
with(des, table(Variety, nitro)) 
with(des, table(Variety, nitro, block)) 
Alex Coppock
  • 2,122
  • 3
  • 15
  • 31
  • This looked really promising, but when I tried it, I got `` for all elements of the columns `Variety_new` and `nitro_new` (and zero counts for the second set of tables in your code). I installed from github as specified, and have package version 0.5.0. --- The latest commits must have broken something, because when I reverted to the one on CRAN (0.4.1), it worked. Thanks for answering, though it's a bit discouraging that it took 4 months before somebody came forward who understood the question. – Russ Lenth May 02 '17 at 19:19
  • Sorry! You may have caught us in the middle of an update. Version 0.6.0 now on github (and accepted by CRAN but it takes a couple of days for the binaries...) – Alex Coppock May 03 '17 at 15:15
  • No need to be sorry! Github is both a blessing and a curse -- a blessing in making the latest developments instantly available -- a curse in making the latest errors instantly available. – Russ Lenth May 03 '17 at 17:43
  • A follow-up, if you don't mind: This is the same design, but constructed differently. Starting with `des <- rev(expand.grid(subplot=1:4, wholeplot=1:3, block=1:6))`, I want to use the **randomizr** package to randomly assign `wholeplot`s to the levels of `Variety` within each `block`, and randomly assign the `subplot`s to the levels of `nitro` within each `block:wholeplot`. Is that straightforward? – Russ Lenth May 03 '17 at 17:47
  • 1
    Is this what you need? library(randomizr) des <- rev(expand.grid(subplot=1:4, wholeplot=1:3, block=1:6)) des <- within(des,{ plot_id <- paste0(block, "_", wholeplot) Variety <- block_and_cluster_ra(block_var = block, clust_var = plot_id, condition_names = c("Golden Rain", "Marvellous", "Victory")) nitro <- block_ra(block_var = plot_id, condition_names = c(0, 0.2, 0.4, 0.6)) }) with(des, table(Variety, block)) with(des, table(Variety, nitro)) with(des, table(Variety, nitro, block)) – Alex Coppock May 03 '17 at 20:59
  • Exactly. And it even works with the previous version of **randomizr**. Thanks. I like this way better as it preserves the overall block and plot structure. Now I'll try to understand this well enough to do things like this on my own. PS - I wonder if you might add this (formatted nicely) to your answer so future viewers can understand it better. – Russ Lenth May 03 '17 at 22:27