データステップでの統計検定の再現

2024年3月22日 大阪SAS勉強会 データステップを用いた統計検定の再現 折村 奈美

目次 1.SASプロシジャと検定 2.Studentのt検定 3.カイ二乗検定 4.Wilcoxonの符号付順位和検定 5.まとめ Copyright©EPS All rights reserved. 2

SASプロシジャと検定 Copyright©EPS All rights reserved. 3

SASプロシジャと検定 SASのプロシジャは検定を行う上で非常に便利である 短いコードを書くだけで 検定を実行してくれる https://go.documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/statug/statug_ttest_examples01.htm 検定の中身を理解しなくてもプロシジャを使えば検定ができてしまう ⇒難しい解析に対応できないかも? 検定を理解してさまざまな解析への対応力を付ける目的で データステップで検定を再現した ◎本発表内のルール 以下のプロシジャは使用してもよい ・proc sort ・proc transpose ・proc summary ・proc rank Copyright©EPS All rights reserved. t検定、カイ二乗検定、 Wilcoxon符号付順位和検定を再現する 4

• https://go.documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/statug/statug_ttest_examples01.htm

Studentのt検定 Copyright©EPS All rights reserved. 5

Studentのt検定 性別間の体重を比較する proc ttest data=sashelp.class; class sex; var weight; run; t値の公式 𝑡= 𝑥ҧ1 − 𝑥ҧ 2 1 1 𝑆2 𝑛 + 𝑛 1 2 𝑆 2（プールした分散）の公式 (𝑛1−1)𝑆12 + (𝑛2−1)𝑆22 2 𝑆 = 𝑛1 + 𝑛2 − 2 これを再現する Copyright©EPS All rights reserved. 6

Studentのt検定 ①S 2とt値の計算に必要な n,平均,不偏分散を求める proc summary data=sashelp.class nway; var weight; class sex; output out=wk1 n= mean= var=/ autoname; run; proc transpose data=wk1 out=n(drop=_:) prefix=n; var Weight_N; run; proc transpose data=wk1 out=mean(drop=_:) prefix=mean; var Weight_Mean; run; proc transpose data=wk1 out=var(drop=_:) prefix=var; var Weight_Var; run; data wk2; set n; set mean; set var; run; Copyright©EPS All rights reserved. 7

Studentのt検定 ②公式に当てはめて自由度、S2 、t値を算出 p値も算出する t値の公式 𝑡= data wk3; set wk2; /*自由度*/ df=n1+n2-2; /*𝑆 2*/ s2=divide(((n1-1)*var1+(n2-1)*var2), df); 𝑥ҧ1 − 𝑥ҧ 2 𝑆2 1 1 + 𝑛1 𝑛2 𝑆 2の公式 𝑆 2 (𝑛1−1)𝑆12 + (𝑛2−1)𝑆22 = 𝑛1 + 𝑛2 − 2 /*t値*/ t=divide((mean1-mean2), sqrt(s2*(1/n1+1/n2))); /*p値*/ if t<0 then p=2*probt(t, df); else p=2*(1-probt(t, df)); run; Copyright©EPS All rights reserved. 一致！ 8

カイ二乗検定 Copyright©EPS All rights reserved. 9

カイ二乗検定 (独立性の検定) TREATとRESPが独立であるかどうか検定 data data2; TREAT="N";RESP="N";_FREQ_=5;output; TREAT="N";RESP="Y";_FREQ_=15;output; TREAT="Y";RESP="N";_FREQ_=8;output; TREAT="Y";RESP="Y";_FREQ_=12;output; run; proc freq data=data2; tables TREAT*RESP/ CHISQ; weight _FREQ_/zeros; run; カイ二乗値の公式 𝑚 𝑛 2 (𝑥 − 𝐸 ) 𝑖𝑗 𝑖𝑗 χ2 = ෍ ෍ 𝐸𝑖𝑗 𝑖=1 𝑗=1 自由度(𝑚 − 1) × (𝑛 − 1)のカイ二乗分布に従う 𝑥𝑖𝑗 : 観測度数, 𝐸𝑖𝑗 : 期待度数 𝑚: 縦の行数, 𝑛: 横の列数 Copyright©EPS All rights reserved. これを再現する 10

カイ二乗検定 (独立性の検定) ①クロス表の合計を計算する proc summary data=data2 nway; var _FREQ_; class TREAT; output out=wk1_treat(drop=_FREQ_) sum=total_treat; run; proc summary data=data2 nway; var _FREQ_; class RESP; output out=wk1_resp(drop=_FREQ_) sum=total_resp; run; proc summary data=data2 nway; var _FREQ_; output out=_wk1_total(keep=total) sum=total; run; data wk1_total; set _wk1_total; _TYPE_=1; run; Copyright©EPS All rights reserved. 11

