| This post was kindly contributed by SAS and R - go there to comment and to read the full post. |
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.
| This post was kindly contributed by SAS and R - go there to comment and to read the full post. |