スライド 1

Size: px

Start display at page:

Download "スライド 1"

ゆうりゅうしばもと
5 years ago
Views:

1 SASによる二項比率における正確な信頼区間の比較原茂恵美子 1) 武藤彬正 1) 宮島育哉 2) 榊原伊織 2) 1) 株式会社タクミインフォメーションテクノロジーシステム開発推進部 2) 株式会社タクミインフォメーションテクノロジービジネスソリューション部 Comparison of Five Exact Confidence Intervals for the Binomial Proportion using SAS Emiko Haramo 1) Akimasa Muto 1) Ikuya Miyajima 2) Iori Sakakibara 2) 1) System Development Department, Takumi Information Technology Inc. 2) Business Solutions Department, Takumi Information Technology Inc.

2 要旨 : SAS に導入されていない手法を含め二項比率の実際の信頼水準が常に名目の信頼水準を下回らない 5 つの正確な方法を紹介比較するとともにそれぞれの SAS プログラムを提供するキーワード : 二項比率信頼区間正確検定 2

3 1. 本研究の目的 2. 導入 3. 信頼区間の構成方法 3.1 Clopper-Pearson の信頼区間 3.2 一般化対数尤度比検定統計量に基づく正確な信頼区間 3.3 Score 検定統計量に基づく正確な信頼区間 3.4 Sterne 検定に基づく正確な信頼区間 3.5 Blaker 信頼区間 4. 結果 4.1 Coverage Probability 4.2 Expected Length 5. まとめ 3

4 1. 本研究の目的 2. 導入 3. 信頼区間の構成方法 3.1 Clopper-Pearson の信頼区間 3.2 一般化対数尤度比検定統計量に基づく正確な信頼区間 3.3 Score 検定統計量に基づく正確な信頼区間 3.4 Sterne 検定に基づく正確な信頼区間 3.5 Blaker 信頼区間 4. 結果 4.1 Coverage Probability 4.2 Expected Length 5. まとめ 4

5 背景二項比率の信頼区間に関する研究は古くからされており今現在も取り上げられている SAS に導入されている Wald 区間は計算が容易であるため二項比率の信頼区間として頻繁に使用されているが実際の信頼水準が名目の信頼水準を下回る事がある実際の信頼水準が名目の信頼水準を下回らない正確な信頼区間として SAS に導入されている Clopper-Pearson の正確な信頼区間が一般的に用いられているが小標本においてはとても保守的な方法である本研究の目的 1 SASに導入されていない方法を含む常に実際の信頼水準が名目の信頼水準を下回らない5つの正確な方法を紹介 2 小標本において信頼区間が真のパラメータを含む実際の確率や信頼区間の幅などを評価基準に比較及び検証 5

6 1. 本研究の目的 2. 導入 3. 信頼区間の構成方法 3.1 Clopper-Pearson の信頼区間 3.2 一般化対数尤度比検定統計量に基づく正確な信頼区間 3.3 Score 検定統計量に基づく正確な信頼区間 3.4 Sterne 検定に基づく正確な信頼区間 3.5 Blaker 信頼区間 4. 結果 4.1 Coverage Probability 4.2 Expected Length 5. まとめ 6

7 Notation 確率変数 X X~Bin(n, p) に従う. ( 但し Bin(n, p) はパラメータ n, p の二項分布とする ) 信頼区間の定義二項比率の100(1-α)% 信頼区間 P( P L p P U ) (1-α) (0 α 1 ) [P L, P U ] ( 但し P L : 信頼下限 P U : 信頼上限とし 1-α : 信頼水準 ( 名目の信頼水準 )) 7

8 評価指標とその定義 1. Coverage Probability ( 以下 C.P.) の定義信頼区間が真のパラメータ p をその中に含む実際の確率を C.P. という. 一般的に離散分布では, 名目の信頼水準に一致しない. P(P L p P ) 1 α C.P. は標本サイズ : 大 C.P.: 名目の信頼水準に近づく C.P. は名目の信頼水準に近い方が良い U C.P n=10, α= Expected Length ( 以下 E.L.) の定義信頼区間の幅の平均 E.L. E P U P L E.L. は幅が狭い方が良い p 8

