The trouble with ‘ntile()’

Notes from a data witch，Danielle Navarro 的数据科学博客，关于统计学、R 语言与数据可视化。关注
You keep using that function. I do not think it does what you think it does

As longtime readers of this blog will no doubt be aware, I am a tidyverse girl. As deeply as I adore the R programming language, there are some wild things that are baked into the behaviour of base R, and a lot of what tidyverse does so very, very well is smooth over those rough edges. Whether it is data visualisation with ggplot2, data wrangling with dplyr, or navigating the hell that is string manipulation with stringr rather than grapple with the stunningly inconsistent base R regex tools, it has been a blessing. I don’t usually bother to write posts pointing out the obvious, though, so I’ve been more likely to write about a few of the less celebrated packages like fs and cli that sit underneath tidyverse and make everyone’s life a little less painful.

In short, I am a fan.

And so it is with a certain level of trepidation and regret that I find myself writing a post about something that a tidyverse package doesn’t do well. After all, who am I to criticise? I have written some truly terrible code over the years and had some very unwise choices end up in packages that people actually use. But of course that’s the point… even the very best tools have problems, even the best programmers make mistakes, etc. It’s foolish to pretend otherwise.

With that in mind, this is a post about dplyr::ntile(). It is a warning that you should not use it for any purpose where statistical accuracy is important. If you try to guess what it does based on the function name – and are one of the fortunate souls never to have worked with SQL – you will get burned when you discover what it actually does. It is emphatically not a tool for binning data into quantile-based groups, and it will absolutely misbehave when you use it for data analysis. Please be careful.

library(dplyr)
library(ggplot2)

A convenient data set

I’ll start by introducing a fictitious dataset er_data that is deliberately designed to slightly exaggerate the problems that show up in the wild when using ntile(). The data set mimics the kind of thing you might encounter when doing an exposure-response (ER) analysis in pharmacometric work. The setup is this, which I’ll over-explain since most people reading my blog don’t work in drug development. We have data from three studies:

  • Study S01 is a phase one dose escalation study with a “single ascending dose” design. The details don’t matter too much for the toy example but the key thing here is that it’s a sequential design. The first batch of subjects receive a very low dose. If there are no adverse events, the next batch gets a higher dose. And so on. Including a study with this kind of design is critical for safety purposes, so if your outcome variable is a safety endpoint there will probably be something like this in your data set. Usually, all the subjects in this study will be male:1 you do not want to risk the possibility of someone in this study being pregnant and unaware of it.
  • Study S02 is a “drug-drug interaction” study, again part of the phase 1 work. This one is specifically related to oral contraceptives: we do not want to risk the possibility that our new drug alters the effectiveness of birth control pills (or vice versa). So there’s a good chance that you’ll see a study like this too: unsurprisingly, all participants in this study are female, and the dose level is fixed.
  • Study S03 is a phase 2 “dose-finding” study, of the sort typically conducted to pin down the efficacy of the drug. There’s a broader range of participants here (men and women are both included), and there’s variation in the dose level also.

In this situation, the tabulation of subject count by study and gender might look a bit like this:

er_data |> count(study_id, phase, description, sex)
# A tibble: 4 × 5
  study_id phase description            sex       n
                          
1 S01          1 SAD dose-escalation    M        31
2 S02          1 DDI oral contraceptive F        24
3 S03          2 Dose-finding           F        60
4 S03          2 Dose-finding           M        60

Since our fictitious data set loosely mirrors an exposure-response analysis scenario, the data set includes some typical measures of drug exposure (e.g., cmax represents peak drug concentration, auc measures the “area under the curve” measure of total drug exposure over some period of time), and some response measure that we are interested in (e.g., a safety measure, an efficacy measure etc). It’s not super-important for the current post, but just to give you a sense of it, this is what the exposure-response relationship looks like in this data:

er_data |> 
  ggplot(aes(auc, response)) + 
  geom_point(aes(color = study_id)) +
  geom_smooth(formula = y ~ x, method = "lm", color = "#222")

The higher the drug exposure, the stronger the response. That’s usually the thing we’re interested in when conducting an exposure-response analysis but it’s not central to the current post so I’ll move along.