カイ二乗検定 (独立性の検定) ②①で計算した合計を結合する data wk2_treat; merge data2 wk1_treat; by TREAT; run; proc sort data=wk2_treat; by RESP; run; data wk2_resp; merge wk2_treat wk1_resp; by RESP; run; data wk2(drop=_TYPE_); merge wk2_resp wk1_total; by _TYPE_; run; 上のクロス表をデータセットにできた proc sort data=wk2; by TREAT RESP; run; Copyright©EPS All rights reserved. 12

カイ二乗検定 (独立性の検定) ③期待度数および観測度数と期待度数との差を求め、χ2値を計算する data wk3; set wk2; /*期待度数*/ ef=divide((total_treat*total_resp), total); /*観測度数と期待度数の差*/ diff=_FREQ_-ef; /*χ2値*/ _x2=divide(diff**2, ef); run; proc summary data=wk3 nway; var _x2; output out=wk4 sum=x2; run; /*p値*/ data wk5; set wk4; prob=1-probchi(x2, 1); run; Copyright©EPS All rights reserved. 𝑚 χ2 𝑛 (𝑥𝑖𝑗 − 𝐸𝑖𝑗 )2 = ෍෍ 𝐸𝑖𝑗 𝑖=1 𝑗=1 𝑥𝑖𝑗 : 観測度数, 𝐸𝑖𝑗 : 期待度数 一致！ 13

Wilcoxon符号付順位和検定 Copyright©EPS All rights reserved. 14

Wilcoxonの符号付順位和検定 BEFOREとAFTERに差があるかを調べる proc univariate data=data1; var DIFF; run; S値の公式 𝑆 = 𝑅+ − 𝑛𝑡 (𝑛𝑡 + 1) 4 𝑅+ : 正の符号付順位の合計 𝑛𝑡 : n数 *DIFF=0のデータは除く これを再現する Copyright©EPS All rights reserved. 15

Wilcoxonの符号付順位和検定 ①差の符号付順位を求める /*DIFFの絶対値*/ data wk1; set data1; DIFF_abs=abs(DIFF); run; /*DIFFの絶対値を用いて順位付け*/ proc rank data=wk1 out=wk2 ties=mean; where DIFF_abs^=0; var DIFF_abs; RANKS rank; run; /*順位に符号を付ける*/ data wk3; set wk2; if 0<DIFF then signed_rank=rank; else if DIFF<0 then signed_rank=-rank; run; Copyright©EPS All rights reserved. 16

Wilcoxonの符号付順位和検定 ②n数、𝑹+、S値を計算 /*n数*/ proc summary data=wk3 nway; var signed_rank; output out=wk4_n n=/ autoname; run; /*𝑅+ */ proc summary data=wk3 nway; where 0<signed_rank; var signed_rank; output out=wk4_sum sum=/ autoname; run; S値の公式 𝑆 = 𝑅+ − 𝑛𝑡 (𝑛𝑡 + 1) 4 𝑅+ : 正の符号付順位の合計 𝑛𝑡 : n数 *DIFF=0のデータは除く data wk4; set wk4_n; set wk4_sum; /*S値*/ S=signed_rank_Sum -divide(signed_rank_N*(signed_rank_N+1), 4); run; Copyright©EPS All rights reserved. S値一致！ 17

Wilcoxonの符号付順位和検定 ④p値を求める N≥25のときp値は以下のように計算する T+ ：符号付順位が正の合計順位 T- ：符号付順位が負の合計順位 T = min(T+ , T- ) 𝑍= 𝑛𝑡 (𝑛𝑡 + 1) 4 𝑛𝑡 (𝑛𝑡 + 1)(2𝑛𝑡 + 1) 24 𝑇− 右の式で標準化してp値を求める /*T+を計算*/ proc summary data=wk3 nway; where 0<signed_rank; var rank; output out=wk5_plus sum=rank_plus; run; /*T-を計算*/ proc summary data=wk3 nway; where signed_rank<0; var rank; output out=wk5D_minus sum=rank_minus; run; Copyright©EPS All rights reserved. 18

Wilcoxonの符号付順位和検定 ④T, z値を計算し、p値を求める data wk6; set wk5_plus; set wk5_minus; set wk4(keep=signed_rank_n); 𝑍= 𝑛𝑡 (𝑛𝑡 + 1) 4 𝑛𝑡 (𝑛𝑡 + 1)(2𝑛𝑡 + 1) 24 𝑇− T=min(rank_plus, rank_minus); v=sqrt(divide(signed_rank_n*(signed_rank_n+1)*(2*signed_rank_n+1), 24)); e=T-divide(signed_rank_n*(signed_rank_n+1), 4); z=divide(e, v); prob=probnorm(z)*2; run; p値が一致しない t分布で近似して計算している? Tests for Location :: Base SAS(R) 9.4 Procedures Guide: Statistical Procedures, Second Edition Copyright©EPS All rights reserved. 19

• https://support.sas.com/documentation/cdl/en/procstat/66703/HTML/default/viewer.htm

Wilcoxonの符号付順位和検定 N≥20のときp値は以下のように計算する 𝑆 1 𝑛 𝑛 + 1 2𝑛 + 1 − 𝑑0 𝑑0 + 1 2𝑑0 + 1 − σ𝑖>0 𝑑𝑖 (𝑑𝑖 + 1)(𝑑𝑖 − 1) 2 𝑛 × 𝑉𝑎𝑟 𝑆 − 𝑉𝑎𝑟 𝑆 = 24 𝑛−1 n: DIFF=0を含むn数 𝑑0 : DIFF=0のn数 𝑑𝑖 : DIFF≠0のn数 上の式で求められた統計量tが自由度n-1のt分布で近似される 𝑡= 𝑆2 /*nを計算*/ proc summary data=data1 nway; var DIFF; output out=wk_n n=n; run; /*𝑑0を計算*/ proc summary data=data1 nway; where DIFF=0; var DIFF; output out=wk_d0 n=d0; run; Copyright©EPS All rights reserved. 20

Wilcoxonの符号付順位和検定 𝑆 1 𝑛 𝑛 + 1 2𝑛 + 1 − 𝑑 𝑑 + 1 2𝑑 + 1 − σ𝑖>0 𝑑𝑖 (𝑑𝑖 + 1)(𝑑𝑖 − 1) 0 0 0 2 𝑛 × 𝑉𝑎𝑟 𝑆 − 𝑆 2 𝑉𝑎𝑟 𝑆 = 24 𝑛−1 n: DIFF=0を含むn数 𝑑0 : DIFF=0のn数 𝑑𝑖 : DIFF≠0のn数 上の式で求められた統計量tが自由度n-1のt分布で近似される 𝑡= data wk5; set wk4; set wk_n; set wk_d0; /*Var(S)*/ var=divide(n*(n+1)*(2*n+1)-d0*(d0+1)*(2*d0+1) -divide(signed_rank_n*(signed_rank_n+1)*(signed_rank_n-1), 2), 24); /*t値*/ _t=sqrt(divide(n*var-S**2, n-1)); t=divide(S, _t); /*自由度*/ df=n-1; /*p値*/ prob=probt(t, df)*2; run; Copyright©EPS All rights reserved. p値一致！ 21

まとめ • Studentのt検定、カイ二乗検定、Wilcoxon符号付順位和検定をデータステップで 再現し、理解を深めることができた • 検定の方法は一つではない 調べるサイトによって記載されている方法が違って結果も異なる • データステップでの検定の再現は、難易度の高い解析に対応する力を養うために有用である と感じた Copyright©EPS All rights reserved. 22

参考文献 [1] 株式会社アイスタット, t検定の概要と手順・結果の見方, 2024-03-07確認 https://istat.co.jp/sk_commentary/t-test/t-test [2] GMO, カイ二乗検定とは？検定手法を解説-GMOリサーチ, 2024-02-15確認 https://gmo-research.jp/research-column/chi-square-test [3] Social Survey Research Information Co., Ltd., 6-4. ノンパラメトリック検定-対応のある2標本の差の検定|統計学の時 間|統計WEB, 2024-02-15確認 https://bellcurve.jp/statistics/course/26159.html [4] SAS Institute, SAS Help Center: Tests for Location, 2024-02-15確認 https://go.documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/procstat/procstat_univariate_details17.htm [5] jmp, jmp Help, Wilcoxonの符号付順位検定の統計的詳細, 2024-02-15確認 https://www.jmp.com/support/help/ja/17.2/index.shtml#page/jmp/statistical-details-for-the-wilcoxon-signedrank-test.shtml Copyright©EPS All rights reserved. 23

• https://istat.co.jp/sk_commentary/t-test/t-test
• https://gmo-research.jp/research-column/chi-square-test
• https://bellcurve.jp/statistics/course/26159.html
• https://go.documentation.sas.com/doc/en/pgmsascdc/9.4_3.5/procstat/procstat_univariate_details17.htm
• https://www.jmp.com/support/help/ja/17.2/index.shtml