De problemen met ‘ntile()’ – Notities van een data-heks

Je blijft die functie gebruiken, maar ik denk niet dat hij doet wat jij denkt dat hij doet.

Als vaste lezers van deze blog ongetwijfeld weten, ben ik een tidyverse-meisje. Hoeveel ik de programmeertaal R ook bewonder, er zitten sommige vreemde eigenschappen in base R. Veel van wat tidyverse zo goed doet, is het gladstrijken van die ruwe randjes. Of het nu gaat om datavisualisatie met ggplot2, data-wrangling met dplyr, of het navigeren door de hel van stringmanipulatie met stringr (in plaats van te worstelen met de verbazingwekkend inconsistente regex-tools van base R), het is een zegen geweest. Normaal gesproken neem ik niet de moeite om posts te schrijven die het voor de hand liggende benoemen; ik schrijf liever over minder gevierde pakketten zoals fs en cli die onder tidyverse liggen en ieders leven iets minder pijnlijk maken.

Kortom, ik ben een fan.

Daarom schrijf ik met een zeker gevoel van angst en spijt deze post over iets wat een tidyverse-pakket niet goed doet. Immers, wie ben ik om kritiek te leveren? Ik heb door de jaren heen zelf echt verschrikkelijke code geschreven en heb onverstandige keuzes gemaakt die terecht zijn gekomen in pakketten die mensen daadwerkelijk gebruiken. Maar dat is precies het punt: zelfs de beste tools hebben problemen, zelfs de beste programmeurs maken fouten. Het is dwaas om te doen alsof dat niet zo is.

Met dat in gedachten gaat deze post over dplyr::ntile(). Het is een waarschuwing dat je dit niet moet gebruiken voor doeleinden waarbij statistische nauwkeurigheid belangrijk is. Als je probeert te raden wat de functie doet op basis van de naam – en je bent een van de gelukkigen die nooit met SQL hebben gewerkt – kom je bedrogen uit wanneer je ontdekt wat het daadwerkelijk doet. Het is absoluut geen instrument om data in groepen te verdelen op basis van kwantielen, en het zal onvoorspelbaar gedrag vertonen wanneer je het gebruikt voor data-analyse. Wees alsjeblieft voorzichtig.

library(dplyr)
library(ggplot2)

Een handige dataset

Ik begin met het introduceren van een fictieve dataset, er_data, die bewust is ontworpen om de problemen die optreden bij het gebruik van ntile() enigszins uit te vergroten. De dataset bootst het soort data na dat je tegenkomt bij een exposure-response (ER) analyse in farmcometrisch werk.

De opzet is als volgt (ik leg dit uitgebreid uit, aangezien de meeste lezers van mijn blog niet in medicijnontwikkeling werken): we hebben data uit drie studies:

  • Studie S01 is een fase 1 dosis-escalatiestudie met een "single ascending dose"-ontwerp. Dit is een sequentieel ontwerp: de eerste groep proefpersonen krijgt een zeer lage dosis. Als er geen bijwerkingen optreden, krijgt de volgende groep een hogere dosis, enzovoort. Dit is cruciaal voor de veiligheid. Meestal zijn alle proefpersonen in deze studie man; men wil het risico vermijden dat iemand in deze studie zwanger is zonder het te weten.
  • Studie S02 is een "drug-drug interaction"-studie, eveneens onderdeel van het fase 1-werk. Deze is specifiek gericht op orale anticonceptiva: men wil voorkomen dat het nieuwe medicijn de effectiviteit van de pil beïnvloedt (of andersom). In deze studie zijn alle deelnemers vrouw en is het dosisniveau vast.
  • Studie S03 is een fase 2 "dose-finding"-studie, bedoeld om de effectiviteit van het medicijn vast te stellen. Hier is een bredere groep deelnemers (zowel mannen als vrouwen) en is er variatie in het dosisniveau.

In deze situatie zou de tabel van het aantal proefpersonen per studie en geslacht er ongeveer zo uitzien:

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

Omdat onze fictieve dataset een scenario van exposure-response weerspiegelt, bevat deze typische maten voor blootstelling aan het medicijn (bijv. cmax voor de piekconcentratie en auc voor de totale blootstelling over een bepaalde periode) en een responsmaat (bijv. veiligheid of effectiviteit).

Het belangrijkste voor deze post is dat onze dataset op een systematische en logische manier is georganiseerd. Elke rij in erdata correspondeert met een specifieke proefpersoon. De rijen zijn niet willekeurig geordend, maar gerangschikt per studyid, vervolgens per dosemg en tot slot per subjectid:

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
...

De kolomnamen betekenen het volgende:

  • row_id: origineel rijnummer voor de administratie.
  • study_id: uit welke studie de data komen.
  • phase: fase 1 of fase 2 studie.
  • description: korte beschrijving van de studie.
  • subject_id: unieke identifier voor elke proefpersoon.
  • sex: geslacht (man of vrouw).
  • dose_mg: gegeven dosis in milligram.
  • wt_kg: lichaamsgewicht in kilogram.
  • auc en cmax: blootstellingsmetrieken.
  • response: de responsvariabele.