What is central to the current post, however, is that our data set is organised in a systematic, sensible way. Each row in er_data corresponds to a specific subject, and each column corresponds to a particular measurement. Because the data programmer who prepared this data set is not pointlessly cruel, the rows are not ordered randomly. Instead they are arranged in a fashion that makes it easy for the analyst to understand: rows are ordered by study_id, then by dose_mg, and then by subject_id:

er_data
# A tibble: 175 × 11
   row_id study_id phase description         subject_id sex   dose_mg wt_kg   auc   cmax response
                                          
 1      1 S01          1 SAD dose-escalation S01-001    M          10  88.5 1.87  0.265      16.8
 2      2 S01          1 SAD dose-escalation S01-002    M          10  82.3 1.71  0.262      29.1
 3      3 S01          1 SAD dose-escalation S01-003    M          10 117   1.21  0.129      16.9
 4      4 S01          1 SAD dose-escalation S01-004    M          10  85.7 0.871 0.0828     19.0
 5      5 S01          1 SAD dose-escalation S01-005    M          10  83.9 1.54  0.147      16.6
 6      6 S01          1 SAD dose-escalation S01-006    M          10  77.2 2.21  0.136      29.3
 7      7 S01          1 SAD dose-escalation S01-007    M          10  77.2 1.31  0.106      19.8
 8      8 S01          1 SAD dose-escalation S01-008    M          30  87.9 4.93  0.426      26.1
 9      9 S01          1 SAD dose-escalation S01-009    M          30  84.1 4.46  0.494      20.3
10     10 S01          1 SAD dose-escalation S01-010    M          30  77.2 4.93  0.411      19.7
# ℹ 165 more rows

The column names tell you exactly what each variable represents:

  • row_id exists for bookkeeping purposes, and contains the original row number
  • study_id indicates which study the data come from
  • phase indicates whether this is a phase 1 study or a phase 2 study
  • description gives a brief description of the study
  • subject_id provides a unique identifier for each subject
  • sex indicates whether the person is male or female
  • dose_mg specifies the dose they were given (in milligrams)
  • wt_kg specifies their body weight (in kilograms)
  • auc and cmax are the two exposure metrics
  • response is the response variable named in the least imaginative way possible

Again, most of this isn’t germane to the point. The critical thing to draw your attention to is the wt_kg variable, which only records a person’s body weight to the nearest 10th of a kilogram. It is rounded to one decimal point because real-world scales generally report weight at that level of precision.2 Consequently, even though weight is “theoretically” a continuously-varying quantity it is not even remotely so in real data sets: in every real life data analysis we end up with quite a few people recorded as having “identical” weights. This will happen most often near the middle of the distribution:

er_data |> filter(wt_kg == median(wt_kg))
# A tibble: 14 × 11
   row_id study_id phase description            subject_id sex   dose_mg wt_kg   auc  cmax response
                                            
 1      6 S01          1 SAD dose-escalation    S01-006    M          10  77.2  2.21 0.136     29.3
 2      7 S01          1 SAD dose-escalation    S01-007    M          10  77.2  1.31 0.106     19.8
 3     10 S01          1 SAD dose-escalation    S01-010    M          30  77.2  4.93 0.411     19.7
 4     12 S01          1 SAD dose-escalation    S01-012    M          30  77.2  4.14 0.252     11.6
 5     23 S01          1 SAD dose-escalation    S01-023    M         100  77.2 16.1  1.34      42.6
 6     38 S02          1 DDI oral contraceptive S02-038    F         100  77.2 14.2  1.65      36.8
 7     40 S02          1 DDI oral contraceptive S02-040    F         100  77.2 26.1  2.62      41.8
 8     41 S02          1 DDI oral contraceptive S02-041    F         100  77.2 13.3  1.42      36.5
 9     48 S02          1 DDI oral contraceptive S02-048    F         100  77.2 19.7  1.68      38.5
10     71 S03          2 Dose-finding           S03-071    F          50  77.2  8.99 0.731     32.1
11     86 S03          2 Dose-finding           S03-086    M          50  77.2 17.3  1.29      34.2
12    146 S03          2 Dose-finding           S03-146    M         200  77.2 31.6  2.79      74.9
13    155 S03          2 Dose-finding           S03-155    F         200  77.2 39.1  2.57      54.5
14    156 S03          2 Dose-finding           S03-156    M         200  77.2 34.0  4.19      57.6

Yep, we have ties.

But wait… 14 tied values at the median? In a data set with only 175 rows? That’s 8% of the data set. That would come as something of a surprise in real life, but it’s entirely to be expected when the author of the post has placed her thumb on the scales and set up her data set in a manner that exaggerates the issue she’s trying to document. These ties are the exact thing that will create the problem I’m about to write about, so I set up the data set to help make it a little easier to see the problem with ntile().

That being said, although my example data are a little contrived, they are not grotesquely unrealistic. While it’s a bit unlikely that you’d see 8% of a real-life data set with tied values at the median, I have absolutely encountered real-world data where 2% of the sample clusters at the median weight. This “clump of tied values at the median” situation is in fact quite common in the wild.

With this as the tediously long preamble, let’s do some data analysis and make the rather unfortunate mistake of applying the ntile() function as part of it…

The trouble with ntile()

Okay, let’s start by doing something fairly ordinary. It is quite common when analysing pharmacometric data to group subjects into weight-based bins and look at exposures separately by weight bin. Something like this perhaps?

er_data |> 
  mutate(wt_bin = factor(ntile(wt_kg, n = 4))) |> 
  ggplot(aes(wt_bin, auc)) + 
  geom_boxplot() + 
  facet_wrap(~study_id)

Again, the details don’t matter much and I’m not even trying to make this super realistic. The point is that it is very common to do something vaguely similar to this: group a “continuous” variable into bins, and then summarise the data separately by bin.

Now we get to the unpleasant part. Because our “continuous” variable wt_kg is not actually continuous and is measured only up to a low level of precision, our data set contains all those tied values. As we saw earlier, there are 14 people who sit “exactly” at the median weight. What happens to those people? Because their weight is precisely at the median, it is perfectly sensible to assign them to bin 2, but equally sensible to assign them to bin 3. When binning “continuous” variables into discrete categories, you have to make some decision about what to do about the tied values that just happen to sit on a break point.

Okay… so what is the “correct” answer to this decision? Well… it depends on who you ask. If, for example, you were to ask a statistician, you would find yourself trapped in verrrrry long and painful discussion about estimating distribution quantiles from order statistics, and there is a good chance that will cause you to burst into tears. So perhaps you decide not to ask a statistician, and instead you ask what ntile() has to say about that topic…

er_data |> 
  mutate(wt_bin = ntile(wt_kg, n = 4)) |> 
  filter(wt_kg == median(wt_kg)) |> 
  split(~wt_bin)




  
  
  
  
  芦苇 - 优质的中文分类社区
  
  
  
  
  
  
  
  



  2`
# A tibble: 4 × 12
  row_id study_id phase description         subject_id sex   dose_mg wt_kg   auc  cmax response wt_bin
                                          
1      6 S01          1 SAD dose-escalation S01-006    M          10  77.2  2.21 0.136     29.3      2
2      7 S01          1 SAD dose-escalation S01-007    M          10  77.2  1.31 0.106     19.8      2
3     10 S01          1 SAD dose-escalation S01-010    M          30  77.2  4.93 0.411     19.7      2
4     12 S01          1 SAD dose-escalation S01-012    M          30  77.2  4.14 0.252     11.6      2





  
  
  
  
  芦苇 - 优质的中文分类社区
  
  
  
  
  
  
  
  



  3`
# A tibble: 10 × 12
   row_id study_id phase description            subject_id sex   dose_mg wt_kg   auc  cmax response wt_bin
                                              
 1     23 S01          1 SAD dose-escalation    S01-023    M         100  77.2 16.1  1.34      42.6      3
 2     38 S02          1 DDI oral contraceptive S02-038    F         100  77.2 14.2  1.65      36.8      3
 3     40 S02          1 DDI oral contraceptive S02-040    F         100  77.2 26.1  2.62      41.8      3
 4     41 S02          1 DDI oral contraceptive S02-041    F         100  77.2 13.3  1.42      36.5      3
 5     48 S02          1 DDI oral contraceptive S02-048    F         100  77.2 19.7  1.68      38.5      3
 6     71 S03          2 Dose-finding           S03-071    F          50  77.2  8.99 0.731     32.1      3
 7     86 S03          2 Dose-finding           S03-086    M          50  77.2 17.3  1.29      34.2      3
 8    146 S03          2 Dose-finding           S03-146    M         200  77.2 31.6  2.79      74.9      3
 9    155 S03          2 Dose-finding           S03-155    F         200  77.2 39.1  2.57      54.5      3
