Example 2014.2: Block randomization

This week I had to block-randomize some units. This is ordinarily the sort of thing I would do in SAS, just because it would be faster for me. But I had already started work on the project R, using knitr/LaTeX to make a PDF, so it made sense to continue the work in R.

R
As is my standard practice now in both languages, I set thing up to make it easy to create a function later. I do this by creating variables with the ingredients to begin with, then call them as variables, rather than as values, in my code. In the example, I assume 40 assignments are required, with a block size of 6.
I generate the blocks themselves with the rep() function, calculating the number of blocks needed to ensure at least N items will be generated. Then I make a data frame with the block numbers and a random variate, as well as the original order of the envelopes. The only possibly confusing part of the sequence is the use of the order() function. What it returns is a vector of integer values with the row numbers of the original data set sorted by the listed variables. So the expression a1[order(a1$block,a1$rand),] translates to "from the a1 data frame, give me the rows ordered by sorting the rand variable within the block variable, and all columns." I assign the arms in a systematic way to the randomly ordered units, then resort them back into their original order.

seed=42
blocksize = 6
N = 40
set.seed(seed)
block = rep(1:ceiling(N/blocksize), each = blocksize)
a1 = data.frame(block, rand=runif(length(block)), envelope= 1: length(block))
a2 = a1[order(a1$block,a1$rand),]
a2$arm = rep(c("Arm 1", "Arm 2"),times = length(block)/2)
assign = a2[order(a2$envelope),]

> head(assign,12)
   block       rand envelope   arm
1      1 0.76450776        1 Arm 1
2      1 0.62361346        2 Arm 2
3      1 0.14844661        3 Arm 2
4      1 0.08026447        4 Arm 1
5      1 0.46406955        5 Arm 1
6      1 0.77936816        6 Arm 2
7      2 0.73352796        7 Arm 2
8      2 0.81723044        8 Arm 1
9      2 0.17016248        9 Arm 2
10     2 0.94472033       10 Arm 2
11     2 0.29362384       11 Arm 1
12     2 0.14907205       12 Arm 1

It's trivial to convert this to a function-- all I have to do is omit the lines where I assign values to the seed, sample size, and block size, and make the same names into parameters of the function.

blockrand = function(seed,blocksize,N){
  set.seed(seed)
  block = rep(1:ceiling(N/blocksize), each = blocksize)
  a1 = data.frame(block, rand=runif(length(block)), envelope= 1: length(block))
  a2 = a1[order(a1$block,a1$rand),]
  a2$arm = rep(c("Arm 1", "Arm 2"),times = length(block)/2)
  assign = a2[order(a2$envelope),]
  return(assign)
}

SAS
This job is also pretty simple in SAS. I use the do loop, twice, to produce the blocks and items (or units) within block, sssign the arm systematically, and generate the random variate which will provide the sort order within block. Then sort on the random order within block, and use the "Obs" (observation number) that's printed with the data as the envelope number.

%let N = 40;
%let blocksize = 6;
%let seed = 42;
data blocks;
call streaminit(&seed);
do block = 1 to ceil(&N/&blocksize);
  do item = 1 to &blocksize;
     if item le &blocksize/2 then arm="Arm 1";
    else arm="Arm 2";
     rand = rand('UNIFORM');
  output;
  end;
end;
run;

proc sort data = blocks; by block rand; run;

proc print data = blocks (obs = 12) obs="Envelope"; run;

Envelope    block    item     arm       rand

    1         1        3     Arm 1    0.13661
    2         1        1     Arm 1    0.51339
    3         1        5     Arm 2    0.72828
    4         1        2     Arm 1    0.74696
    5         1        4     Arm 2    0.75284
    6         1        6     Arm 2    0.90095
    7         2        2     Arm 1    0.04539
    8         2        6     Arm 2    0.15949
    9         2        4     Arm 2    0.21871
   10         2        1     Arm 1    0.66036
   11         2        5     Arm 2    0.85673
   12         2        3     Arm 1    0.98189

