The trouble with ‘ntile()’ – Notes from a data witch

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 grappling 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 <chr> <dbl> <chr> <chr> <int> 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 <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> 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 <dbl> <chr> <dbl> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> 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))) |>