10    156 S03          2 Dose-finding           S03-156    M         200  77.2 34.0  4.19      57.6      3

Oh.

Oh no.

That doesn’t look at all right, does it?

It’s not very well documented on the package website, but the ntile() function breaks ties based on the row number. Computationally I imagine this is quite a convenient thing to do, but unfortunately what that means in the real world is that your “quantile groups” are not based on any recognised estimator of distribution quantiles. For a realistic data set where observations are arranged in a structured, meaningful way – i.e., every properly formatted data set – ntile() breaks ties in a highly non-random way.

In our scenario, 4 of the 14 median-tied values end up allocated to Q2, and all four of those come from the phase 1 dose-escalation study. They are all male. And even within the dose-escalation study, shall we take a look at which observation didn’t get allocated to Q2? It’s the one from the 100 mg dose condition. Remember, study S01 is the dose-escalation study and it is a sequential design. The data set very reasonably orders the subjects within this study by dose level, so the decision about which of the S01 subjects is labelled as “below median weight” is guaranteed to be the one given the highest dose of the drug. Yikes.

In other words, while you might be tempted to believe that using ntile() has partitioned our data set into quantile-based weight bins, it has not done so. Our wt_bin variable is something that looks like a weight quartile group, but it is something else entirely. Because our data set has tied values at the median weight, and because our data are organised in a sensible fashion, the most accurate way to describe what we have done in a report that we send to the US FDA would look something like this:

We assigned subjects to “weight” quartiles that were mostly based on baseline weight, but also a little bit based on study id, and also a tiny bit based on the dose

I feel quite certain that neither the FDA nor any other regulatory authority on the planet would consider this to be an acceptable way of doing the analysis. Sure, if all you’re doing is preparing a quick-and-dirty dashboard for a C-suite executive who doesn’t really care much about the details, go for it… use ntile() if you like. I’m not going to stop you, and quite likely neither will your boss. But if you’re doing an analysis that will be used to decide if a new medication should be approved for humans? Choose something else.

What should I do instead?

Okay, Danielle, I agree with you that ntile() is not fit for purpose here, but what should I do instead? Well, Other Danielle, first thank you for a great question, and second I think it helps to recognise that there are two related but genuinely distinct concerns with what ntile() does when binning observations.

  • The tie-breaking rule is statistically indefensible. The core of my objection above is that when ntile() encounters tied values at a break point, it applies a tie-breaking rule (break ties based on row number) that has appalling statistical properties in real life.3 If you are going to have a tie-breaking procedure built into your classification rule, it simply cannot be something that is both (a) barely-documented and (b) highly non-random. That’s really, really bad.
  • If you want your groups to be based on quantiles, you should not break ties. The second point is a little more subtle. In most real world analyses our goal is not actually to create “groups of equal size”. The goal is to split the data based on the quantiles calculated from the sample. If you have large samples and genuinely-continuous data, it’s a distinction without much of a difference. But when you have quasi-continuous data like the er_data dataset I’m using in this blog post, these two things diverge. A quantile-based split is not the same thing as an equal-size split.

Either one of these is a good reason not to use ntile() but they differ in the severity of the problems they cause. The “non-random tie-breaking” problem is deeply dangerous, and automatically disqualifies ntile() as a tool for statistical analysis. But suppose ntile() were to instead break-ties randomly? Would that truly be so bad to have an “equal-size” classification rule as long as tie-breaking were handled in a sensible way?

