RとSASのパーセント点の 計算方法の違い

R と SAS のパーセント点の 計算方法の違い 国立循環器病研究センター データサイエンス部 川端 孝典 1

とある日．． • 臨床試験の患者背景（R ユーザー）をダブルチェックしよう（SAS ユーザー） • 検査値 X（連続変数）の要約統計量を計算しよう • 25％点や75％点の値が微妙に一致しない．何故？ ▸ データは全く同じであることは確認済み 4

とある日．． • 臨床試験の患者背景（R ユーザー）をダブルチェックしよう（SAS ユーザー） • 検査値 X（連続変数）の要約統計量を計算しよう • 25％点や75％点の値が微妙に一致しない．何故？ ▸ データは全く同じであることは確認済み • R と SAS でデフォルトのパーセント点の計算方法が違っている！ ▸ しかも，思ったより計算方法のバリエーションがある！ ▸ まとめました！ ▸ R : stats::quantile ▸ SAS : UNIVARIATE procedure ＊今年のSASユーザー会の発表内容と一部被ってしまいました 5

R でパーセント点を計算するオプション ＜USAGE＞ quantile(x, …) # S3 method for default quantile(x, probs = seq(0, 1, 0.25), na.rm = FALSE, names = TRUE, type = 7, …) quantileでは９種類のパーセント点の定義 がある base::summary の内部で quantile 関数を 使っており，quantile.type オプションで 定義を変更できる． https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/quantile graphics::boxplot でも quantile 関数を採 用しているが，オプションの変更はでき なさそう 6

• https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/quantile

SAS でパーセント点を計算するオプション ＜Syntax example＞ proc univariate data=dataX PCLTDEF = 5; var valueX ; output out=summary_pctl pctlpts = 25 50 75 pctlpre = P_; run; SASでは５種類のパーセント点の定義がある PROC SUMMARY, PROC MEANS でも同じ パーセント点の定義を採用している （QNTLDEF=オプション） 引用：Base SAS® 9.3プロシジャガイド：統計プロシジャ https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#pro cstat_univariate_sect028.htm https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#pro cstat_univariate_sect008.htm 7

• https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#pro

R と SAS のパーセント点のオプションの対応 base::quantile PROC UNIVARIATE type = 2 PCTLDEF = 5 type = 4 PCTLDEF = 1 type = 1 type = 3 type = 5 type = 6 type = 7 type = 8 type = 9 PCTLDEF = 3 PCTLDEF = 2 対応なし PCTLDEF = 4 対応なし 対応なし 対応なし • R の quantile 関数で type = 2 とすると SAS のデフォルトに一致するが，R のデ フォルトは SAS で対応していない． ▸ バリデーションが目的なら，基本的 には R が SAS に合わせるのが無難 • Hyndman and Fan (1996) は統計的な性質 の良さの観点から type = 8 を推奨 • SAS のデフォルトについて紹介 8

SAS デフォルト（PCTLDEF=5）のパーセント点の計算アルゴリズム 1. サンプルサイズ 𝑛𝑛 とパーセント点 𝑝𝑝 を決めて何番目のデータを使うか計算する ▸ 𝑗𝑗 番目の順序統計量を 𝑥𝑥(𝑗𝑗) とする ▸ 例：𝑛𝑛 = 1018, 𝑝𝑝 = 0.25 → 𝑛𝑛𝑛𝑛 = 254.5 ▸ 𝟐𝟐𝟐𝟐𝟐𝟐. 𝟓𝟓 番目 PCTLDEF = 5 の定義式 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 𝑗𝑗 は 𝑛𝑛𝑛𝑛 の整数部分，𝑔𝑔 は 𝑛𝑛𝑛𝑛 の小数部分 9

SAS デフォルト（PCTLDEF=5）のパーセント点の計算アルゴリズム 1. サンプルサイズ 𝑛𝑛 とパーセント点 𝑝𝑝 を決めて何番目のデータを使うか計算する ▸ 𝑗𝑗 番目の順序統計量を 𝑥𝑥(𝑗𝑗) とする ▸ 例：𝑛𝑛 = 1018, 𝑝𝑝 = 0.25 → 𝑛𝑛𝑛𝑛 = 254.5 ▸ 𝟐𝟐𝟐𝟐𝟐𝟐. 𝟓𝟓 番目 2. 𝑛𝑛𝑛𝑛 を整数部分と小数部分にわける ▸ 整数部分：𝑗𝑗 = 254，小数部分：𝑔𝑔 = 0.5 PCTLDEF = 5 の定義式 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 𝑗𝑗 は 𝑛𝑛𝑛𝑛 の整数部分，𝑔𝑔 は 𝑛𝑛𝑛𝑛 の小数部分 10