It's also fairly trivial to make this into a macro in SAS.

%macro blockrand(N, blocksize, seed); 
data blocks;
call streaminit(&seed);
do block = 1 to ceil(&N/&blocksize);
  do item = 1 to &blocksize;
     if item le &blocksize/2 then arm="Arm 1";
    else arm="Arm 2";
     rand = rand('UNIFORM');
  output;
  end;
end;
run;

proc sort data = blocks; by block rand; run;
%mend blockrand;

An unrelated note about aggregators: We love aggregators! Aggregators collect blogs that have similar coverage for the convenience of readers, and for blog authors they offer a way to reach new audiences. SAS and R is aggregated by R-bloggers, PROC-X, and statsblogs with our permission, and by at least 2 other aggregating services which have never contacted us. If you read this on an aggregator that does not credit the blogs it incorporates, please come visit us at SAS and R. We answer comments there and offer direct subscriptions if you like our content. In addition, no one is allowed to profit by this work under our license; if you see advertisements on this page, the aggregator is violating the terms by which we publish our work.

Example 9.18: Constructing the fastest relay team via enumeration

In competitive swimming, the medley relay is a team event in which four different swimmers each swim one of the four strokes: freestyle, breaststroke, backstroke, and butterfly. In general, every swimmer might be able swim any given stroke. How can we choose the fastest relay team? Here we solve this by enumerating all possible teams, though a more efficient routine likely exists.

Some example practice times can be seen on this Google Spreadsheet. Also, using the steps outlined here, the same spreadsheet is available as a CSV file here. (FTR, these are actual practice times for 100 yards for mostly 12-13 year-old swimmers; the names have been changed.)


SAS
We first read the data from the URL, using the technique outlined in section 1.1.6. Note that if you cut-and-paste this, you'll need to get the whole URL onto one line-- we break it up here for display only.

filename swimurl url 'https://docs.google.com/spreadsheet
/pub?key=0AvJKgZUzMYLYdE5xTHlEWkNUM3NoOHB1ZTJoTFMzUUE&
single=true&gid=0&output=csv';

proc import datafile=swimurl out=swim dbms=csv;
run;

Next, we use the point= option in nested set statements to generate a single data set with all the possible combinations of names and times. Meanwhile we change the names of the variables so they don't get overwritten in the next set statement. Note the use of the nobs option to find the number of rows in the data set.

data permute 
    (keep=free freetime fly flytime back backtime  breast breasttime);
set swim (rename = (swimmer=free freestyle=freetime)) nobs=nobs;
do i = 1 to nobs;
  set swim(rename = (swimmer=fly butterfly=flytime)) point=i;
  do j = 1 to nobs;
    set swim(rename = (swimmer=back backstroke=backtime)) point=j;
    do k = 1 to nobs;
      set swim(rename = (swimmer=breast breaststroke=breasttime)) 
        point = k;
      output;
 end;
  end;
end;
run;

The resulting data set has 12^4 rows, and includes rosters with the same swimmer swimming all four legs. In fact, a quick glance will show that Anna has the best time in each stroke, and thus the best "team" based on these practice times would use her for each stroke. This is against the rules, and also probably isn't reflective of performance in a race. We'll remove illegal line-ups using a where statement (section 1.5.1) and also calculate the total team time.

data prep;
set permute;
where free ne back and free ne breast and free ne fly and 
  back ne breast and back ne fly and breast ne fly;
time = sum(freetime, flytime, backtime, breasttime);
run;

The resulting data set has (12 permute 4) lines. To find the best team, we just sort by the total time and look at the first line. Here the first 10 lines (10 best teams) are shown.