In all honesty, I think… probably not? There’s no fundamental law of statistics that says that “thou shalt send quantile-splits to the FDA and not equal-size splits”. Any time you carve a continuous variable into discrete categories you’re necessarily losing some degree of fidelity, so in a sense it seems silly to worry too much about the difference between a quantile-based split and a (properly constructed) equal-size split.4 It is mostly by convention that we use quantile-splits rather than equal-size splits. Nevertheless, if we intend to follow the convention then we actually do need to follow the convention: you cannot perform an equal-size split and then report to the FDA that you performed a quantile-split. That’s not okay.

Implementing quantile-splits

Okay, so let’s assume we’ve decided to implement quantile-based groups, and as such we start with the principle that we should not break ties at all. Because of this, we need to make a decision at the beginning about whether edge cases “fall upwards” or “fall downwards”.5 It doesn’t matter which direction we choose, but we do have to choose. If two people are recorded as having the same weight we must not allocate them to different weight bins. Yes, that means that our quartile-split will not contain the same number of subjects, and in some cases the group sizes can be highly uneven, but that is perfectly fine.

In practice then, we might write a helper function that looks something like this:

cut_quantile <- function(x, n = 4, right_closed = TRUE) {
  breaks <- quantile(x, probs = (0:n)/n, na.rm = TRUE)
  bin_num <- as.numeric(cut(x, breaks, labels = 1:n, include.lowest = TRUE, right = right_closed))
  bin_fct <- factor(bin_num, levels = 1:n, labels = paste0("Q", 1:n)) 
  return(bin_fct)
}

In the code above we have a clear division of responsibilities: quantile() is responsible for calculating the empirical quantiles of the data set, and cut() is responsible for assigning each observation to one of the bins based on its value. This division of roles is important:

  • quantile() ensures that there is always a clearly-defined rule that defines the bins, based solely on the distributional properties of the data and an estimator that the analyst trusts
  • cut() ensures that the rule is applied in exactly the same way to every observation in the data set, based only on its value and not upon any extraneous properties

When applied to er_data where we have 175 subjects, quantile-binning will absolutely not produce four groups of almost-equal size. The sheer number of tied values means we will almost certainly not end up with a 44/44/44/43 split. The purpose of quantile-binning is not to create equal sized groups in the sample, it is to define groups based on our best guess about what an equal-split would look like in the population. It is fundamentally an inferential statistic.6 As a consequence, here’s what happens when we use the cut_quantile() function to bin the data by weight quartile:

er_data |> 
  mutate(wt_bin = cut_quantile(wt_kg)) |> 
  summarise(.by = wt_bin, n_subj = n(), min_wt = min(wt_kg), max_wt = max(wt_kg))
# A tibble: 4 × 4
  wt_bin n_subj min_wt max_wt
         
1 Q4         44   85.6  117  
2 Q3         33   77.5   85.5
3 Q2         54   67.6   77.2
4 Q1         44   52.3   67.4

The key thing here is to notice the minimum and maximum value within each bin. The maximum value of every weight bin is always strictly lower than the minimum value of the next highest weight bin. It must be so: even though our data set has a lot of subjects recorded at identical weights, we do not break ties, and so any tied values at a break always move in the same direction as one another. In this example they always move down because cut() defaults to “left-open right-closed” intervals: all 14 subjects with weight 77.2 kg are assigned to “Q2” and none of them are assigned to “Q3”. Consequently, the minimum weight in the “Q3” group is slightly higher, at 77.5 kg.

You should always expect this pattern for quantile-based bins: it follows as a logical requirement that there be a deterministic rule that maps the value of an observation to the bin to which it is assigned. Because we are doing quantile-based binning, we are not permitted to break ties and you cannot see the same number appear twice in that table.

By contrast, let’s look at what ntile() does:

er_data |> 
  mutate(wt_bin = factor(paste0("Q", ntile(wt_kg, 4)))) |> 
  summarise(.by = wt_bin, n_subj = n(), min_wt = min(wt_kg), max_wt = max(wt_kg))
# A tibble: 4 × 4
  wt_bin n_subj min_wt max_wt
         
1 Q4         43   85.7  117  
2 Q3         44   77.2   85.6
3 Q2         44   67.6   77.2
4 Q1         44   52.3   67.4