9 1. 本研究の目的 2. 導入 3. 信頼区間の構成方法 3.1 Clopper-Pearson の信頼区間 3.2 一般化対数尤度比検定統計量に基づく正確な信頼区間 3.3 Score 検定統計量に基づく正確な信頼区間 3.4 Sterne 検定に基づく正確な信頼区間 3.5 Blaker 信頼区間 4. 結果 4.1 Coverage Probability 4.2 Expected Length 5. まとめ 9

10 Clopper-Pearson の信頼区間 2) / ( )/ ( 1 1 2) / ( 1) ( ) 1),2( 2(, 1),2 ( 2 x n x x x n F x n x x F x n x x 1), 1,2, ( n x Clopper andpearson(1934) パーセント点分布の上側の自由度 2 / 100 ), ( : 2) / (, F a b F b a 2 / ) (1 ) Pr( k n L n x k k L k n L p p C p x X 2 / ) (1 ) Pr( 0 k n U x k k U k n U p p C p x X

11 PGM Clopper-Pearson の信頼区間 SAS のプロシジャでも用意されている proc freq data = 入力データ ; weight 変数 1 / zeros ; table 変数 2 / binomial alpha = 有意水準 ; output out = 出力データ binomial ; run ; (1) FREQ プロシジャを用いる (2) TABLE ステートメントで BINOMIAL オプションを指定するデフォルト : 漸近標準誤差 Wald 信頼限界 Clopper-Pearson の信頼区間二項比率に関する漸近同等性検定 (3) OUTPUT ステートメントに BINOMIAL オプションを指定 11

12 3.2 一般化対数尤度比検定統計量に基づく正確な信頼区間帰無仮説 H0 : p P 0 とした時の一般化対数尤度比統計量は以下となる GLLR( P0 x) xln( x ( np0 )) ( n x)ln(( n x)/( n(1 P0 ))) 観測値 x が与えられたもとで GLLR(P 0 x i ) 集合を t とすると p-value はとなる従って p-value(x, P 0 ) Fleiss(2003) GLLR(P 0 x) を満たす X i の p value( x, P 0) Bin ( t n, P0 ) t α となる全ての P 0 が信頼区間となる 12

13 PGM 一般化対数尤度比検定統計量に基づく正確な信頼区間 (1) 一般化対数尤度比検定統計量 (GLLR) の算出 xとpの範囲で場合分け 1 0 < p < 1 のとき GLLR = ( x * log( x / ( n * p ) ) ) + ( ( n - x ) * log( ( n - x ) / ( n * ( 1 - p ) ) ) ) 2 x = n and 0 < p <= 1 のとき GLLR = ( x * log( x / ( n * p ) ) ) 3 1 と 2 を除き x = 0 and 0 <= p < 1 のとき GLLR = ( ( n - x ) * log( ( n - x ) / ( n * ( 1 - p ) ) ) ) 4 ( x = 0 and p = 1 ) または ( x = n and p = 0 ) のとき GLLR =. t が 0 から n の時の X i の GLLRi を計算 13

14 PGM 一般化対数尤度比検定統計量に基づく正確な信頼区間 (2) それぞれの t の時の二項比率を算出 if t = 0 then Bin = probbnml( p, n, t ) ; else Bin = probbnml( p, n, t ) - probbnml( p, n, t - 1 ) ; (3) GLLR(P 0 x i ) GLLR(P 0 x) の指定 if n( GLLR, GLLRi ) = 2 and GLLRi >= GLLR then output ; (4) (3) を満たす t の集合から p-value(x, P 0 ) α の抽出 proc means data = 入力データ nway noprint ; output out = 出力データ ( where = ( p-valueの変数 >= alpha ) ) sum = p-valueの変数 ; 14

15 3.3 Score 検定統計量に基づく正確な信頼区間帰無仮説 H : p 0 P 0 とした時の Score 統計量は以下となる 2 W( P0 x) = ( x-p0 ) / ( np0 (1-P0 ) ) Hirji(2006) 観測値 x が与えられたもとで W(P 0 x i ) W(P 0 x) を満たすX i の集合を t とすると p-value はとなる従って p-value(x, P 0 ) p value( x, P 0) Bin ( t n, P0 ) t α となる全ての P 0 が信頼区間となる 15

16 PGM Score 検定統計量に基づく正確な信頼区間 (1) Score 検定統計量 (SCORE) の算出 xとpの範囲で場合分け 1 0 =< x =< n and 0 < p < 1 のとき SCORE = ( ( x - ( n * p ) ) ** 2 ) / ( n * p * ( 1 - p ) ) 2 ( x = 0 and p = 0 ) または ( x = n and p = 1 ) のとき SCORE = 0 3 ( x = 0 and p = 1 ) または ( x = n and p = 0 ) のとき SCORE =. tが0 からnの時のX i のSCOREiを計算 SCOREとSCOREiを比較信頼区間算出 16

17 3.4 Sterne 検定に基づく正確な信頼区間両側有意水準 αの仮説検定において確率関数 f(x n, p) の高い値から足し合わせ初めて1- αを超えるところまでを採択域とするある観測値 xに対して採択域にあるpの集合を信頼区間とする ( ただし f(x n, p) はパラメータn, pの二項分布の確率関数とする ) Reiczigel (2003) n = 6, α = 0.1 X P i

PGM Sterne 検定に基づく正確な信頼区間 (1) 二項比率の確率関数の算出 Prob=PDF( 'binom', X, p, n ) ; 0.

18 PGM Sterne 検定に基づく正確な信頼区間 (1) 二項比率の確率関数の算出 Prob=PDF( 'binom', X, p, n ) ; 初めて有意水準を超える (2) 確率関数の合計を算出 Prob_AD + Prob ; (3) 確率関数の合計が有意水準を超えるポイントにフラグを立てる BProb_AD = lag1( Prob_AD ) ; if (BProb_AD < 0.9 ) then FLG = 0 ; 採択域 else FLG = 1 ; 作成したフラグから上限と下限を出力して信頼区間とする 18

19 3.5 Blaker 信頼区間 Sterne 信頼区間に似ている区間である Acceptability function ( 以下 ACC) を算出する ACC(x) - α 0 が採択域となる x n ACC ( x) = minf ( i, p), f ( j, p) i 0 j x Blaker (2000) 以下の表は ACC(x) - α にて採択域を示している n = 6, α = 0.1 X P i ACC(x)- 19

20 PGM Blaker 信頼区間 data work.out1 ; length STAT $ 100 ; STAT = "Blaker" ; do p_ = 0 to 1000 ; do x = 0 to &n ; p = p_ / 1000 ; if x = 0 then p1 = 1 ; else p1 = 1 - probbnml( p, &n, x - 1 ) ; p2 = probbnml( p, &n, x ) ; output ; end ; end ; run ; data work.out2 ; set work.out1; do u = &n to 0 by -1.0 ; px1 = probbnml( p, &n, u ) ; if px1 >= p1 then x1 = u ; if px1 >= ( 1 - p2 ) then x2 = u ; end ; if x1 = 0 then a1 = p1 ; else a1 = p1 + probbnml( p, &n, x1-1 ) ; a2 = p probbnml( p, &n, x2 ) ; accept = min( a1, a2 ) ; if ( accept - &alpha ) >= 0 ; run ; Blaker(2000) に記載されているプログラムから SAS 化 20

21 1. 本研究の目的 2. 導入 3. 信頼区間の構成方法 3.1 Clopper-Pearson の信頼区間 3.2 一般化対数尤度比検定統計量に基づく正確な信頼区間 3.3 Score 検定統計量に基づく正確な信頼区間 3.4 Sterne 検定に基づく正確な信頼区間 3.5 Blaker 信頼区間 4. 結果 4.1 Coverage Probability 4.2 Expected Length 5. まとめ 21

22 4.1 Coverage Probability C.P. が名目の信頼水準に近い方が良い信頼区間検証方法それぞれの正確な信頼区間による C.P. のグラフで比較設定 95% 信頼区間二項比率 : p = シミュレーション : 10 万回 22

23 4.1.1 Coverage Probability (α = 0.05, n = 5) Clopper-Pearson GLLR Score p Sterne p p Blaker p p Clopper-Pearson は明らかに保守的である GLLRとScore Blaker はpの値により高くなることがある 23

24 4.1.1 Coverage Probability (α = 0.05, n = 10) Clopper-Pearson GLLR Score p Sterne p p Blaker p p 全体的に名目の信頼水準に近づいた Clopper-Pearson はまだ保守的である 24

25 4.1.2 Coverage Probability (α = 0.05, p = 0.25, n = 5-95) 1.00 Exact Exact Clopper-Pearson GLLR Score n n n Sterne 1.00 Blaker n Clopper-Pearson Score Blakerはnの値によりばらつきがある GLLRとSterneはばらつきが少ない n 25

26 4.1.2 Coverage Probability (α = 0.05, p = 0.50, n = 5-95) 1.00 Clopper-Pearson Exact GLLR Exact Score n Sterne 1.00 n Blaker n n p =0.25 と同様の結果となった n n Clopper-Pearson Score Blaker は n の値によりばらつきがある GLLR と Sterne はばらつきが少ない n 26

27 4.2 Expected Length E.L. は幅が狭い方が良い信頼区間検証方法それぞれの正確な信頼区間による E.L. のグラフで比較設定 95% 信頼区間標本サイズ :n = 5, 10, 15, 20 二項比率 :p = 0.001~0.500 シミュレーション : 10 万回 27

28 4.2 Expected Length (α = 0.05, p = ) n = 5 n = 10 n = 15 n = 20 p p 28

29 4.2 Expected Length のまとめ標本数が大きくなるにつれ全体的に E.L. は小さくなり手法ごとの違いも小さくなった Clopper-Pearson はすべての n において最も幅が広く保守的な値となった n = 5 GLLR は n が 10 以上のとき p = 0 付近では最も幅が狭くなっているが p = 0.5 では広い n = 10 29

30 1. 本研究の目的 2. 導入 3. 信頼区間の構成方法 3.1 Clopper-Pearson の信頼区間 3.2 一般化対数尤度比検定統計量に基づく正確な信頼区間 3.3 Score 検定統計量に基づく正確な信頼区間 3.4 Sterne 検定に基づく正確な信頼区間 3.5 Blaker 信頼区間 4. 結果 4.1 Coverage Probability 4.2 Expected Length 5. まとめ 30

31 まとめ信頼区間の構成方法 Clopper-Pearson の信頼区間一般化対数尤度比検定統計量 (GLLR) に基づく正確な信頼区間 Score 検定統計量に基づく正確な信頼区間 Sterne 検定に基づく正確な信頼区間 Blaker 信頼区間 Coverage Probability 全体的に高い n の値によりばらつきがある p の値により高い n の値によるばらつきが少ない p の値により高い n の値によりばらつきがある名目の有意水準に近い n の値によるばらつきが少ない名目の有意水準に近い n の値によりばらつきがある Expected Length 明らかに幅が広い p の値により幅が広い狭い狭い狭い小標本において SAS に導入されていない Sterne 検定に基づく正確な信頼区間が比較的良い方法 31

32 マクロプログラム 32

33 %macro Five_Exact ( alpha =, n =, x = ) ; %* 1. Clopper-Pearson Confidence Interval *; data work.out1 ; i = 0 ; val = &n - &x ; output ; i = 1 ; val = &x ; output ; run ; proc freq data = work.out1 noprint ; weight val / zeros ; table i / binomial alpha = &alpha ; output out = work.out2 binomial ; run ; Freq プロシジャにて Clopper-Pearson の信頼区間を算出 data work.c_p ; set work.out2 ; length STAT $ 100 ; STAT = "Clopper-Pearson" ; Low = round( XL_BIN, ) ; Up = round( XU_BIN, ) ; keep STAT Low Up ; run ; 33

34 %* 2. Exact GLLR Confidence Interval *; data work.out1 ( drop = p_ ) ; length STAT $ 100 ; STAT = "GLLR" ; do p_ = 0 to 1000 by 1 ; p = p_ / 1000 ; if &x = 0 and 0 <= P < 1 then GLLR = ( ( &n - &x ) * log( ( &n - &x ) / ( &n * ( 1 - p ) ) ) ) ; else if &x = &n and 0 < P <= 1 then GLLR = ( &x * log( &x / ( &n * p ) ) ) ; else if 0 < P < 1 then GLLR = ( &x * log( &x / ( &n * p ) ) ) + ( ( &n - &x ) * log( ( &n - &x ) / ( &n * ( 1 - p ) ) ) ) ; else GLLR =. ; do t = 0 to &n ; if t = 0 and 0 <= P < 1 then GLLRi = ( ( &n - t ) * log( ( &n - t ) / ( &n * ( 1 - p ) ) ) ) ; else if t = &n and 0 < P <= 1 then GLLRi = ( t * log( t / ( &n * p ) ) ) ; else if 0 < P < 1 GLLR を算出 GLLRi を算出 then GLLRi = ( t * log( t / ( &n * p ) ) ) + ( ( &n - t ) * log( ( &n - t ) / ( &n * ( 1 - p ) ) ) ) ; else GLLRi =. ; 34

35 if t = 0 then Bin = probbnml( p, &n, t ) ; else Bin = probbnml( p, &n, t ) - probbnml( p, &n, t - 1 ) ; if n( GLLR, GLLRi ) = 2 and GLLRi >= GLLR then output ; end ; end ; run ; proc means data = work.out1 nway noprint ; class STAT p ; var Bin ; output out = work.out2 ( where = ( pval >= &alpha ) ) sum = pval ; run ; GLLRi >= GLLR の判定を行う proc means data = work.out2 nway noprint ; class STAT ; var p ; output out = work.gllr MIN = Low MAX = Up ; run ; p-value(x, P 0 ) >= の抽出 35

36 %* 3. Exact Score Confidence Interval *; data work.out1 ( drop = p_ ) ; length STAT $ 100 ; STAT = "Score" ; do p_ = 0 to 1000 by 1 ; p = p_ / 1000 ; if &x = 0 and P = 0 then SCORE = 0 ; else if 0 =< &x =< &n and 0 < P < 1 then SCORE = ( ( &x - ( &n * p ) ) ** 2 ) / ( &n * p * ( 1 - p ) ) ; else if &x = &n and P = 1 then SCORE = 0 ; else SCORE =. ; do t = 0 to &n ; if t = 0 and P = 0 then SCOREi = 0 ; else if 0 =< t =< &n and 0 < P < 1 then SCOREi = ( ( t - ( &n * p ) ) ** 2 ) / ( &n * p * ( 1 - p ) ) ; else if t = &n and P = 1 then SCOREi = 0 ; else SCOREi =. ; SCORE を算出 SCOREi を算出 36

37 if t = 0 then Bin = probbnml( p, &n, t ) ; else Bin = probbnml( p, &n, t ) - probbnml( p, &n, t - 1 ) ; if n( SCORE, SCOREi ) = 2 and SCOREi >= SCORE then output ; end ; end ; run ; proc means data = work.out1 nway noprint ; class STAT p ; var Bin ; output out = work.out2 ( where = ( pval >= &alpha ) ) sum = pval ; run ; SCOREi >= SCORE の判定を行う proc means data = work.out2 nway noprint ; class STAT ; var p ; output out = work.score MIN = Low MAX = Up ; run ; p-value(x, P 0 ) >= の抽出 37

38 %* 4. Sterne Confidence Interval *; data work.out1; length STAT $ 100 ; STAT = "Sterne" ; do x = 0 to &n ; do p_ = 0 to 1000 ; p = p_ / 1000 ; Prob = pdf( 'binom', x, p, &n ) ; output ; end ; end ; run ; 確率関数の算出 proc sort data = work.out1 ; by p descending Prob ; run ; 38

39 data work.out2 ; set work.out1 ; by p ; if first.p then do ; Prob_AD = 0 ; BProb_AD =. ; FLG = 0 ; end ; Prob_AD + prob ; BProb_AD = lag1( Prob_AD ) ; if ( BProb_AD < (1 - &alpha) ) then FLG = 0 ; else FLG = 1 ; run ; 採択域を算出しフラグを立てる proc means data = work.out2 nway noprint ; where x = &x and FLG = 0 ; class STAT ; var p ; output out = work.sterne MIN = Low MAX = Up ; run ; 39

40 %* 5. Blaker Confidence Interval *; data work.out1 ; length STAT $ 100 ; STAT = "Blaker" ; do p_ = 0 to 1000 ; do x = 0 to &n ; p = p_ / 1000 ; if x = 0 then p1 = 1 ; else p1 = 1 - probbnml( p, &n, x - 1 ) ; p2 = probbnml( p, &n, x ) ; output ; end ; end ; run ; Blaker(2000) に記載されているプログラムから SAS 化 data work.out2 ; set work.out1; do u = &n to 0 by -1.0 ; if u >= 0.0 ; px1 = probbnml( p, &n, u ) ; if px1 >= p1 then x1 = u ; if px1 >= ( 1 - p2 ) then x2 = u ; end ; if x1 = 0 then a1 = p1 ; else a1 = p1 + probbnml( p, &n, x1-1 ) ; a2 = p probbnml( p, &n, x2 ) ; accept = min( a1, a2 ) ; if ( accept - &alpha ) >= 0 ; run ; proc means data = work.out2 nway noprint ; where x = &x ; class STAT ; var p ; output out = work.blaker MIN = Low MAX = Up ; run ; 40

41 Title "Exact Confidence Interval ( N = &N, X = &x )" ; data work.outds ; set work.c_p work.gllr work.score work.sterne work.blaker ; run ; 全ての方法の信頼区間を合わせる proc print label noobs ; var STAT Low Up ; format Low Up 8.3 ; label STAT = Exact Confidence Interval Low = Lower Confidence Interval Up = Upper Confidence Interval ; run ; %mend ; %Five_Exact( alpha = 0.05, n = 10, x = 3 ) ; 41

42 出力結果例 42

43 Reference 1. Agresti, A. and Coull, B. A. (1998). Approximate is better than exact for interval estimation of binomial proportions. The American Statistician 52, Newcombe, R.G. (1998). Two-sided confidence intervals for the single proportion: comparison of seven methods. Statistics in Medicine 17, Clopper, C.J. and Pearson, E.S. (1934). The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika 26, Blaker, H. (2000). Confidence curves and improved exact confidence intervals for discrete distributions. The Canadian Journal of Statistics 28, Hirji, K.(2006). Exact analysis of discrete data. Chapman and Hall/CRC: New York. pp Reiczigel, J. (2003). Condence intervals for the binomial parameter: some new considerations. Statistics in Medicine 22,

44 ご清聴ありがとうございました 44

スライド 1

スライド 1 SAS による二項比率の差の非劣性検定の正確な方法について武藤彬正宮島育哉榊原伊織株式会社タクミインフォメーションテクノロジー Eact method of non-inferiority test for two binomial proportions using SAS Akimasa Muto Ikuya Miyajima Iori Sakakibara Takumi Information