<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/balancing-on-multiple-factors-when-the-sample-is-too-small-to-stratify/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.rdatagen.net/post/balancing-when-sample-is-too-small-to-stratify/"> ouR data generation</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Ideally, a study that uses randomization provides a balance of characteristics that might be associated with the outcome being studied. This way, we can be more confident that any differences in outcomes between the groups are due to the group assignments and not to differences in characteristics. Unfortunately, randomization does not <em>guarantee</em> balance, especially with smaller sample sizes. If we want to be certain that groups are balanced with respect to a particular characteristic, we need to do something like stratified randomization.</p>
<p>When the sample size is small and we want to guarantee balance across <em>multiple</em> characteristics, the task is a bit more challenging. Say we have 20 schools that we are randomizing to two groups, 10 in each, and want to make sure the groups are balanced with respect to 4 characteristics: language, poverty, location, and size. Simple stratification may not work so well. If we assume that these four characteristics are binary (e.g. either “yes” or “no”), there are 16 possible combinations. One or more of these combinations could easily be represented by a single school – so it would be impossible to randomize within each of the 16 combinations. What to do?</p>
<p>One possible approach is to generate all possible randomization schemes of the 20 schools, and keep only those schemes that are balanced with respect to the four characteristics. Once we have a list of acceptable randomization schemes, we can just pick one of <em>those</em> at random. (Of course, it is preferable if each school has close to a 50% chance of being assigned to either intervention group.)</p>
<div>
<h2>Simulate school-level data</h2>
<p>To start, we generate data for our 20 hypothetical schools using <code>simstudy</code> functions:</p>
<pre><code>library(simstudy)
set.seed(125)

# define data characteristics for schools
ddef &lt;- defData(varname = &quot;language&quot;, formula = .3, dist = &quot;binary&quot;)
ddef &lt;- defData(ddef, &quot;poverty&quot;, formula = .2, dist = &quot;binary&quot;)
ddef &lt;- defData(ddef, &quot;location&quot;, formula = .5, dist = &quot;binary&quot;)
ddef &lt;- defData(ddef, &quot;size&quot;, formula = .5, dist = &quot;binary&quot;)
ddef</code></pre>
<pre><code>##     varname formula variance   dist     link
## 1: language     0.3        0 binary identity
## 2:  poverty     0.2        0 binary identity
## 3: location     0.5        0 binary identity
## 4:     size     0.5        0 binary identity</code></pre>
<pre><code># generate schools
dt &lt;- genData(20, ddef)

# number of schools in each combination
dt[, .N, keyby = .(language,poverty,location,size)]</code></pre>
<pre><code>##    language poverty location size N
## 1:        0       0        0    1 5
## 2:        0       0        1    0 1
## 3:        0       0        1    1 5
## 4:        0       1        0    0 1
## 5:        0       1        1    0 2
## 6:        1       0        0    0 2
## 7:        1       0        0    1 1
## 8:        1       0        1    0 1
## 9:        1       0        1    1 2</code></pre>
<p>In this case, we have nine different combinations of the four characteristics, four of which include only a single school (rows 2, 4, 7, and 8). Stratification wouldn’t work necessarily work here if our goal was balance across all four characteristics.</p>
</div>
<div>
<h2>Create randomization scenarios to assess for balance</h2>
<p>Ideally, we would generate all possible randomization combinations and check them all for balance. If the number of total units (e.g. schools) is small, this does not pose a challenge (e.g. if N=4, then we only have six possible randomization schemes: TTCC, TCTC, TCCT, CTTC, CTCT, CCTT). However, with N=20, then there are 184,756 possible randomization schemes. Depending on the efficiency of the algorithm, it may be impractical to evaluate all the schemes. So, an alternative is to sample a subset of the schemes and evaluate those. For illustration purposes (so that you can understand what I am doing), I am using some very inefficient <code>R</code> code (using a loops). As a result, I cannot evaluate all possible schemes in a reasonable period of time to get this post out; I decided to sample instead to evaluate 1000 possible randomizations. (At the end of this post, I show results using much more efficient code that uses data.table and Rcpp code much more effectively – so that we can quickly evaluate millions of randomization schemes.)</p>
<p>To start, I create all combinations of randomization schemes:</p>
<pre><code>totalSchools = 20
rxSchools = 10

xRx &lt;- t(combn(totalSchools, rxSchools)) 

# show 5 randomly sampled combinations

sampleRows &lt;- sample(nrow(xRx), 5, replace = FALSE)
xRx[sampleRows,]</code></pre>
<pre><code>##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    2    3    5    6    7    8   10   12   14    19
## [2,]    5    6    7    8   13   14   15   16   17    18
## [3,]    1    3    4    5    7    9   12   15   17    20
## [4,]    2    3    4    5    9   11   14   15   19    20
## [5,]    3    5    6    7    8   10   11   12   15    16</code></pre>
<p>Below is a function (which I chose to do in Rcpp) that converts the <code>xRx</code> matrix of school ids to a 20-column matrix of 1’s and 0’s indicating whether or not a school is randomized to the intervention in a particular scenario:</p>
<pre><code>#include &lt;RcppArmadillo.h>

// [[Rcpp::depends(RcppArmadillo)]]

using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]

NumericMatrix convert01(NumericMatrix xmat, int tcols) {
  
  int xrows = xmat.nrow();
  int xcols = xmat.ncol();
  
  NumericMatrix pmat(xrows, tcols);
  
  for (int i=0; i &lt; xrows; i++) {
    for (int j=0; j &lt; xcols; j++)  {
      pmat(i, xmat(i,j) - 1) = 1; 
    }
  } 
  return(pmat);
}</code></pre>
<pre><code>x01 &lt;- convert01(xRx, totalSchools)

# show some rows

x01[sampleRows,]</code></pre>
<pre><code>##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,]    0    1    1    0    1    1    1    1    0     1     0     1
## [2,]    0    0    0    0    1    1    1    1    0     0     0     0
## [3,]    1    0    1    1    1    0    1    0    1     0     0     1
## [4,]    0    1    1    1    1    0    0    0    1     0     1     0
## [5,]    0    0    1    0    1    1    1    1    0     1     1     1
##      [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,]     0     1     0     0     0     0     1     0
## [2,]     1     1     1     1     1     1     0     0
## [3,]     0     0     1     0     1     0     0     1
## [4,]     0     1     1     0     0     0     1     1
## [5,]     0     0     1     1     0     0     0     0</code></pre>
<p>Because the evaluation code is so inefficient, I draw 1,000 rows at random from this “intervention” matrix <code>x01</code> (after converting it to a data.table).</p>
<pre><code># convert matrix to data.table
d01 &lt;- data.table(x01)
d01[, id := .I]

ids &lt;- sample(nrow(d01), 1000, replace = FALSE)
sampleD01 &lt;- d01[id %in% ids]</code></pre>
<p>Now we are ready to evaluate each of the 1,000 schemes. As I mentioned before, this approach is highly inefficient as the algorithm requires us to literally loop through each each combination to find the balanced ones. I have sacrificed efficiency and speed for clarity of code (I hope).</p>
<pre><code>for (i in 1:1000) {
  
  dt[, grp:= t(sampleD01[i,1:20])]
  
  dx &lt;- dt[ , .N, keyby = .(language, grp)]
  dc &lt;- dcast(dx, language ~ grp, fill = 0, value.var = &quot;N&quot; )
  dc[, diff := abs(`1` - `0`)]
  
  # we declare a scheme balanced if counts differ by
  # no more than 1 school
  
  sampleD01[i, language := (sum(dc[, diff > 1]) == 0)]
  
  dx &lt;- dt[ , .N, keyby = .(poverty, grp)]
  dc &lt;- dcast(dx, poverty ~ grp, fill = 0, value.var = &quot;N&quot; )
  dc[, diff := abs(`1` - `0`)]
  
  sampleD01[i, poverty := (sum(dc[, diff > 1]) == 0)]
  
  dx &lt;- dt[ , .N, keyby = .(location, grp)]
  dc &lt;- dcast(dx, location ~ grp, fill = 0, value.var = &quot;N&quot; )
  dc[, diff := abs(`1` - `0`)]
  
  sampleD01[i, location := (sum(dc[, diff > 1]) == 0)]
  
  dx &lt;- dt[ , .N, keyby = .(size, grp)]
  dc &lt;- dcast(dx, size ~ grp, fill = 0, value.var = &quot;N&quot; )
  dc[, diff := abs(`1` - `0`)]
  
  sampleD01[i, size := (sum(dc[, diff > 1]) == 0)]
  
}</code></pre>
<p>The final determination of balance is made if a scheme is balanced across all four characteristics. In this case, 136 of the 1,000 schemes were balanced based on this criterion:</p>
<pre><code>sampleD01[, balanced := all(language, poverty, location, size), keyby = id]

# proportion of sampled combinations that are balanced ...

sampleD01[,mean(balanced)]</code></pre>
<pre><code>## [1] 0.136</code></pre>
<p>And let’s inspect the actual balance of two randomly selected schemes – one which is balanced, and one which is not:</p>
<pre><code>sTrue &lt;- sampleD01[balanced == TRUE]
sFalse &lt;- sampleD01[balanced == FALSE]</code></pre>
<div>
<h3>A balanced scheme</h3>
<pre><code>dtAssigned &lt;- copy(dt)
dtAssigned[, group := as.vector(t(sTrue[sample(.N, 1), 1:20]))]

dtAssigned[, .N, keyby=.(language, group)]</code></pre>
<pre><code>##    language group N
## 1:        0     0 7
## 2:        0     1 7
## 3:        1     0 3
## 4:        1     1 3</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(poverty, group)]</code></pre>
<pre><code>##    poverty group N
## 1:       0     0 9
## 2:       0     1 8
## 3:       1     0 1
## 4:       1     1 2</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(location, group)]</code></pre>
<pre><code>##    location group N
## 1:        0     0 4
## 2:        0     1 5
## 3:        1     0 6
## 4:        1     1 5</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(size, group)]</code></pre>
<pre><code>##    size group N
## 1:    0     0 3
## 2:    0     1 4
## 3:    1     0 7
## 4:    1     1 6</code></pre>
</div>
<div>
<h3>An unbalanced scheme</h3>
<p>In this case, language and location are imbalanced, though size and poverty are fine.</p>
<pre><code>dtAssigned &lt;- copy(dt)
dtAssigned[, group := as.vector(t(sFalse[sample(.N, 1), 1:20]))]

