Here is the post:

I almost exclusively follow the R tag–it’s my current “playground”–and I have been noticing many questions phrased along the way of “what is an elegant way to do XYZ” or answers in the form of “this may not be the most elegant solution, but it does what you want.” (I know many is subjective, but search the R tag for elegant, and you will find at least one question every other day or so with the word in it.)

My “question” is, what exactly is elegant code?

Is it:

• Code that is easily readable?
• Code that does everything in one line (I see a lot of congratulations for these so-called “one-liners”)?
• Code that is faster?
• Something totally different?

I’ve seen in some documents (for example A (Not So) Short Introduction to S4, Section 1.6) what I feel are much more understandable terms for describing code: good, clean, bad, dangerous, dirty…. As I write this, @Andrie is schooling (politely–”Summer of Love” fashion) a new user on how to “clean your code and make life easier in the long term” (emphasis added). And, what comes to mind when I try to write code are things like consistency, readability, efficiency, not producing unexpected output, not breaking other functions. These descriptions for code seem to fit with such terms, but not necessarily with what I would call elegant.

I’m also often confused when an answerer posts what I consider an excellent answer, perhaps incrementally building towards the desired output, but modestly states it is not elegant. If the answer does exactly what is expected to do, isn’t that the whole point of the code?

Some further background

My initial questions regarding this started when I read a line in a question by @AriB.Friedman about making base R reshape easier to use that read:

It is a truth universally acknowledged that R’s base reshape command is speedy and powerful but has miserable syntax.

The word elegant was not used here, but this is an example of where I might use it. Ari was trying to make an incredibly useful function easier to use by making the user input more consistent, much like some of Hadley Wickham‘s packages are wrappers around base R functions, often combining common steps that a user might go through. The result of these new functions is not always more efficient, but the code that is written is often more easy to decipher, and thus provide what I might call a more elegant approach.

(Note: I fully expect that this question would be closed as “off topic”, “not a real question”, or “not constructive”, since it is not directly programming related; but while it stays open, hopefully there are some helpful answers.)

Don't make me squirt my milk bottle at you!

I also had a few alternatives–and I’m still not sure which one is my favorite.

Originally, the slogans were “we two, ours two”, or at least that was the catchy English version–regional languages usually had a slogan more along the lines of “one family, two children”. And, the change to the new slogan led to at least one humorous math discussion with an auto driver who commented that, “Earlier, it was ‘we two, ours two’; now, it is ‘we two, ours one’. What’s next? ‘We two, ours half?’”

Anyway, keen observers might have noticed the following new addition to selected trucks:

Get the function at http://news.mrdwab.com/stratified. Usage notes in the head of the function.

Here’s the function:

stratified = function(df, id, group, size, seed="NULL", ...) {
  #  USE: * Specify your data frame, ID variable (as column number),
  #   and grouping variable (as column number) as the first three
  #   arguments.
  #    * Decide on your sample size. For a sample proportional to
  #   the population, enter "size" as a decimal. For an equal number
  #   of samples from each group, enter "size" as a whole number.
  #    * Decide on if you want to use a seed or not. If not, leave
  #   blank or type "NULL" (with quotes).
  #
  #  Example 1: To sample 10% of each group from a data frame named
  #    "z", where the ID variable is the first variable, the
  #    grouping variable is the fourth variable, and the desired
  #    seed is "1", use:
  #
  #     > stratified(z, 1, 4, .1, 1)
  #
  #  Example 2: Run the same sample as above but without a seed:
  #
  #     > stratified(z, 1, 4, .1)
  #
  #  Example 3: To sample 5 from each group from a data frame named
  #    "z", where the ID variable is the first variable, the
  #    grouping variable is the third variable, and the desired
  #    seed is 2, use:
  #
  #     > stratified(z, 1, 3, 5, 2)
  #
  #  NOTE: Not tested on datasets with LOTS of groups or with HUGE
  #  differences in group sizes. Probably INCREDIBLY inefficient.
  k = unstack(data.frame(as.vector(df[id]), as.vector(df[group])))
  l = length(k)
  results = vector("list", l)
  if (seed == "NULL" & size < 1) {
    for (i in 1:length(k)) {
      N = k[[i]]
      n = round(length(N)*size)
      results[[i]] = list(sample(N, n, ...))
    }
  } else if (seed == "NULL" & size >= 1) {
    for (i in 1:length(k)) {
      N = k[[i]]
      results[[i]] = list(sample(N, size, ...))
    }
  } else if (size < 1) {
    for (i in 1:length(k)) {
      set.seed(seed)
      N = k[[i]]
      n = round(length(N)*size)
      results[[i]] = list(sample(N, n, ...))
    }
  } else if (size >= 1) {
    for (i in 1:length(k)) {
      set.seed(seed)
      N = k[[i]]
      results[[i]] = list(sample(N, size, ...))
    }
  }
  z = data.frame(c(unlist(results)))
  names(z) = names(df[id])
  w = merge(df, z)
  w[order(w[group]), ]
}