proc sort data=prep; by time; run;
proc print data=prep (obs=10); run;

                                                        b
                                                        r
           f                             b              e
           r              f              a              a
           e              l              c      b       s
           e              y              k      r       t
   f       t              t      b       t      e       t      t
   r       i      f       i      a       i      a       i      i
   e       m      l       m      c       m      s       m      m
   e       e      y       e      k       e      t       e      e

 Kara    109.3  Dora    126.8  Lara    117.7  Anna    126.9  480.7
 Anna    102.8  Dora    126.8  Lara    117.7  Beth    134.6  481.9
 Kara    109.3  Anna    120.5  Lara    117.7  Beth    134.6  482.1
 Anna    102.8  Dora    126.8  Lara    117.7  Honora  136.4  483.7
 Kara    109.3  Jane    129.8  Lara    117.7  Anna    126.9  483.7
 Kara    109.3  Anna    120.5  Lara    117.7  Dora    136.4  483.9
 Kara    109.3  Anna    120.5  Lara    117.7  Honora  136.4  483.9
 Kara    109.3  Lara    123.1  Jane    124.7  Anna    126.9  484.0
 Anna    102.8  Dora    126.8  Lara    117.7  Inez    136.8  484.1
 Carrie  112.7  Dora    126.8  Lara    117.7  Anna    126.9  484.1

The best time shaves a whole second off the predicted time using the second-best team.