SAS デフォルト（PCTLDEF=5）のパーセント点の計算アルゴリズム 1. サンプルサイズ 𝑛𝑛 とパーセント点 𝑝𝑝 を決めて何番目のデータを使うか計算する ▸ 𝑗𝑗 番目の順序統計量を 𝑥𝑥(𝑗𝑗) とする ▸ 例：𝑛𝑛 = 1018, 𝑝𝑝 = 0.25 → 𝑛𝑛𝑛𝑛 = 254.5 ▸ 𝟐𝟐𝟐𝟐𝟐𝟐. 𝟓𝟓 番目 2. PCTLDEF = 5 の定義式 𝑛𝑛𝑛𝑛 を整数部分と小数部分にわける ▸ 整数部分：𝑗𝑗 = 254，小数部分：𝑔𝑔 = 0.5 3. 定義式に沿ってパーセント点に対応する点を計算する ▸ 小数部分は 𝑔𝑔 = 0.5 > 0 ▸ 𝑥𝑥255 を採用！ 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 𝑗𝑗 は 𝑛𝑛𝑛𝑛 の整数部分，𝑔𝑔 は 𝑛𝑛𝑛𝑛 の小数部分 • 日本語訳：𝑛𝑛𝑛𝑛 がぴったり整数になれば間（平均）を取る．ぴったり整数でなければ一つ上を取る． 11

ところで．． PCTLDEF = 5 の定義式 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 𝑗𝑗 は 𝑛𝑛𝑛𝑛 の整数部分，𝑔𝑔 は 𝑛𝑛𝑛𝑛 の小数部分 12

ところで．． PCTLDEF = 5 の定義式 第一四分位，中央値，第 三四分位の計算方法って こんな定義で勉強されま した？？ 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 𝑗𝑗 は 𝑛𝑛𝑛𝑛 の整数部分，𝑔𝑔 は 𝑛𝑛𝑛𝑛 の小数部分 13

ところで．． PCTLDEF = 5 の定義式 第一四分位，中央値，第 三四分位の計算方法って こんな定義で勉強されま した？？ 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 𝑗𝑗 は 𝑛𝑛𝑛𝑛 の整数部分，𝑔𝑔 は 𝑛𝑛𝑛𝑛 の小数部分 勉強してきたもの，思っているものと異なる パーセント点の定義がデフォルトで使われている！ 14

SAS でパーセント点を計算するオプション（再掲） ＜Syntax example＞ proc univariate data=dataX PCLTDEF = 5; var valueX ; output out=summary_pctl pctlpts = 25 50 75 pctlpre = P_; run; SASでは５種類のパーセント点の定義がある ⇧平均化された経験分布関数 empirical distribution function with averaging. PROC SUMMARY, PROC MEANS でも同じ パーセント点の定義を採用している （QNTLDEF=オプション） 引用：Base SAS® 9.3プロシジャガイド：統計プロシジャ https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#pro cstat_univariate_sect028.htm https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#pro cstat_univariate_sect008.htm 15

• https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#pro

平均化された経験分布関数 𝑥𝑥(𝑗𝑗) , if 𝑔𝑔 = 0 𝑦𝑦 = � 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 > 0 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 Hyndman and Fan (1996)  𝑗𝑗 番目と 𝑗𝑗 + 1 番目の間を取る（平均化）することによって（分位点としての整 合性を保ちながら）経験分布関数の逆関数の連続性を考慮（補正）している 16

我々が習ったパーセント点は？（代表例） 1. 分位点がデータとデータの間に落ちると きは，中央値と同様にその２つのデータ の平均値を分位点とする ▸ 統計学入門（東京大学出版会） ▸ R でも SAS でも対応なし ▸ SAS で重み付きパーセント点を使え ばいけそう https://www.utp.or.jp/book/b300857.html https://product.rakuten.co.jp/product//f388f0adf4489916eb626e39c26e8a62/ 2. 線形補間 ▸ 𝑋𝑋𝑝𝑝 = 1 − 𝑔𝑔 𝑋𝑋(𝑗𝑗) + 𝑔𝑔𝑋𝑋(𝑗𝑗+1) （𝑗𝑗 は 𝑛𝑛𝑛𝑛 の整数部分，𝑔𝑔 は 𝑛𝑛𝑛𝑛 の小数部分） ▸ 医学への統計学【第3版】（朝倉書 店） ▸ PCTLDEF = 1, type = 4 17

• https://www.utp.or.jp/book/b300857.html
• https://product.rakuten.co.jp/product/-

まとめと結論 • R の quantile 関数で type = 2 とすると SAS のデフォルトに一致するが，R のデフォルトは SAS で対 応していない． ▸ バリデーションが目的なら，基本的には R が SAS に合わせるのが無難 • パーセント点のデフォルトの定義は，統計学の教科書の定義と異なっていることがほとんどなので， 使用前にはマニュアルを読んで確認する． • 個人的には，教科書で勉強した（イメージしている）パーセント点と同じで，オプションを指定する だけで再現できる線形補間が使いやすいだろう（PCTLDEF = 1, type = 4） 18

参考 1. Rdocumentation. Stats (version3.6.2). Quantile: Sample Quantiles. https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/quantile (2024/10/18) 2. Base SAS(R) 9.3プロシジャガイド：統計プロシジャ. https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#procstat_uni variate_sect028.htm (2024/10/18) 3. Hyndman, R. J., & Fan, Y. (1996). Sample Quantiles in Statistical Packages. The American Statistician, 50(4), 361–365. https://doi.org/10.2307/2684934. 4. 東京大学教養学部統計学教室 (1991). 統計学入門. 東京大学出版会. P33. 5. 古川俊之 (2013). 医学への統計学【第３版】. 朝倉書店. P18. 20

• https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/quantile
• https://www.sas.com/offices/asiapacific/japan/service/help/webdoc/procstat/viewer.htm#procstat_univariate_sect028.htm