769 Views
October 21, 24
スライド概要
[第10回大阪sas勉強会]
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では9種類のパーセント点の定義 がある base::summary の内部で quantile 関数を 使っており,quantile.type オプションで 定義を変更できる. https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/quantile graphics::boxplot でも quantile 関数を採 用しているが,オプションの変更はでき なさそう 6
SAS でパーセント点を計算するオプション <Syntax example> proc univariate data=dataX PCLTDEF = 5; var valueX ; output out=summary_pctl pctlpts = 25 50 75 pctlpre = P_; run; SASでは5種類のパーセント点の定義がある 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
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では5種類のパーセント点の定義がある ⇧平均化された経験分布関数 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
平均化された経験分布関数 𝑥𝑥(𝑗𝑗) , if 𝑔𝑔 = 0 𝑦𝑦 = � 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 > 0 1 𝑥𝑥 + 𝑥𝑥(𝑗𝑗+1) , 𝑦𝑦 = �2 (𝑗𝑗) 𝑥𝑥(𝑗𝑗+1) , if 𝑔𝑔 = 0 if 𝑔𝑔 > 0 Hyndman and Fan (1996) 𝑗𝑗 番目と 𝑗𝑗 + 1 番目の間を取る(平均化)することによって(分位点としての整 合性を保ちながら)経験分布関数の逆関数の連続性を考慮(補正)している 16
我々が習ったパーセント点は?(代表例) 1. 分位点がデータとデータの間に落ちると きは,中央値と同様にその2つのデータ の平均値を分位点とする ▸ 統計学入門(東京大学出版会) ▸ 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
まとめと結論 • 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). 医学への統計学【第3版】. 朝倉書店. P18. 20