R
Since published Google Spreadsheets use https rather than http, we use the RCurl package and its getURL() function. (Note that if you cut-and-paste this, you'll need to get the whole URL onto one line-- we break it up here for display only.) Then we can read the data with read.csv() and textConnection().

library(RCurl)
swim = getURL("https://docs.google.com/spreadsheet
/pub?key=0AvJKgZUzMYLYdE5xTHlEWkNUM3NoOHB1ZTJoTFMzUUE&
single=true&gid=0&output=csv")

swim2=read.csv(textConnection(swim))

To make an object with the combinations of names, we use the expand.grid() function highlighted in Example 7.22, providing as arguments the swimmers names four times. As in the SAS example, the result has has 12^4 rows. The combn() function might be a better fit here, but was more difficult to use.

test2 = expand.grid(swim2$Swimmer,swim2$Swimmer, swim2$Swimmer, swim2$Swimmer)  

It'll be useful to assign these copies of the names to each of the strokes. We'll do that with the rename() function available in the reshape package. (This approach is mentioned in section 1.3.4.). Then we can remove the rows where the same name appears twice using some logic. The logic is nested in the with() function to save some keystrokes and is generally preferable to attach()ing the test2 object.

library(reshape)
test2 = rename(test2, c("Var1" = "free", "Var2" = "fly", 
  "Var3" = "back", "Var4" = "breast"))
test3 = with(test2, test2[(free != breast) & (free != fly) 
  & (free != back) & (breast != fly) & (breast != back) 
  & (fly != back) ,])

Finally, we can use the which.min() function to pick the best team.

> bestteam =   
+ test3[which.min(swim2$Freestyle[test3$free]+swim2$Breaststroke[test3$breast] +
+ swim2$Butterfly[test3$fly] + swim2$Backstroke[test3$back]),]
>bestteam
     free  fly back breast
1631 Kara Dora Lara   Anna

For new users of R, this may look very peculiar-- it uses elegant but powerful features of the R language that may be challenging for new users to grasp. Essentially, in swim2$Freestyle[test3$free] we say: from the "freestyle" times in the swim2 object, take the time from the row that has the name in a row of "free" names in the test3 object. The which.min() function replicates this request for every row in the test3 object (which has all of the permutations) in it, returning the row number with that minimum sum. The outer test3[rows,columns] syntax grabs the values in this row. (The number 1631 is the row number, for some reason showing the row in the test2 object created by expand.grid().)

Now, we might also want the actual times associated with the best team. We can find them by calling the correct rows (names from the best team) and columns (stroke associated with that name) from the original data set.

> times = c(swim2[swim2$Swimmer == bestteam$free,2], 
+      swim2[swim2$Swimmer == bestteam$fly,3], 
+      swim2[swim2$Swimmer == bestteam$back,4], 
+      swim2[swim2$Swimmer == bestteam$breast,5])
> times
[1] 109.3 126.8 117.7 126.9

If instead, one wanted to list the times in order, one approach would be to add columns to the test3 object with the time for each stroke, calculate their sum, and sort on the sum.

Example 7.19: find the closest pair of observations

Suppose we need to find the closest pair of observations on some variable x. For example, we might be concerned that some data had been accidentally duplicated. We return the ID's of the two closest observations, and their distance from each other. In both languages, we'll first create the data, then sort it, recognizing that the smallest difference must come between two adjacent observations, once they're sorted.

SAS

First, we'll generate the data (see section 1.11.1) and and sort it (section 1.5.6).

data ds;
do id = 1 to 10;
  x = normal(0);
  output;
  end;
run;

proc sort data=ds; by x; run;

Then we'll find the closest pairing and find their IDs. There's relatively complicated stuff going on here, so we annotate within the code rather than try to explain in text.

data smallest;
set ds end=eof;
retain smalldist last lastid smallid largeid;  
   /* variables we keep between processing observations 
      in the ds dataset */
if _n_ gt 1 then do;  /* no diff for the smallest obs */
   smalldist = min(smalldist, x - last);  
                 /* is it smaller now than last time? */
   if smalldist = x - last then do;       /* if so */
      smallid = lastid;                    
                 /* the current smallest id is the one 
                                    from the last row */
      largeid = id;                        
                 /* and the current largest id is */
   end;                                 /* the one in this row */
end;
last = x;           /* if this row turns out to be the  */
lastid = id;       /* smaller end of the smallest pair  */
                     /* we'll need to know x and the id */
if eof then output;   /* if we're at the last row, stop */
run;

The value of smalldist is missing for the first execution of the min function, which returns the minimum of the non-missing observations as discussed in section 1.4.14. Note that the end option to the set statement makes a temporary variable which = 1 for the last row of the data set. An output statement is the implicit last line of all data steps. When it is explicit, lines are output only when the condition is true.

The result is a one-line data set with the two IDs as smallid and largeid and their distance as smalldist.

Here is an example of the generated input data and the resulting smallest dataset.

proc print data = ds; run;

Obs    id        x

  1     8    -2.14041
  2     2    -1.17532
  3    10    -0.80683
  4     4    -0.77437
  5     9    -0.06806
  6     5     0.14425
  7     6     0.40988
  8     7     0.66603
  9     1     0.68565
 10     3     1.86831

proc print data = smallest;
var smalldist smallid largeid; 
run;

Obs    smalldist    smallid    largeid

 1      0.019622       7          1

R

First, generate data (1.11.1). Next, sort the observations (1.5.6), then find the adjacent differences.

id <- 1:10
x <- rnorm(10)
sortx <- x[order(x)]
oneless <- sortx[2:length(x)]
diff <- oneless - sortx[1:length(x)-1]

Again, the next step is complex, but the use of indexing makes intra-code comments awkward, so we attempt to explain it here. The ID of the smaller of the two observations with the smallest distance is the value in the ID vector that's in the place where x equals the value of the sorted x vector that's in the same place as the smallest difference. The larger ID can be found the same way, using the shifted vector.

smallid <- id[x == sortx[which.min(diff)]]
largeid <- id[x == oneless[which.min(diff)]]
smalldist <-  min(diff)

> options(digits=2)
> x
 [1]  1.412 -0.518  1.298  0.351  2.123 -1.388  
 [7]  0.431  1.268  0.658 -0.014
> sortx
 [1] -1.388 -0.518 -0.014  0.351  0.431  0.658  
 [7]  1.268  1.298  1.412  2.123
> oneless
 [1] -0.518 -0.014  0.351  0.431  0.658  1.268  
 [8]  1.298  1.412  2.123
> diff
 [1] 0.87 0.50 0.36 0.08 0.23 0.61 0.03 0.11 0.71
> smallid
 [1] 8
> largeid
 [1] 3
> smalldist
 [1] 0.03

Example 7.11: Plot an empirical cumulative distribution function from scratch

In example 7.8, we used built-in functions to produce an empirical CDF plot. But the empirical cumulative distribution function (CDF) is simple to calculate directly, and it might be useful to have more control over its appearance than is afforded by the direct method employed in example 7.8. In a later post, we'll use that control to improve upon example 7.8.

SAS

In SAS, we first read the data from a locally stored Stata data set. Then, we calculate the empirical CDF for pcs. We do this by first sorting (1.5.6) on pcs. Then the order statistic of the observation divided by the total number of obervations is the empirical probability that X ≤ x. We find the number of observations using the nobs data set option (A.6.1) and the observation number using the _n_ implied variable (1.4.15). Then we can plot the empirical CDF using proc gplot (section 5.1.1) to make a scatterplot and a symbol statement with the j interpolation (as in section 1.13.5) to connect the points.

proc import datafile="C:\book\help.dta"
   out=help dbms=dta;
run;

proc sort data=help;
   by pcs;
run;

data help_a;
set help nobs=totalobs;
   ecdf_pcs = _n_ / totalobs;
run;

symbol1 i=j v=none c=blue;
proc gplot data=help_a;
   plot ecdf_pcs * pcs;
run;
quit;

R
Now, R. First, make the ECDF values in a similar fashion to SAS, then make an empty frame (5.1.1), then fill the frame with lines connecting the calculated values (5.2.1).

library(foreign)
ds <- read.dta
  ("http://www.math.smith.edu/sasr/datasets/help.dta")
attach(ds)
sortpcs <- sort(pcs)
ecdfpcs <- (1:length(sortpcs))/length(sortpcs)
plot(sortpcs, ecdfpcs, type="n")
lines(sortpcs, ecdfpcs)

Example 7.10: Get data from R into SAS

In our previous entry, we described how to generate a dataset from SAS that could be used for analyses in R. Alternatively, someone primarily using R might want to test the new ”statistical graphics” procedures available starting with SAS 9.2. Here we show how one might create a long data set in R, export it, read it in SAS, and generate a plot similar to that shown in example 7.9.

R

In R, we use the reshape() function (section 1.5.3) to make the ”long” data set. We make a data frame (section B.4.5) from the data we’ll need in SAS. We check the data, using the order() function (section 1.5.6) to organize the new data frame by subject instead of visit number. Finally, we write out the data set in dBase format (section 1.2.2).

> ds <- read.csv
   ("http://www.math.smith.edu/sasr/datasets/helpmiss.csv")
> attach(ds)
> long <- reshape(ds, idvar="id", 
    varying=list(c("cesd1","cesd2","cesd3","cesd4")),
    v.names="cesdtv", timevar="visit", direction="long")
> attach(long)
> tosas <- data.frame(id, visit, cesdtv)
> head(tosas[order(id, visit),])

     id visit cesdtv
1     1     1      7
471   1     2     NA
941   1     3      8
1411  1     4      5
2     2     1     11
472   2     2     NA

> library(foreign)
> write.dbf(tosas,"c:\\book\\tosas.dbf")

SAS
In SAS, we read in the data from the dBase format file (section 1.1.5), proc sort (section 1.5.6) it in place, and check the data using proc print(section 1.2.4).

proc import datafile="C:\book\tosas.dbf" 
   out=fromr dbms=dbf;
run;

proc sort data=fromr; by id visit; run;

proc print data=fromr (obs=5); run;

 Obs        id    visit    cesdtv

   1         1        1         7
   2         1        2         .
   3         1        3         8
   4         1        4         5
   5         2        1        11

Finally, we generate the desired plot with proc sgpanel (section 5.1.11), using the where statement (section A.6.3) to select the first 20 subjects.

proc sgpanel data=fromr;
   where id le 20;
   panelby id / rows=4 columns=5;
   scatter x=visit y=cesdtv;
run;