And here are some examples of the function in action:

source("http://news.mrdwab.com/stratified")
# Make up some data
dat = data.frame(A = 1:100,
                 B = sample(c("AA", "BB", "CC", "DD", "EE"), 100,
                            replace=T),
                 C = rnorm(100),
                 D = abs(round(rnorm(100), digits=1)),
                 E = sample(c("CA", "NY", "TX"), 100, replace=T))
# view the first few rows
head(dat)
# A  B     C   D  E
# 1 1 CC -0.07870439 0.6 NY
# 2 2 CC -0.65048634 0.3 TX
# 3 3 EE  1.02703616 1.3 NY
# 4 4 BB -1.08696775 0.4 TX
# 5 5 CC  0.56741795 0.2 CA
# 6 6 AA -0.46448941 0.5 TX
# Sample 10% from each group from variable B, no seed
stratified(dat, 1, 2, .1)
# A  B     C   D  E
# 2   6 AA -0.46448941 0.5 TX
# 7  71 AA  1.98128479 2.1 CA
# 5  53 BB  1.00539398 0.7 NY
# 10 97 BB  0.68252675 1.9 NY
# 1   1 CC -0.07870439 0.6 NY
# 4  42 CC -2.00256854 0.3 TX
# 8  76 DD -0.84151459 0.2 NY
# 9  95 DD -0.47276142 0.3 CA
# 11 99 DD  1.05173419 2.1 TX
# 3  10 EE -0.69079473 1.1 TX
# 6  57 EE -0.38210921 1.5 CA
# Sample 10% from each group from variable E, seed of 1
stratified(dat, 1, 5, .1, 1)
# A  B    C   D  E
# 4  33 AA  1.6105099 0.5 CA
# 7  48 AA  0.3128274 0.6 CA
# 9  62 DD  0.4673061 0.0 CA
# 10 86 EE  0.4047880 1.6 CA
# 3  28 AA -1.6815553 0.3 NY
# 5  36 AA  0.3307508 0.3 NY
# 8  53 BB  1.0053940 0.7 NY
# 1  21 DD  0.5229282 1.2 TX
# 2  27 BB  0.8678977 0.7 TX
# 6  44 DD -0.5790353 0.9 TX
# You can also be verbose if it helps you remember what you're doing
stratified(df=dat, id=1, group=5, size=.1, seed=1)
# A  B    C   D  E
# 4  33 AA  1.6105099 0.5 CA
# 7  48 AA  0.3128274 0.6 CA
# 9  62 DD  0.4673061 0.0 CA
# 10 86 EE  0.4047880 1.6 CA
# 3  28 AA -1.6815553 0.3 NY
# 5  36 AA  0.3307508 0.3 NY
# 8  53 BB  1.0053940 0.7 NY
# 1  21 DD  0.5229282 1.2 TX
# 2  27 BB  0.8678977 0.7 TX
# 6  44 DD -0.5790353 0.9 TX

IMPORTANT: This is here mostly to remind me of how I solved my problem. You should read Stratified random sampling in R from a data frame if you really want to use this function.

I know that sampling is quite complex, and I will admit that I know very little about its complexities. Fortunately, software like R lets you draw simple random samples pretty easily, either either with or without replacement. Unfortunately, I could not find any feature to allow me to do simple stratified random sampling, at least not with the features I was looking for. Fortunately again, with a little bit of experimenting, it can be pretty easy to learn how to write functions in R when a direct solution does not present itself.

This post shares my initial “work-in-progress” on writing an R function for stratified sampling.

The problem…

Here’s the minimum that I was hoping for:

• I wanted to be able to draw both a proportional sample (which is more common, for it allows you to make generalizations about the population as a whole) as well as a fixed-size sample (which less common, but it is useful for making comparisons across groups).
• I often use a seed when sampling, so I wanted that to be a part of the function.
• I wanted the output to be the same as if I were to sample from each group individually.
  • I was hoping that my output could be stored as a new object that I could then reuse (either a list or a data frame, preferably the latter).

My initial searches directed me to Yihui Xie’s page on stratified sampling using tapply(). However, this option did not satisfy my needs. As far as I could figure, it only allowed me to take a fixed sample size. Also, I wasn’t totally satisfied with the output.

Consider the following. In Yihui Xie’s example, there is a difference between the results one would get if they sampled from each group separately, but using the same seed.

dat = data.frame(x = 1:15, stratum = gl(3, 5))
dat
##     x stratum
## 1   1       1
## 2   2       1
## 3   3       1
## 4   4       1
## 5   5       1
## 6   6       2
## 7   7       2
## 8   8       2
## 9   9       2
## 10 10       2
## 11 11       3
## 12 12       3
## 13 13       3
## 14 14       3
## 15 15       3
set.seed(1)
tapply(dat$x, dat$stratum, sample, size = 3)
## $`1`
## [1] 2 5 4
##
## $`2`
## [1] 10  6  8
##
## $`3`
## [1] 15 13 12
##
# Compare with what we get when we sample individually:
set.seed(1)
sample(1:5, 3)
## [1] 2 5 4
set.seed(1)
sample(6:10, 3)
## [1]  7 10  9
set.seed(1)
sample(11:15, 3)
## [1] 12 15 14
# But notice this
set.seed(1)
sample(1:5, 3)
## [1] 2 5 4
sample(6:10, 3)
## [1] 10  6  8
sample(11:15, 3)
## [1] 15 13 12

Essentially, the tapply version and the set.seed once version both set the seed for that session or until a new seed is provided. This is fine, but if one does an intermediate step that involves random number generation you won’t get the same results. For example:

set.seed(1)
sample(1:5, 3)
## [1] 2 5 4
sample(6:10, 3)
## [1] 10  6  8
rnorm(10)
##  [1]  1.5953  0.3295 -0.8205  0.4874  0.7383  0.5758 -0.3054  1.5118
##  [9]  0.3898 -0.6212
sample(11:15, 3)
## [1] 11 12 13

Stratified sampling, Mr. DWAB style

The solution I arrived at is to use “unstack()” and a few conditional loops to take the samples.

And, without more rambling, here’s what I came up with.