Okay, yes. This time we get the 44/44/44/43 split that we might naively have expected, and even I will concede it is quite satisfying to see that nice little tabulation. But it really is the wrong thing to do here if the goal is to implement quantile-binning. Instead of looking at the counts in n_subj and sighing with pleasure, you should be looking at the min_wt and max_wt values within each “quantile bin” and recoiling in horror. Because when you look at these columns, the first thing you notice is that the heaviest subject in “Q2” has… exactly the same weight as the lightest subject in “Q3”? It’s not a rounding error, it’s a genuine collision. Two subjects with precisely the same measured weight have been allocated to different “weight bins”. That’s nonsensical. I mean… it’s simply wrong. Those nice-looking bin counts that we obtain by using ntile() are masking a genuine statistical error: this is not a quantile split.7

On convention

At this point in the post there’s not a lot left to say except unpack some small details and then end with a very snarky conclusion. Let’s do the small details first, yes?

Because at this point we have been clear that our intention is to define quantile-bins rather than equal-split bins, and we now have a tool for doing so, there is only one substantive decision to make: should tied values at a break point “fall upwards” (left-closed, right-open intervals) or should they “fall downwards” (left-open, right-closed intervals)? As you can probably imagine, it really really does not matter, but it helps a bit to unpick the hidden “conventions” that make absolutely no scientific difference whatsoever but can easily lead you astray if you start mistaking a convention for something meaningful.

Our thoroughly tedious story begins with the conventions for defining a “distribution”. By convention we define the cumulative distribution function (CDF) at for a random variate to be the probability that the observed value for that variable does not exceed . Applied to the empirical distribution of weights within the er_data data set, we can easily plot the empirical CDF with the assistance of the rather handy little base R function factory ecdf(). Here it is:

# use the in-built ecdf() function to build a custom function
# for the empirical CDF of the weights in our data set
ecdf_wt <- ecdf(er_data$wt_kg)

# a little bit of data wrangling to make nice plottable data
plot_df <- bind_rows(
  in_data = distinct(er_data, wt_kg), 
  extends = tibble(wt_kg = c(45, 125)),
  .id = "type"
) |> 
  arrange(wt_kg) |> 
  mutate(
    wt_cdf = ecdf_wt(wt_kg),
    label = wt_kg %in% c(77.1, 77.2, 77.5),
    label_x = case_when(
      wt_kg == 77.1 ~ 77.7,
      wt_kg == 77.2 ~ 76.8,
      wt_kg == 77.5 ~ 78.0,
      .default = NA
    )
  )

# show the empirical CDF for the weights data, with the locations of the 
# classification-relevant probabilities highlighted
base_plt <- plot_df |> 
  ggplot(aes(wt_kg, wt_cdf)) + 
  geom_hline(yintercept = 0:4/4, linetype = "dashed", color = "red") + 
  geom_point(data = plot_df |> filter(type == "in_data"), size = 2) +
  labs(
    x = "Weight (kg)", 
    y = "Empirical CDF",
  )
base_plt + 
  geom_step(direction = "hv") + 
  annotate(
    geom = "label", 
    x = 50, 
    y = c(.125, .375, .625, .875), 
    label = paste0("Q", 1:4)
  )

Just looking at that plot you can easily see where the author has “placed her thumb on the scale” when constructing the toy data set. The sharp jump in the ECDF function that occurs at 77.2 kg – which just so happens to be the sample median – is not natural. That’s the part I rigged in order to exaggerate the problems with ntile(). In real life you can indeed encounter large-ish jumps in the ECDF at that location but… yeah, not that big. The artificiality of the data set shines through in this plot. But that is not the important thing.

Let’s zoom in on the critical part of the ECDF plot, the region where the ECDF crosses the median value of 77.2 kg. This is what our plot looks like. Each dot in this plot depicts a weight value that is actually observed in the data set. Notice that the definition of the CDF specifically refers to . Since it is obviously the case that , and we observe that 95 of the 175 observations are 77.2 kg or lower, the more precise statement about that critical value is that 56% of people have weight 77.2 kg or lower. That is to say, every single person with a weight of 77.2 kg in our data set is at the 56th weight percentile. Because of that, the little black dot corresponding to 77.2 kg sits above the median…

zoom_plt <- base_plt + 
  geom_label(
    data = plot_df |> filter(label == TRUE), 
    mapping = aes(x = label_x, label = wt_kg)
  ) + 
  coord_cartesian(xlim = c(74, 81), ylim = c(.4, .6))

zoom_plt + geom_step(direction = "hv")

If we agree with this line of reasoning, tied observations at the break point should all “fall upwards”: everybody with 77.2 kg weight should be assigned to “Q3”, not “Q2”. That’s the behaviour we should expect if our quantile-binning rule uses left-closed right-open intervals. Unfortunately that is not the default behaviour of the cut() function in R, which uses right-closed intervals by default. Notice the difference between this…

# use right-closed intervals, matching the cut() default
er_data |> 
  mutate(wt_bin = cut_quantile(wt_kg, right_closed = TRUE)) |> 
  summarise(.by = wt_bin, n_subj = n(), min_wt = min(wt_kg), max_wt = max(wt_kg))
# A tibble: 4 × 4
  wt_bin n_subj min_wt max_wt
         
1 Q4         44   85.6  117  
2 Q3         33   77.5   85.5
3 Q2         54   67.6   77.2
4 Q1         44   52.3   67.4

and this…

# use right-open intervals, matching the implied ecdf() behaviour and also
# the default in santoku::chop_quantiles()  
er_data |> 
  mutate(wt_bin = cut_quantile(wt_kg, right_closed = FALSE)) |> 
  summarise(.by = wt_bin, n_subj = n(), min_wt = min(wt_kg), max_wt = max(wt_kg))
# A tibble: 4 × 4
  wt_bin n_subj min_wt max_wt
         
1 Q4         44   85.6  117  
2 Q3         47   77.2   85.5
3 Q2         40   67.6   77.1
4 Q1         44   52.3   67.4

Based on the way that the black dots have been plotted in the figures above, the second version is the quantile binning that we obtain.

Okay… so does that mean there is an objectively correct answer here? We should use left-open right-closed intervals when defining quantile splits for data analyses to be submitted to regulatory agencies, end of story?

LOL.

No.

LMAO, even.

All of this is convention. For example, if by historical accident we had decided that the definition of the CDF should be rather than , literally every aspect of the previous argument would reverse. In that universe the meaning of the CDF would be subtly changed, but it would still be describing the same mathematical object. Neither universe is morally or statistically superior. It genuinely does not matter if you use right-open or right-closed intervals to define your quartile bins.

To highlight just how arbitrary it all is, notice in the plotting code above I set direction = "hv" as the argument to geom_step(), meaning that the way a step function is drawn is to move horizontally first, and vertically second. If I adopt that convention, the step function representing the cumulative distribution of weights crosses the 50th percentile at 77.2 kg. But no law of gods or man says that direction = "hv" expresses any deep truth about the world.

zoom_plt + geom_step(direction = "vh")

Behold! The weight distribution function now crosses the median at 77.1 kg. Such is the power of changing a into a or a hv into a vh. Which is to say… not very much. Pick one convention and stick with it. You’ll be fine.

Epilogue: You can blame SQL

Looking back on this little exercise, and reflecting about my initial horror at the discovery of what dplyr::ntile() actually does, it seems to me that my complaint is not so much with its behaviour as it is with its name. To anyone with statistical training, the name strongly implies that it is a function for constructing quantile-based bins. What it actually does, however, is mimic the equally-stupid behaviour of the NTILE function in SQL. Just like its SQL counterpart, ntile() forces identical values to be assigned into different groups if that would produce a more “even” split of the data.8 In SQL though, this behaviour leads to even stranger outcomes: because databases don’t really have any inherent concept of row order9 the tie-breaking behaviour in NTILE is simply declared to be “nondeterministic”. All you’re really told about its behaviour is this:

If the number of rows in a partition isn’t divisible by integer_expression, this causes groups of two sizes that differ by one member. Larger groups come before smaller groups in the order specified by the OVER clause. For example, if the total number of rows is 53 and the number of groups is five, the first three groups have 11 rows and the two remaining groups have 10 rows each. If on the other hand the total number of rows is divisible by the number of groups, the rows are evenly distributed among the groups. For example, if the total number of rows is 50, and there are five groups, each bucket contains 10 rows.