dtAssigned[, .N, keyby=.(language, group)]</code></pre>
<pre><code>##    language group N
## 1:        0     0 8
## 2:        0     1 6
## 3:        1     0 2
## 4:        1     1 4</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(poverty, group)]</code></pre>
<pre><code>##    poverty group N
## 1:       0     0 8
## 2:       0     1 9
## 3:       1     0 2
## 4:       1     1 1</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(location, group)]</code></pre>
<pre><code>##    location group N
## 1:        0     0 3
## 2:        0     1 6
## 3:        1     0 7
## 4:        1     1 4</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(size, group)]</code></pre>
<pre><code>##    size group N
## 1:    0     0 4
## 2:    0     1 3
## 3:    1     0 6
## 4:    1     1 7</code></pre>
</div>
</div>
<div>
<h2>Fast implementation with data.table and Rcpp</h2>
<p>As I alluded to before, if we want to implement this in the real world, it would be preferable to use code that does not bog down when we want to search 100,000+ possible randomization schemes. I have written a set of <code>R</code> and <code>Rcpp</code> functions the facilitate this. (Code is available <a rel="nofollow" href="https://github.com/kgoldfeld/RDataGenBlog/tree/master/static/img/post-balance" target="_blank">here</a>.)</p>
<pre><code># generate all possible schemes

xperm &lt;- xPerms(totalSchools, rxSchools, N=NULL) 

nrow(xperm)</code></pre>
<pre><code>## [1] 184756</code></pre>
<pre><code>xperm[sample(nrow(xperm), 5, replace = FALSE)]</code></pre>
<pre><code>##    V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 V14 V15 V16 V17 V18 V19
## 1:  0  1  1  1  1  0  1  1  0   0   0   1   0   1   0   0   0   1   0
## 2:  1  1  0  1  1  1  1  1  0   1   0   0   0   0   1   0   1   0   0
## 3:  1  0  1  0  0  1  1  1  1   1   1   0   1   0   0   1   0   0   0
## 4:  1  1  1  0  0  1  0  1  0   1   1   1   1   0   0   1   0   0   0
## 5:  1  1  0  0  1  0  0  1  1   1   1   0   1   0   0   0   1   1   0
##    V20    id
## 1:   1 94784
## 2:   0 19535
## 3:   0 61644
## 4:   0 14633
## 5:   0 35651</code></pre>
<pre><code># prepare data for evaluation

dtMat &lt;- as.matrix(dt[,-1])
cc &lt;- parse(text=attr(xperm, &quot;varlist&quot;))
cc</code></pre>
<pre><code>## expression(c(V1, V2, V3, V4, V5, V6, V7, V8, V9, V10, V11, V12, 
##     V13, V14, V15, V16, V17, V18, V19, V20))</code></pre>
<pre><code># evaluate each combination

sF &lt;-  xperm[, cppChk(eval(cc), dtMat), keyby = id]
sF[sample(nrow(sF), 5, replace = FALSE)]</code></pre>
<pre><code>##        id    V1
## 1:  15924 FALSE
## 2:  68284 FALSE
## 3: 149360 FALSE
## 4:  62924 FALSE
## 5:  14009  TRUE</code></pre>
<pre><code># keep only the balanced schemes

sFinal &lt;- xperm[sF$V1]
nrow(sFinal)</code></pre>
<pre><code>## [1] 7742</code></pre>
<pre><code># randomize from the balanced schemes

selectRow &lt;- sample(nrow(sFinal), 1)

# check balance of randomized scheme

dtAssigned &lt;- copy(dt)
dtAssigned[, group := as.vector(t(sFinal[selectRow, -&quot;id&quot;]))]

dtAssigned[, .N, keyby=.(language, group)]</code></pre>
<pre><code>##    language group N
## 1:        0     0 7
## 2:        0     1 7
## 3:        1     0 3
## 4:        1     1 3</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(poverty, group)]</code></pre>
<pre><code>##    poverty group N
## 1:       0     0 9
## 2:       0     1 8
## 3:       1     0 1
## 4:       1     1 2</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(location, group)]</code></pre>
<pre><code>##    location group N
## 1:        0     0 5
## 2:        0     1 4
## 3:        1     0 5
## 4:        1     1 6</code></pre>
<pre><code>dtAssigned[, .N, keyby=.(size, group)]</code></pre>
<pre><code>##    size group N
## 1:    0     0 3
## 2:    0     1 4
## 3:    1     0 7
## 4:    1     1 6</code></pre>
</div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://www.rdatagen.net/post/balancing-when-sample-is-too-small-to-stratify/"> ouR data generation</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/smrLcBQ5CzE" height="1"></a></div>