stratified = function(df, size, seed, dframe=FALSE, ...) {
  # USE: Start with a data frame with your cases in one
  # column and your groups in another column. Decide on
  # if you want to use a seed or not. If not, seed should
  # be "NO" (with quotes). Decide on if you want your
  # output as a data frame or not; by default, dframe is
  # set to "FALSE". To take a sample proportional to the
  # population size in each group, enter "size" as a
  # decimal. Otherwise, enter size as a whole number.
  #
  # Example 1a: To sample 10% of each group from a data
  # frame named "z" and using a seed of "1", use:
  #   > stratified(z, .1, 1)
  # Example 1b: To run the same sample as above but display
  # the result as a data.frame, use:
  #   > stratified(z, .1, 1, T)
  # Example 2: To sample 10% of each group from a data.frame
  # named "z" and using no seed, use:
  #   > stratified(z, .1, "NO")
  # Example 3: To sample 5 from each group from a data.frame
  # named "z" and using a seed of 30, use:
  #   > stratified(z, 5, 30)
  #
  # NOTE: Not recommended for datasets with LOTS of groups
  # or with HUGE differences in group sizes.
  k = unstack(df)
  if (dframe == FALSE) {
    if (seed == "NO" & size < 1) {
      for (i in 1:length(k)) {
        N = k[[i]]
        n = round(length(N)*size)
        pre = structure(list("Group" = names(k[i]),
                             "Population Size" = length(k[[i]]),
                             "Sample Size" = n,
                             Seed = seed,
                             Sample = sample(N, n, ...)),
                        class = "power.htest")
        print(pre)
      }
    } else if (seed == "NO" & size >= 1) {
      for (i in 1:length(k)) {
        N = k[[i]]
        pre = structure(list("Group" = names(k[i]),
                             "Population Size" = length(k[[i]]),
                             "Sample Size" = size,
                             Seed = seed,
                             Sample = sample(N, size, ...)),
                        class = "power.htest")
        print(pre)
      }
    } else if (size < 1) {
      for (i in 1:length(k)) {
        set.seed(seed)
        N = k[[i]]
        n = round(length(N)*size)
        pre = structure(list("Group" = names(k[i]),
                             "Population Size" = length(k[[i]]),
                             "Sample Size" = n,
                             Seed = seed,
                             Sample = sample(N, n, ...)),
                        class = "power.htest")
        print(pre)
      }
    } else if (size >= 1) {
      for (i in 1:length(k)) {
        set.seed(seed)
        N = k[[i]]
        pre = structure(list("Group" = names(k[i]),
                             "Population Size" = length(k[[i]]),
                             "Sample Size" = size,
                             Seed = seed,
                             Sample = sample(N, size, ...)),
                        class = "power.htest")
        print(pre)
      }
    }
  } else {
    if (seed == "NO" & size < 1) {
      for (i in 1:length(k)) {
        N = k[[i]]
        n = round(length(N)*size)
        res = data.frame(names(k[i]), sample(N, n, ...))
        names(res) = c("Group", "Samples")
        print(res)
      }
    } else if (seed == "NO" & size >= 1) {
      for (i in 1:length(k)) {
        N = k[[i]]
        res = data.frame(names(k[i]), sample(N, size, ...))
        names(res) = c("Group", "Samples")
        print(res)
      }
    } else if (size < 1) {
      for (i in 1:length(k)) {
        set.seed(seed)
        N = k[[i]]
        n = round(length(N)*size)
        res = data.frame(names(k[i]), sample(N, n, ...))
        names(res) = c("Group", "Samples")
        print(res)
      }
    } else if (size >= 1) {
      for (i in 1:length(k)) {
        set.seed(seed)
        N = k[[i]]
        res = data.frame(names(k[i]), sample(N, size, ...))
        names(res) = c("Group", "Samples")
        print(res)
      }
    }
  }
}

You can load the function by typing:

source("http://news.mrdwab.com/stratified-beta")

And now, to test it…

Let’s generate some dummy data and see what we can come up with. The function takes the following arguments (in the following order):

• df: The source data frame, with the first column being the IDs and the second column being the groups.
• size: The sample size you want, either as a percentage (for proportional sampling–expressed as a decimal) or as a whole number.
• seed: The seed you want to use. If you don’t want to use a seed, enter “NO“.
• dframe: What format you want the output in, either a list or a data frame. Defaults to a list (dframe=FALSE), which is better at the moment since the data frame option is not working the way I expect it to yet.

Here we go!

