// HACKER NEWS — CYBERSECURITY
The trouble with 'ntile()'
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.
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.
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:
In this situation, the tabulation of subject count by study and gender might look a bit like this:
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:
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:
The column names tell you exactly what each variable represents: