プログラミング
'Ntile()' の問題点
The Trouble with 'Ntile()' (blog.djnavarro.net)
要約
R言語のdplyrパッケージに含まれるntile()関数は、データを指定された数のグループに分割する際に、統計的な精度が重要な場合には使用すべきではないと警告しています。この関数は、SQLのNTILE()関数とは異なり、データの正確な分位点に基づいたビン分割を行わず、予期しない動作をすることがあります。記事では、薬剤開発の曝露-反応分析を模倣したデータセットを用いて、ntile()がどのように誤解を招く結果を生むかを具体的に解説しています。
全文翻訳
あなたは、その関数を使い続けている。
あなたが思っているような働きはしていないと思う。
このブログの長年の読者なら、私がタイディバース(tidyverse)派であることは間違いなくご存知でしょう。Rプログラミング言語を深く愛していますが、ベースRの挙動にはいくつかのワイルドな側面があり、タイディバースの多くの優れた点はそれらの粗いエッジを滑らかにすることです。ggplot2でのデータ可視化、dplyrでのデータ操作、あるいはstringrでの文字列操作の地獄をナビゲートすること(ベースRの驚くほど一貫性のない正規表現ツールと格闘するのではなく)であれ、それは祝福でした。
しかし、私は通常、明白なことを指摘する投稿はしないので、タイディバースパッケージがあまりうまく機能しないものについて書く傾向があります。結局のところ、私が批判する資格があるでしょうか?私は長年、本当にひどいコードを書いてきましたし、人々が実際に使用するパッケージに、非常に賢明でない選択肢が入り込んでしまったこともあります。しかし、もちろんそれがポイントです…たとえ最高のツールでさえ問題があり、最高のプログラマーでさえ間違いを犯す、など。それを偽ることは愚かです。
それを念頭に置いて、これはdplyr::ntile()に関する投稿です。統計的な精度が重要な目的には使用しないようにという警告です。もしあなたが関数名からその働きを推測しようとし(そしてSQLに一度も携わったことのない幸運な魂の一人であれば)、それが実際に何をするのかを知ったときに痛い目を見るでしょう。それは、データを分位点ベースのグループにビン分割するためのツールでは断じてなく、データ分析に使用すると絶対に誤った動作をします。どうか注意してください。
library(dplyr)
library(ggplot2)
便利なデータセット
まず、n tile()を使用する際に発生する問題をわずかに誇張するように意図的に設計された架空のデータセットer_dataを紹介します。このデータセットは、薬物動態学(pharmacometric)の曝露-反応(ER)分析で遭遇する可能性のあるものを模倣しています。セットアップは次のとおりです。ブログの読者のほとんどは薬物開発に携わっていないため、詳しく説明します。
私たちは3つの研究からのデータを持っています:
研究S01は、フェーズ1の用量漸増試験で、「単回漸増投与」デザインです。詳細はトイ例ではあまり重要ではありませんが、重要なのはそれが逐次デザインであるということです。最初の被験者バッチは非常に低い用量を受け取ります。有害事象がなければ、次のバッチはより高い用量を受け取ります。そして、そのように続きます。この種のデザインの研究を含めることは、安全性の目的で非常に重要です。したがって、もしあなたの結果変数(outcome variable)が安全性エンドポイントであれば、データセットにこのようなものが含まれている可能性が高いです。通常、この研究のすべての被験者は男性です1。妊娠している可能性があり、それに気づいていない人がいるリスクを冒したくありません。
研究S02は、「薬物間相互作用」研究で、これもフェーズ1の作業の一部です。これは特に経口避妊薬に関連しています。新しい薬が避妊薬の効果を変える(またはその逆)リスクを冒したくありません。したがって、このような研究を目にする可能性は十分にあります。驚くことではありませんが、この研究のすべての参加者は女性であり、用量レベルは固定されています。
研究S03は、フェーズ2の「用量決定」研究で、通常、薬の有効性を特定するために実施されるものです。ここでは、より広範な参加者(男性と女性の両方が含まれる)がおり、用量レベルにもばらつきがあります。
この状況では、研究と性別ごとの被験者数の集計は、次のようになります。
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
私たちの架空のデータセットは、曝露-反応分析シナリオを緩やかに模倣しているため、データセットには薬物曝露の典型的な測定値(例:cmaxはピーク薬物濃度を表し、aucは一定期間の総薬物曝露の「曲線下面積」測定値)と、関心のある応答測定値(例:安全性測定値、有効性測定値など)が含まれています。現在の投稿にはあまり重要ではありませんが、感覚をつかむために、これはデータセットの曝露-反応関係の見た目です。
er_data |> ggplot(aes(auc, response)) + geom_point(aes(color = study_id)) + geom_smooth(formula = y ~ x, method = "lm", color = "#222")
薬物曝露が高いほど、応答は強くなります。これは通常、曝露-反応分析を行う際に私たちが関心を持つことですが、現在の投稿の中心ではないため、先に進みます。しかし、現在の投稿の中心は、私たちのデータセットが体系的で合理的な方法で編成されていることです。er_dataの各行は特定の被験者に対応し、各列は特定の測定値に対応します。このデータセットを準備したデータプログラマーは無意味に意地悪ではないため、行はランダムに並べられていません。代わりに、アナリストが理解しやすいように並べられています:行はstudy_id、次にdose_mg、そして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
列名は、各変数が何を表すかを正確に示しています。
row_idは記録保持のために存在し、元の行番号を含んでいます。
study_idはデータがどの研究から来たかを示します。
phaseは、フェーズ1研究かフェーズ2研究かを示します。
descriptionは、研究の簡単な説明を提供します。
subject_idは、各被験者のユニークな識別子を提供します。
sexは、その人が男性か女性かを示します。
dose_mgは、彼らに投与された用量(ミリグラム単位)を指定します。
wt_kgは、彼らの体重(キログラム単位)を指定します。
aucとcmaxは2つの曝露メトリックです。
responseは、最も想像力のない方法で名付けられた応答変数です。
繰り返しになりますが、これらのほとんどはポイントとは関係ありません。現在の投稿で中心となるのは、wt_kg変数です。これは、人の体重を10分の1キログラム単位でしか記録していません。実際のスケールは通常その精度で重量を報告するため、小数点以下1桁に丸められています2。その結果、体重は「理論上」連続的に変化する量であるにもかかわらず、実際のデータセットではまったくそうではなくなっています。すべての実生活のデータ分析で、かなりの数の人々が「同一の」体重を記録していることになります。これは分布の中央付近で最も頻繁に発生します。
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