# Generate some data
a = 1:100
set.seed(123)
b = sample(c("a", "b", "c", "d"), 100, replace = T)
z = data.frame(a, b)
# Check how big each group is
table(z$b)
##
##  a  b  c  d
## 26 27 20 27
# Make sure the function is loaded before you continue!
# source("http://news.mrdwab.com/stratified-beta")
# Take a 15% sample and use a seed of 1
stratified(z, .15, 1)
##
##
##
##           Group = a
## Population Size = 26
##     Sample Size = 4
##            Seed = 1
##          Sample = 38, 45, 54, 81
##
##
##
##
##           Group = b
## Population Size = 27
##     Sample Size = 4
##            Seed = 1
##          Sample = 39, 43, 60, 79
##
##
##
##
##           Group = c
## Population Size = 20
##     Sample Size = 3
##            Seed = 1
##          Sample = 23, 26, 33
##
##
##
##
##           Group = d
## Population Size = 27
##     Sample Size = 4
##            Seed = 1
##          Sample = 21, 31, 53, 71
##
# Take a sample of 5 from each group and use a seed of 1
stratified(z, 5, 1)
##
##
##
##           Group = a
## Population Size = 26
##     Sample Size = 5
##            Seed = 1
##          Sample = 38, 45, 54, 81, 30
##
##
##
##
##           Group = b
## Population Size = 27
##     Sample Size = 5
##            Seed = 1
##          Sample = 39, 43, 60, 79, 19
##
##
##
##
##           Group = c
## Population Size = 20
##     Sample Size = 5
##            Seed = 1
##          Sample = 23, 26, 33, 78, 14
##
##
##
##
##           Group = d
## Population Size = 27
##     Sample Size = 5
##            Seed = 1
##          Sample = 21, 31, 53, 71, 11
##
# Take a sample of 15 from each group, with replacement, and a seed of 1
stratified(z, 15, 1, replace=T)
##
##
##
##           Group = a
## Population Size = 26
##     Sample Size = 15
##            Seed = 1
##          Sample = 38, 45, 56, 91, 35, 91, 96, 74, 62, 15, 35, 30, 74, 45, 81
##
##
##
##
##           Group = b
## Population Size = 27
##     Sample Size = 15
##            Seed = 1
##          Sample = 39, 44, 63, 93, 29, 93, 95, 66, 64, 3, 29, 19, 70, 44, 77
##
##
##
##
##           Group = c
## Population Size = 20
##     Sample Size = 15
##            Seed = 1
##          Sample = 23, 26, 55, 94, 22, 92, 94, 72, 61, 9, 22, 14, 72, 26, 78
##
##
##
##
##           Group = d
## Population Size = 27
##     Sample Size = 15
##            Seed = 1
##          Sample = 21, 32, 58, 88, 16, 88, 89, 65, 59, 4, 16, 11, 67, 32, 69
##
# Take a sample of 10% from each group, using a seed of 1,
# and display the output as a data frame
stratified(z, .1, 1, dframe=T)
##   Group Samples
## 1     a      38
## 2     a      45
## 3     a      54
##   Group Samples
## 1     b      39
## 2     b      43
## 3     b      60
##   Group Samples
## 1     c      23
## 2     c      26
##   Group Samples
## 1     d      21
## 2     d      31
## 3     d      53

Everything seems OK.

Replicating the results from tapply()

I mentioned earlier that the results are different from what you would get if you were to use the tapply() function. However, it is easy to get the same results using this stratified function–simply move your “seed” outside of the function (enter seed as "NO" [with quotes] and instead, use set.seed() as you normally would).

# See what tapply() gives us
set.seed(1); tapply(z$a, z$b, sample, size = 4)
## $a
## [1] 38 45 54 81
##
## $b
## [1] 29 86 95 63
##
## $c
## [1] 61  9 14 92
##
## $d
## [1] 67 31 68 34
##
# The normal usage for the stratified function
stratified(z, 4, 1, dframe = T)
##   Group Samples
## 1     a      38
## 2     a      45
## 3     a      54
## 4     a      81
##   Group Samples
## 1     b      39
## 2     b      43
## 3     b      60
## 4     b      79
##   Group Samples
## 1     c      23
## 2     c      26
## 3     c      33
## 4     c      78
##   Group Samples
## 1     d      21
## 2     d      31
## 3     d      53
## 4     d      71
# Getting the same results as tapply()
# Set the seed before using the function,
# and set the seed for the function as "NO"
set.seed(1); stratified(z, 4, "NO", dframe = T)
##   Group Samples
## 1     a      38
## 2     a      45
## 3     a      54
## 4     a      81
##   Group Samples
## 1     b      29
## 2     b      86
## 3     b      95
## 4     b      63
##   Group Samples
## 1     c      61
## 2     c       9
## 3     c      14
## 4     c      92
##   Group Samples
## 1     d      67
## 2     d      31
## 3     d      68
## 4     d      34

So, set your seed according to what you want to replicate.

The unfortunate…

There are some advantages to each of the output formats. I’ve set up the list to be quite verbose, which is useful with the proportionate sampling since it shows us how many samples have been taken from each group. The data frame output format, on the other hand, is quite compact.

What I still need to figure out, though, is why R won’t store my output. I suspect that it has something to do with how my loops are set up. I assume that somewhere, I need to add something like an rbind command.

When the time is right, I will be sure to post what I’ve found.