Let vooral op de variabele wt_kg. Deze registreert het gewicht slechts tot op één decimaal nauwkeurig, omdat echte weegschalen meestal op dat niveau rapporteren. Hoewel gewicht theoretisch een continu variërende grootheid is, is dat in echte datasets niet het geval: we komen vaak proefpersonen tegen met "identieke" gewichten. Dit gebeurt het vaakst in het midden van de distributie:

er_data |> filter(wt_kg == median(wt_kg))
# A tibble: 14 × 11
...

Ja, we hebben ties (gelijke waarden). In een dataset van slechts 175 rijen zijn er 14 waarden gelijk aan de mediaan. Dat is 8% van de dataset. In het echt zou dat verrassend zijn, maar ik heb dit bewust zo ingericht om het probleem met ntile() makkelijker zichtbaar te maken.

Het probleem met ntile()

Laten we beginnen met iets heel gewoons. Bij het analyseren van farmcometrische data is het gebruikelijk om proefpersonen in groepen (bins) in te delen op basis van gewicht, om vervolgens de blootstelling per gewichtsgroep te bekijken. Bijvoorbeeld:

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

Het doel is om een "continu" variabel in bins te verdelen en de data per bin samen te vatten. Nu komt het onprettige deel. Omdat wt_kg niet echt continu is, hebben we die gelijke waarden. Er zijn 14 mensen die "exact" op het mediane gewicht zitten. Wat gebeurt er met hen? Omdat hun gewicht precies op de mediaan ligt, is het logisch om hen in bin 2 te plaatsen, maar even logisch om hen in bin 3 te plaatsen.

Wat is het "correcte" antwoord? Als je het aan een statisticus vraagt, beland je in een pijnlijke discussie over het schatten van distributiekwantielen uit orde-statistieken. Laten we in plaats daarvan kijken wat ntile() doet:

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

Het resultaat is dat sommige mensen met het mediane gewicht in bin 2 terechtkomen en anderen in bin 3. Dat ziet er niet goed uit.

De ntile()-functie lost ties op basis van het rijnummer. Computationeel is dit handig, maar in de echte wereld betekent dit dat je "kwantielgroepen" niet gebaseerd zijn op een erkende schatter van distributiekwantielen. In een goed geformatteerde dataset, waar observaties gestructureerd zijn, lost ntile() ties op een zeer niet-willekeurige manier op.

In ons scenario komen 4 van de 14 mediane waarden in Q2 terecht. Deze vier komen allemaal uit de fase 1 dosis-escalatiestudie en zijn allemaal man. Binnen die studie is de beslissing over wie "boven het mediane gewicht" is, gegarandeerd degene die de hoogste dosis van het medicijn kreeg, omdat de data gesorteerd zijn op dosis.

Met andere woorden: hoewel het lijkt alsof ntile() de data in kwantielgroepen heeft verdeeld, is dat niet gebeurd. De meest accurate manier om dit te beschrijven in een rapport voor de FDA zou zijn: "We hebben proefpersonen toegewezen aan 'gewichts'-kwartielen die grotendeels gebaseerd waren op het basisgewicht, maar ook een beetje op het studie-id en een heel klein beetje op de dosis."

Geen enkele regelgevende instantie zou dit accepteren. Voor een snelle dashboard voor een manager is het prima, maar voor een analyse die bepaalt of een medicijn goedgekeurd wordt voor mensen, moet je iets anders kiezen.

Wat kan ik in plaats daarvan doen?

Er zijn twee verschillende problemen met wat ntile() doet:

  1. De regel voor het oplossen van ties is statistisch onverdedigbaar. Het feit dat ties worden opgelost op basis van het rijnummer is in de praktijk verschrikkelijk. Een classificatieregel mag niet gebaseerd zijn op een nauwelijks gedocumenteerde, niet-willekeurige eigenschap.
  2. Als je groepen wilt baseren op kwantielen, moet je ties niet oplossen. In de meeste analyses is het doel niet om "groepen van gelijke grootte" te maken, maar om de data te splitsen op basis van de kwantielen van de steekproef. Bij quasi-continue data zijn een kwantiel-split en een gelijke-grootte-split niet hetzelfde.

Hoewel het conventie is om kwantiel-splitsingen te gebruiken, is het essentieel dat we die conventie ook echt volgen. Je kunt niet een gelijke-grootte-split uitvoeren en vervolgens rapporteren dat je een kwantiel-split hebt gedaan.

Kwantiel-splitsingen implementeren

Als we besluiten kwantiel-gebaseerde groepen te implementeren, hanteren we het principe dat we ties niet mogen oplossen. Dit betekent dat we vooraf moeten beslissen of grensgevallen "omhoog" of "omlaag" vallen. Als twee mensen hetzelfde gewicht hebben, mogen ze niet in verschillende bins worden geplaatst. Dit betekent dat de groepen niet precies even groot zullen zijn, maar dat is volkomen acceptabel.