NTILE prescribes the size of the groups, but not their membership or the rule that maps values to categories. Despite the painfully inappropriate similarity of name, NTILE is not in any meaningful sense a tool for quantile binning.10 Indeed, if the ORDER BY clause in an NTILE expression produces ties at the breaks… well, who even knows what it will do? It’s “nondeterministic”, which most certainly is not the same thing as “random” in this context. Breaking ties at random would be a lot less unhinged, since the statistical properties of random allocation are well understood. But breaking ties using a rule that is unknown but probably systematically related to variables that matter for your analysis? That is horrifically bad. There are perhaps some contexts in which the “quick and dirty carve up” that NTILE and ntile() provide is useful, but serious data analysis is not one of them.

And without question, neither one is suitable for regulatory work in the pharmaceuticals industry.

Footnotes

  1. Blah blah blah, sex and gender, sex and gender. This is not the post in which I will revisit that topic.↩︎
  2. Besides, even if the data set did report body weight at higher precision, it wouldn’t be meaningful: body weight is time-varying and can change by modest amounts in the course of a few hours. Reporting a subject having a body weight of 72.23482 kilograms is absurd even if they were weighed under fasting conditions while completely naked, and the scale really were that precise.↩︎
  3. In the epilogue, when I connect this back to the true culprit – the SQL NTILE function – it will turn out that what SQL does is even worse.↩︎
  4. At one point in the post I did consider the idea that quantile splits are inherently superior because they have a clearly defined rule that perfectly maps every data value onto exactly one category. And I do still like this as a point of difference. But again, are they really that different? An equal-size rule with random tie breaks has a deterministic rule defined for every data value except the internal break points, and for those it has a properly specified probabilistic procedure for allocating observations; so the difference really is that for the equal-size split you have to pay attention to the RNG state. All in all… it feels like this is a reason to prefer quantile-splits over equal-size splits, but it doesn’t feel like a very strong one.↩︎
  5. Or to be fancy and use the correct language, decide in advance whether the intervals that define your bins are “left-closed right-open” intervals (edge cases move up) or “left-open right-closed” intervals (edge cases move down).↩︎
  6. More on this later. The intuition that people often have about continuous data is that “empirical quantiles” should actually exist, and that these quantiles should be purely descriptive properties of the sample. Real data are not so kind. Our “continuous” data are not actually continuous, and the presence of ties means you need to lean on some statistical theory to work out how to act.↩︎
  7. To foreshadow the epilogue, it turns out that there is a historical reason why ntile() is so bad. The dplyr package is heavily inspired by SQL syntax, and it inherits both the good and the bad parts of SQL. SQL is a terribly elegant language for defining data manipulation operations; but it is completely devoid of any statistical foundation. Though SQL is often referred to as a tool for data analytics it is actually very bad at the analysis part of data analysis. It is precisely the opposite of R, which usually has a very solid statistical justification for its behaviour, but is often very shaky on the programming fundamentals. The dplyr package ends up in this strange situation with regards to ntile() because its foundations in SQL are inappropriate for this particular task.↩︎
  8. I mean, I don’t know why the SQL version does this either, but I am entirely unsurprised that it does something statistically absurd. It was written by programmers, after all. You can’t expect much from them.↩︎
  9. Yes, they have indexes, but that’s a different thing.↩︎
  10. With some amusement, I note that programmers complain endlessly about R being a programming language designed by statisticians who aren’t very good at programming. Well, it turns out that the reverse is also true… SQL supplies a data analytics framework designed by programmers who aren’t very good at data analysis. Call it NTILE or ntile(), it doesn’t matter. It’s statistically batshit and does not belong in any serious data analysis.↩︎

Reuse

CC BY 4.0

Citation

BibTeX citation:

@online{navarro2026,
  author = {Navarro, Danielle},
  title = {The Trouble with “Ntile()”},
  date = {2026-09-20},
  url = {https://blog.djnavarro.net/posts/2026-09-20_the-trouble-with-ntile/},
  langid = {en}
}

For attribution, please cite this work as: Navarro, Danielle. 2026. “The Trouble with ‘Ntile()’.” September 20. .

添加评论
点赞收藏
点踩分享查看原文
评论
?
参与讨论