Comments
Description
Transcript
“n ≥ 30” ルールを捨てるときがきた
“n ≥ 30” ルールを捨てるときがきた 2012 年 6 月 10 日版 Tim Hesterberg∗ 日本語訳 荒木 孝治 † 2013 年 12 月 10 日 1 はじめに 正規近似と t 統計量に基づく信頼区間や仮説検定は，統計学で非常に良く利用されている．これら は，ある正則条件が成立するとき，標本平均や回帰係数といった推定量の分布はサンプルの大きさが大 きくなると正規分布で近似することができ，対応する t 統計量は t 分布で近似できるという漸近的な結 果に基づいている． 有限サンプルに対しては，共通の大まかなルールに基づいて行われる．たとえば，独立で同一の分布 （i.i.d）に従うデータの平均値に対しては，サンプルの大きさが 30 以上あり，サンプルがあまり歪んで いない場合，正規分布に基づく推測を行ってよいというものである． しかし，“あまり歪んでいない” とはどういうことか．平均以外の統計量の場合には，どのような基 準を適用すればよいのか．そして，サンプルの大きさが小さい，あるいは，サンプルがかなり歪んでい るときにどうするのか．とにかく正規分布に基づく推測を行っておこう！（例外はあるとしても）そし て，ロジスティック回帰といった他の状況においては，どのような診断測度を利用すればよいのか？ 本稿で私は次を主張する． “n ≥ 30 かつ，あまり歪んでいない” という規準を，ブートストラップに基づくより効果的な診 断でおきかえるべきである． ブートストラップに基づく診断手法は適用が簡単であり， t 分布に従う推測がどれほど不正確かということを示すと，多くの統計学者は驚くだろう． t 分布に基づく推測より良い代替法がある． また，次に示す 2 つの関連する点を指摘する． 1000 個のブートストラップサンプルは高品質の診断を行うには不十分であり， 大きなサンプルの場合より良い推測が可能であるが，小さなサンプルの場合には，さらなる研究 が必要である． ∗ † Google, 651 N. 34th St., Seattle WA 98103 翻訳に関する意見等は [email protected] まで． 1 2 ブートストラップ診断 ここでは，ブートストラップ法を簡単に説明しておく．より詳細な入門には，Hesterberg et al.(2003)， Efron and Tibshirani (1993)，Davison and Hinkley (1997) を参照． 母集団分布 F（多変量であってもよい）からの独立で同一の分布に従う（i.i.d）サンプルを X1 , . . . , Xn とし、パラメータ θ を推定したいとする．θb を推定量とし，経験分布関数 Fbn （各データポイントに確 率 1/n を与える関数）のみを通じてデータに依存しているとする． 通常のノンパラメトリックなブートストラップでは，経験分布からサンプルを抽出，つまり， X1∗ , . . . , Xn∗ から復元抽出を行い，これに基づいて対応する統計量 θb∗ を計算する．これを多数回，たと ∗ を得る．これを用いてブー えば，B = 1000 回繰り返すと，B 個のブートストラップ統計量 θb1∗ , . . . , θbB トストラップ分布を構成し，標準誤差や信頼区間の計算や診断を行うことができる． 2.1 TV の平均の例 学生の Barrett Rogers が，ケーブル TV の基本契約と拡張契約（追加コストが必要）の 30 分当たり のコマーシャル時間（分）を調べた例を取り上げる．基本チャネルのコマーシャルの平均時間は 9.21 分で，これは 10 回の観測データに基づく（学生は貧しいので，TV のランダムな 30 分を 20 回のみ観 測できた）．データは，基本チャネルで 7.0, 10.0, 10.6, 10.2, 8.6, 7.6, 8.2, 10.4, 11.0, 8.5，拡張チャネル で 3.4, 7.8, 9.4, 4.7, 5.4, 7.6, 5.0, 8.0, 7.8, 9.6 であった．基本チャネルの時間の平均のブートストラップ 分布を図 1 に示す． 図 1: TV データに対するブートストラップ分布．1 行目は基本テレビのコマーシャル時間のブー トストラップ平均のヒストグラムと正規 QQ プロット．2 行目は，基本および拡張チャネルの平均 の差のヒストグラムと正規 QQ プロット． 2 2 標本問題の場合，2 つのサンプルから独立にサンプルを抽出し，ブートストラップサンプルの各ペ アに対して，関心のある統計量，例えば，平均の差やハザード比を計算する．基本と拡張の間の平均の 差のブートストラップ分布を図 1 に示す． この場合，サンプルの大きさは非常に小さいが，ブートストラップ法により，1 標本平均と平均の差 に対する標本分布はほぼ正規分布であることがわかる． 2.2 Verizon 平均の例 次の例では，ブートストラップは異なる状況を示す．データを，図 2 に示す．“ILEC” サンプルは， 1664 個のデータである（平均 8.4 時間の修理時間）; より小さな “CLEC” サンプルは，23 個のデータ で，平均 16.5 である（Hesterberg et al., 2003）．これらは顧客の 2 つのグループに対する修理時間で あり，平均が片側の有意水準 0.01 で異なるかどうかを知りたい． ブートストラップ分布を図 2 の下側に示す．大きなサンプル（n = 1664）に対する平均のブートス トラップ分布は，ほぼ正規分布と考えることができる．しかし，小さなサンプルに対しては，かなり歪 みがある． 歪みがこれくらいになると，問題となる．これは，データを QQ プロットで見るという多くの読者の 直観と矛盾する．このブートストラップ分布は，生データではなく，標本分布に対応する．これは，中 心極限定理が適用された状況の後なので，正規分布からの逸脱は，推定における誤差に変換されてい る．歪みが信頼区間に及ぼす影響の大きさの定量化は第 3 節で行う．そこでは，ブートストラップ t 分 布の文脈で評価する．その前に，平均以外の統計量の例を見ておく． 2.3 野生動物の肉データの回帰の例 Brashares et al. (2004) は，Ghana の野生動物の肉（Bushmeat）の取得を目的とする狩りによる野 生生物の減少と魚の供給との間の関係を論じた．図 3 は，一人当たりの魚の供給量，全生物量の推定値 （国立公園にいる 30 種の個体数から推定）の 30 年のデータと，魚の供給量と生物猟の相対変化との散 布図である．魚の供給量が減少すると，生物量が大きく減少していることは明確である．これ（および 他の証拠もサポートしているが）は，魚の供給量が少ないとき，野生動物の狩りが支配的となることを 示している． 図 3 の左下は，20 個のブートストラップサンプルの回帰直線を示す．関心のある量の 1 つに x 軸切 片がある．これは，野生生物の減少が 0 となる魚の供給量の推定値を与える．右下の図は，x 軸切片の ブートストラップ分布を示す．分布は強い正の歪みを持ち，正規近似は使えない． 2.4 Kyphosis ロジスティック回帰の例 Kyphosis（脊柱後弯症）データセット（Chambers and Hastie, 1992）は 4 変数に対する 81 の観測 値から構成される．目的変数 “Kyphosis” は術後に変形が生じたかどうかという 2 値を取る変数で，説 明変数は，Age（年齢），Number（手術を行った脊椎骨の数），Start（手術を行った最初の脊椎骨の番 号）である. Kyphosis を共変量に回帰させるロジスティック回帰を実行する． 図 4 に，Start（共変量で最も効果があるもの）に対する Kyphosis のひまわりプロットと，ロジス ティック回帰からの予測値（Age と Number については，それらのメディアンに固定している）を追加 したものを示す．右上の図は，この予測のための 20 本のブートストラップ曲線である．かなり変動し ていることは明白である．関心のある量の一つは，共変量の固定値に対する予測値の標本分布である； ここでは，Start の値が大きくなると，他の変数を中央値に固定したとき，予測のブートストラップ分 3 図 2: Verizon データセットの修理時間．ILEC と CLEC グループ（それぞれデータ数は， n = 1664，n = 23）．データはトップパネル．各グループの平均のブートストラップ分布は下部の パネル． 布は負の値を取ることはなく，強く正の歪みを持つ． 下部の 2 枚のパネルは，Age と Start に対する回帰係数のブートストラップ分布を示す；ブートスト ラップ分布は強く歪んでいるため，正規近似は適切でない． ある統計パッケージ（S-PLUS）からのロジスティック回帰の出力には賞賛に値する制限がある．そ れは，係数の推定値，標準誤差と t 統計量を与えるが，t 統計量に対応する P 値を与えないというもの である．t 統計が t 分布に従わないので，これは適切な所作である．残念なことに，すべてのパッケー ジがそのような制約を与えているわけではない． これらの例，および，標本平均以外の統計統計量の他の例における共通のテーマは，標本分布が本質 4 図 3: 30 年間の魚の供給量と野生の生物量．下のパネルはブートストラップ直線と x 軸切片の ブートストラップ分布． 的に歪んでいるということである．この場合，平均の場合以上に中心極限定理を信頼するのは慎重でな ければならない． 3 ブートストラップ t t 統計量のブートストラップ分布の話に戻る．古典的な統計理論では，正規分布を仮定するとき，標 √ 本平均 X と標本標準偏差 s は独立であり，t 統計量 t = (X − µ)/(s/ n) は t 分布に従う．ブートスト ラップ法を用いて，実際の分布がどれくらい t 分布に近いかを診断することができる．一般に，θb を推 定量，sθb を θb の標準誤差とし， t = (θb − θ)/sθb (1) を t 統計量とする．これのブートストラップ版は，次のようになる． b ∗ t∗ = (θb∗ − θ)/s θb (2) ここで s∗b は，ブートストラップ標本から求めた標準誤差である． θ 図 5 は Verizon データセットの t 統計量に関連するブートストラップ診断の結果である．左上の図 ∗ は，大きな (n = 1664) ILEC データセットに対する X 対 s∗ の散布図である．データの歪みのため， 5 図 4: Kyphosis データとブートストラップ分布．最初のパネルは応答変数の Kyphosis と 3 つの 共編量の中で最も関係の強い Start とのひまわりプロットにロジスティック回帰による予測値（他 の共変量はメジアン値で固定）を記入したもの．2 つ目のパネルは，20 ブートストラップサンプ ルからの予測値を示している．下のパネルは，ロジスティック回帰における 4 つの係数のうちの 2 つのブートストラップ分布である． これらの量は独立ではなく，強い相関を持つ．なぜなら，比較的大きな平均を持つブートストラップサ ンプルは長い右の裾からの観測値を含む傾向にあり，そのため，標準偏差が大きくなる傾向にあるから である．次の 2 つのパネルは，t 統計量のブートストラップ分布のヒストグラムである．小さな CLEC に対する分布は，強い負の歪みを持ち，二山形である（大きな観測値の出現したり，出現しなかったり したために）． 大きな ILEC サンプルの分布はより正規分布に近いが，少し負の歪みを持つ，最後のパネルは，ILEC データセットのブートストラップ分布に対する正規 QQ プロットである．これも，少し負の歪みを持っ ている． 平均の分布は正の歪みを持つが，t 統計量は負の歪みを持つことに注意．これは，平均が小さいとき， 式 (2) の分母は，負になり，分子は比較的 0 に近くなる傾向になるので，結果は相対的に大きな負の値 を取ることになる． 3.1 数値による診断 図より，n = 1664 に対しては，ブートストラップ t 分布は少し歪んでいるだけであることがわかる． しかし，最後のパネルは，図による確認により私たちが感じる安心感を完全に覆してしまう数値結果を 示している．ブートストラップ t 分布が真に t 分布なら，−t1663,.025 を下回ったり t1663,.025 を上回っ たりするブートストラップ標本は 2.5 ％であることが期待できる．しかし，実際の比率は，3.94 ％と 6 図 5: Verizon データに対する t 統計量のブートストラップ分布．上左パネルは，大きな（n = 1664） ∗ データセット ILEC の X vs s∗ の散布図．他のパネルは，CLEC と ILEC データセットに対す る t 統計量のブートスタラップ分布のヒストグラムと正規 QQ プロットである． 2.12 ％である． 危険率が 3.94 ％というのは，名目の 2.5 ％より 58 ％大きい1 ．これほどの大きさの誤差は許容でき ない．このような誤差は，財務システムを破壊する可能性がある．（うわー．） もし 1664 個の観測数でも t 推測にとって充分ではないのなら，古い t1663,.025 ルールを窓から捨てる ことができる． しかし，状況はそんなに悪いのか？ 4 t 統計量の収束と他の手法 非常に良くない．母集団分布が歪んでいるとき，t 推定が満足できる程度によくなるには永遠の時間 が必要である．中心極限定理は，地理学的なスケールでしか機能しない． この 2 つの要素である “満足できる程度” と歪度を考えよう． “満足できる正確さを持つ信頼区間” を，実際の危険率が両側で名目値の 10 ％以内であると定義す る．95 ％信頼区間では，名目的には各側で 2.5 ％のミスが生じる．すると，満足できる信頼区間は，名 目値の 90 ％から 110 ％の間でミスを生じる．これは各側で 2.25 ％から 2.75 ％の危険率に相当する． しかし重要なのは，各側における危険率であり，全体のそれではないことに注意．実践的な統計問題 で，真に両側であることはほとんどない．片側でのミスは，もう片側でのミスとは異なる効果を持つと いうのが一般的である．さらに，不均衡な区間（片側では 4 ％，もう片側では 1 ％といったもの）は真 1 訳注．(3.94 − 2.5)/2.5 × 100 7 のパラメータもそうであるといった偏った印象を与える．両側の誤差が相殺されることはない．全体の 正確さは，各サイドの絶対誤差の和で計測すべきである． b と t の分布の正規分布と t 分布への収束は， 歪度に戻ろう．母集団分布が歪んでいるとき，X √ O(1/ n) という比較的遅いオーダーになる． 議論のために，指数分布を考える．なぜなら，それはよく知られており，シミュレーションに便利だ からである．指数分布は非常に歪んでいるように見えるが，その歪度は 2 であり，それほど大きくな い．指数分布の右裾は，指数関数的に急激に減少する．分母の自由度が小さい F 分布といった多項式 の裾を持つ分布は，正規分布に近いように見えるが，高い歪度を持つ．実際には，多くの分布は，右裾 が多項式的な振る舞いを示す． 図 6 は，母集団分布が指数分布のときの，t 信頼区間と他の手法とのシミュレーション結果であ る．最初のパネルは，t 信頼区間の各側，および Johnson の歪度調整 t 統計量（Johnson, 1978）， x+ √s (tα n + s √3 (1 6 n + 2t2α )，ここで，s3 は標本歪度，の収束に関する一般的な結果を示す．Johnson 調整区間は，t 区間と比べて早く名目値に収束する．2 番目のパネルは尺度変換バージョンで，調整な しの t 統計量はかなり正確に収束するが，それは n が 5000 を超えてからである．3 番目のパネルは， Johnson の手法はかなり正確な値に収束するが，それは n = 220 のあたりでそれほど早いとは言え ないが，5000 よりはかなりよい．最後のパネルは，これらおよび追加的な 2 つの手法，つまり，ブー トストラップ t 区間（Efron, 1981, 1982; Efron and Tibshirani, 1993）と ABC 区間（DiCiccio and Efron, 1992），ブートストラップ BCa 区間への非サンプリング近似（Efron, 1987）の収束を示す．こ れらおよびいくつかの他のブートストラップ，非ブートストラップ信頼区間は，“2 次の精度”，つまり， O(1/n) のスピードで収束する（Efron and Tibshirani, 1993）．これに対して t 分布による区間は，“1 次の精度” しか持たない． t による信頼区間は正しい信頼率に収束するが，そのスピードは非常に遅い． 5 ブートストラップ標本の大きさ–サンプルの大きさ 64000 の問題 Verizon データのブートストラップ t 分布は，大きさ 104 のブートストラップ標本に基づく．片側の √ 危険率は，0.0394 であると推定される; その推定値の標準誤差は 0.0394 ∗ (1 − 0.0394)/104 = 0.0019 なので，0.0394 は名目値から 7 標準誤差以上離れている2 ． 高度に正確な数値診断に必要なブートストラップ標本の大きさは，伝統的な推奨値よりはるかに大き いことに注意しよう．Efron and Tibshirani (1993) は B = 200 を提案した．図 6 より，t 分布の区間 の信頼率対するシミュレーション結果は，標準誤差の推定には B = 25 以下で充分であり，信頼区間に 対しては，B = 1000 で十分であることを示している． 危険率が 2.5 ％で，それの 10 分の 1 の “満足のできる正確さ” の範囲内での高精度が必要である状 況では，正確さのマージンの 4 分の 1 に標準誤差を減らしたいだろう（0.025/10/4）．この場合，およ そ B = 64000 個のブートストラップ標本の大きさが必要となる．実際には，64, 000 ではなく 62, 400 といった感じである． 6 要約 ルール “n ≥ 30” はコンピュータのない時代の名残に過ぎない．今や，それを捨て去る時代である． 中心極限定理は地質学的時間スケールでしか機能しない． 2 訳注．(0.0394 − 0.025)/0.0019 = 7.58 8 図 6: t 信頼区間とその他の手法の収束に関するシミュレーション結果 なぜこうした問題が知られていないのか．第 1 の理由として，診断にブーツストラップ法を利用する ことを知らないからであると推測する． 別の理由もある．最大の理由は計算に関するものである．1981 年は，ブートストラップの初期の時 代であり，Johnson の修正法が開発された少し後の頃であるが，この頃に図 6 を作成するために必要な 計算時間は，20,000 時間くらいであった．今でも，64,000 個のブートストラップサンプルは多いが，日 常的に実行するには大きすぎるというほどではない．1000 個でも近似的に診断を行うのに十分である． 第 2 に，私たちは良い代替法を必要とする．単純平均の場合，サンプルの大きさがある程度あると Johnson の手法はかなり良いが，サンプルの大きさが小さいとき，歪度を正確に推定できないので良く ない．おそらく必要なのは semi-Bayesian アプローチであり，ここでは推定した歪度を 0 に向けて縮小 する．サンプルの大きさが小さいときには，より大きく縮小する． 比較のために，t 信頼区間は，poin mass を 0 に持つ歪んだ事前分布を持つベイジアンである． 平均以外の統計量に対して，ABC 区間や他のブートストラップおよび非ブートストラップ区間 は，2 次の正確性を持つが，小標本に対しては成立しない．それらの多くは，母集団分布が正規で √ b/ n となる．ここで， データが対称であるとき，一つの平均に適用したとき，近似的に x ± zα/2 σ σ b2 = n−1 ∑ (xi − x)2 である．これらは非常に短く，t 信頼区間より短い．それは 2 つの原因がある． それは，tα/2,n−1 の代わりに zα/2 を，s の代わりに σ b を用いているからである．これらの項は O(1/n) 9 であり，漸近 2 次の正確性に対して影響を与えないが，n が小さいとき，問題となる．シミュレーショ ンによる比較とこれらの項の修正に関しては，Hesterberg (1999) を参照． 歪度が問題とならない状況についても述べておく．例えば，2 つの平均の比較において，サンプルの 大きさが等しく，母集団の標準偏差 s および歪度が等しいとき，歪度はキャンセルされ，t 分布は対称 となる． そのような構造的な解消が生じない状況を考える．サンプルサイズが大きい（例えば，n > 100）が， 極端に大きくないとき（例えば，n < 5000），あるいは，歪度が大きくなったり，本質的に歪んだ統計 量でない限りにおいて，2 次のオーダーで正確な手法を用いるべきである． 中・小標本に対しては，おそらく 2 次の正確さを持つ手法を利用するべきであるが，このとき注意が 必要であり，それらを短くしすぎる因子に対する調整，おそらく歪度を 0 に向けて調整する縮小であ る．こうした推奨方法を定量化するには，さらなる研究が必要である． 最後の理由は，慣性である．人は簡単な診断を利用するものである．だから，統計に関する教育が必 要となる，私たちは，古いルールではなく，ブートストラップ診断を教えるべきである．これには別の 利点がある．ブートストラッピングを行い，ブートストラップ分布の図を見ることは，標本分布や中心 極限定理，標準誤差，バイアス，P 値といった統計概念をよりよく理解する助けとなる． 参考文献 [1] Brashares, J. S., Arcese, P., Sam, M. K., Coppolillo, P. B., Sinclair, A. R. E. and Balmford, A. (2004) Bushmeat hunting, wildlife declines, and fish supply in west africa. Science, 306, 1180-1183. [2] Chambers, J. and Hastie, T. (1992) Statistical Models in S. Wadsworth, California. [3] Davison, A. and Hinkley, D. (1997) Bootstrap Methods and their Applications. Cambridge University Press. [4] DiCiccio, T. and Efron, B. (1992) More accurate confidence intervals in exponential families. Biometrika, 79, 231-245. [5] Efron, B. (1981) Nonparametric standard errors and confidence intervals. Canadian Journal of Statistics, 9, 139 - 172. [6] Efron, B. (1982) The Jackknife, the Bootstrap and Other Resampling Plans. National Science Foundation - Conference Board of the Mathematical Sciences Monograph 38. Philadelphia: Society for Industrial and Applied Mathematics. [7] Efron, B. (1987) Better bootstrap confidence intervals (with discussion). Journal of the American Statistical Association, 82, 171 - 200. [8] Efron, B. and Tibshirani, R. J. (1993) An Introduction to the Bootstrap. Chapman and Hall. [9] Hesterberg, T., Monaghan, S., Moore, D. S., Clipson, A. and Epstein, R. (2003) Bootstrap Methods and Permutation Tests. W. H. Freeman. URL http://bcs.whfreeman.com/ ips5e/content/cat_080/pdf/moore14.pdf. Chapter for The Practice of Business Statistics by Moore, McCabe, Duckworth, and Sclove. [10] Hesterberg, T. C. (1999) Bootstrap tilting confidence intervals. Research Department Technical Report 84, MathSoft, Inc. URL http://www.insightful.com/Hesterberg/articles/ tech84-tiltingCI.pdf. [11] Johnson, N. J. (1978) Modified t tests and confidence intervals for asymmetrical populations. Journal of the American Statistical Association, 73, 536-544 10