In de praktijk kunnen we een helperfunctie schrijven:

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)
}

Hier is er een duidelijke verdeling van taken:

  • quantile() berekent de empirische kwantielen op basis van de distributie van de data.
  • cut() wijst elke observatie toe aan een bin op basis van de waarde, en niet op basis van externe eigenschappen (zoals rijnummer).

Wanneer we dit toepassen op er_data, krijgen we geen groepen van exact gelijke grootte. Dat is ook niet het doel. Het doel van kwantiel-binning is om groepen te definiëren op basis van onze beste schatting van hoe een gelijke split in de populatie eruit zou zien. Het is fundamenteel een inferentiële statistiek.

Kijk naar het verschil in resultaten:

Met cut_quantile(): De maximale waarde van elke bin is altijd strikt lager dan de minimale waarde van de volgende bin. Er is geen overlap. Alle 14 proefpersonen met een gewicht van 77,2 kg worden in dezelfde groep (bijv. Q2) geplaatst.

Met ntile(): Hier krijgen we wel de mooie 44/44/44/43 split, maar als je naar de minimum- en maximumwaarden kijkt, zie je iets verschrikkelijks: de zwaarste persoon in Q2 heeft exact hetzelfde gewicht als de lichtste persoon in Q3. Twee mensen met hetzelfde gewicht zijn in verschillende bins geplaatst. Dat is statistisch gezien onjuist.

Over conventies

Er blijft nu nog één beslissing over: moeten gelijke waarden op een breekpunt "omhoog" vallen (links-gesloten, rechts-open intervallen) of "omlaag" (links-open, rechts-gesloten intervallen)?

Dit is puur een kwestie van conventie. De cumulatieve distributiefunctie (CDF) wordt conventioneel gedefinieerd als $\text{Pr}(X \leq x)$. Als we dit strikt volgen, zouden waarden op het breekpunt "omhoog" moeten vallen. In R gebruikt de cut() functie standaard rechts-gesloten intervallen (omlaag).

Het is belangrijk om te beseffen dat dit geen wetenschappelijk verschil maakt. Als we de definitie van de CDF zouden veranderen naar $\text{Pr}(X < x)$, zou alles omdraaien. Geen van beide universums is moreel of statistisch superieur. Kies één conventie en houd je eraan.

Epiloog: Je kunt SQL de schuld geven

Mijn klacht gaat eigenlijk niet over het gedrag van dplyr::ntile(), maar over de naam. Voor iedereen met een statistische opleiding suggereert de naam dat het een functie is voor kwantiel-bins. In werkelijkheid bootst het het gedrag van de NTILE functie in SQL na.

Net als in SQL dwingt ntile() identieke waarden in verschillende groepen als dat nodig is om de groepen "gelijk" te maken. In SQL is dit gedrag zelfs "niet-deterministisch", omdat databases geen inherent concept van rijvolgorde hebben.

NTILE schrijft de grootte van de groepen voor, maar niet de samenstelling of de regel die waarden aan categorieën koppelt. Ondanks de naam is NTILE geen instrument voor kwantiel-binning.

Er zijn contexten waarin deze "snelle en vuile verdeling" nuttig is, maar serieuze data-analyse is daar niet één van. En zonder twijfel is het ongeschikt voor regelgevend werk in de farmaceutische industrie.

***

Bijlagen en voetnoten

Voetnoten:

  1. Over het onderscheid tussen sekse en gender: dat is niet het onderwerp van deze post.
  2. Zelfs bij een hogere precisie zou lichaamsgewicht niet betekenisvol zijn; het varieert gedurende de dag. Een gewicht van 72,23482 kg rapporteren is absurd.
  3. Zie de epiloog: wat SQL doet is nog erger.
  4. Het verschil tussen een kwantiel-split en een (correct geconstrueerde) gelijke-grootte-split is klein, maar belangrijk voor de rapportage.
  5. Beslis of intervallen "links-gesloten rechts-open" of "links-open rechts-gesloten" zijn.
  6. De keuze van de estimator in quantile() kan invloed hebben op de exacte waarde van het kwantiel, maar niet op het fundamentele probleem van de tie-handling.
  7. Een NTILE-stijl gelijke-grootte binning is per definitie geen kwantiel-split.
  8. Afhankelijk van de estimator kan de mediaan variëren (bijv. tussen 77,2 en 77,5), maar alle gangbare estimators in quantile() zijn het erover eens dat de mediaan hier 77,2 kg is.
  9. dplyr is zwaar geïnspireerd door SQL-syntaxis en erft zowel de goede als de slechte eigenschappen. SQL is elegant voor datamanipulatie, maar mist een statistische basis.
  10. Databases hebben geen inherente rijvolgorde zonder een ORDER BY clausule.
  11. Het is ironisch: programmeurs klagen dat R is ontworpen door statistici die niet goed kunnen programmeren. Het omgekeerde is ook waar: SQL is ontworpen door programmeurs die niet goed zijn in data-analyse.

Referenties: