標本サイズの選択
この例では、統計的検定の実行に必要な標本や観測の件数を決定する方法を示します。ここでは、単純な問題に対する標本サイズの計算を示した後、関数 sampsizepwr
を使用して、2 つのより現実的な問題に対して検出力と標本サイズを計算する方法を示します。最後に、Statistics and Machine Learning Toolbox™ の関数を使用して、関数 sampsizepwr
がサポートしていない検定に必要な標本サイズを計算する方法を示します。
既知の標準偏差を使用した正規平均の片側検定
概念を紹介するために、非現実的な単純な例を使用して、平均を検定し、標準偏差を確認してみましょう。データは連続しており、 正規分布を使用してモデル化できます。ここでは標本サイズ N を判断し、100 という平均と 110 という平均を区別できるようにする必要があります。標準偏差は 20 であることがわかっています。
統計的検定を実施するとき、通常は対立仮説に対して帰無仮説を検定します。検定統計量 T を求め、帰無仮説のもとでその分布を確認します。たとえば帰無仮説が真の場合の発生確率が 5% 未満であるなどの異常な値が見つかった場合は、対立仮説を優先して帰無仮説を棄却します (5% の確率は、検定の有意水準と呼ばれます)。値が異常でない場合、帰無仮説は棄却しません。
この場合、検定統計量 T はサンプル平均です。帰無仮説のもとでは、平均は 100 で、標準偏差は 20/sqrt(N) です。まず、N=16 の固定標本サイズを見てみましょう。Tが影の領域、つまり分布の上裾にある場合は、帰無仮説を棄却します。T が下裾にある場合には棄却しないため、これは片側検定になります。この影の領域のカットオフは、平均より上の 1.6 標準偏差になります。
rng(0,'twister'); mu0 = 100; sig = 20; N = 16; alpha = 0.05; conf = 1-alpha; cutoff = norminv(conf, mu0, sig/sqrt(N)); x = [linspace(90,cutoff), linspace(cutoff,127)]; y = normpdf(x,mu0,sig/sqrt(N)); h1 = plot(x,y); xhi = [cutoff, x(x>=cutoff)]; yhi = [0, y(x>=cutoff)]; patch(xhi,yhi,'b'); title('Distribution of sample mean, N=16'); xlabel('Sample mean'); ylabel('Density'); text(96,.01,sprintf('Reject if mean>%.4g\nProb = 0.05',cutoff),'Color','b');
これは、帰無仮説のものでの T の動作を示していますが、対立仮説のもとではどうでしょうか。対立分布の平均は、赤の曲線で示されるように 110 です。
mu1 = 110; y2 = normpdf(x,mu1,sig/sqrt(N)); h2 = line(x,y2,'Color','r'); yhi = [0, y2(x>=cutoff)]; patch(xhi,yhi,'r','FaceAlpha',0.25); P = 1 - normcdf(cutoff,mu1,sig/sqrt(N)); text(115,.06,sprintf('Reject if T>%.4g\nProb = %.2g',cutoff,P),'Color',[1 0 0]); legend([h1 h2],'Null hypothesis','Alternative hypothesis');
対立仮説が真の場合に帰無仮説を棄却する確率は高くなっています。これは、望むとおりの現象です。確率密度関数 (pdf) ではなく累積分布関数 (cdf) を見てみると、さらにわかりやすくなります。確率は、面積を計算しなくても、グラフから直接読み取ることができます。
ynull = normcdf(x,mu0,sig/sqrt(N)); yalt = normcdf(x,mu1,sig/sqrt(N)); h12 = plot(x,ynull,'b-',x,yalt,'r-'); zval = norminv(conf); cutoff = mu0 + zval * sig / sqrt(N); line([90,cutoff,cutoff],[conf, conf, 0],'LineStyle',':'); msg = sprintf(' Cutoff = 100 + %.2g \\times 20 / \\surd{n}',zval); text(cutoff,.15,msg,'Color','b'); text(min(x),conf,sprintf(' %g%% test',100*alpha),'Color','b',... 'verticalalignment','top') palt = normcdf(cutoff,mu1,sig/sqrt(N)); line([90,cutoff],[palt,palt],'Color','r','LineStyle',':'); text(91,palt+.02,sprintf(' Power is 1-%.2g = %.2g',palt,1-palt),'Color',[1 0 0]); legend(h12,'Null hypothesis','Alternative hypothesis');
このグラフは、N=16 である場合に 2 つの異なる mu 値に対して有意な統計値を得る (帰無仮説を棄却する) 確率を示しています。
検出力関数は、対立仮説が真の場合に帰無仮説を棄却する確率として定義されています。これは、対立仮説の値と標本サイズに依存します。対立仮説値を 110 に固定し、N の関数として検出力 (1 から cdf を減算した値) をグラフにします。検出力が 80% になるように N を選択します。グラフは、約 N=25 が必要であることを示しています。
DesiredPower = 0.80; Nvec = 1:30; cutoff = mu0 + norminv(conf)*sig./sqrt(Nvec); power = 1 - normcdf(cutoff, mu1, sig./sqrt(Nvec)); plot(Nvec,power,'bo-',[0 30],[DesiredPower DesiredPower],'k:'); xlabel('N = sample size') ylabel('Power') title('Power function for the alternative: \mu = 110')
この非常に単純な例では、必要な値を直接計算して 80% の検出力を得るための公式を使用します。
mudiff = (mu1 - mu0) / sig; N80 = ceil(((norminv(1-DesiredPower)-norminv(conf)) / mudiff)^2)
N80 = 25
これが機能することを確認するため、モンテカルロ シミュレーションを行い、平均が 100 であるという帰無仮説と平均が 110 であるという対立仮説に基づいて、サイズが 25 の標本を 400 個生成してみましょう。平均が 100 であるかどうかを調べるために各標本を検定する場合は、1 番目のグループの約 5% と 2 番目のグループの約 80% が有意になる必要があります。
nsamples = 400; samplenum = 1:nsamples; N = 25; h0 = zeros(1,nsamples); h1 = h0; for j = 1:nsamples Z0 = normrnd(mu0,sig,N,1); h0(j) = ztest(Z0,mu0,sig,alpha,'right'); Z1 = normrnd(mu1,sig,N,1); h1(j) = ztest(Z1,mu0,sig,alpha,'right'); end p0 = cumsum(h0) ./ samplenum; p1 = cumsum(h1) ./ samplenum; plot(samplenum,p0,'b-',samplenum,p1,'r-') xlabel('Sample number') ylabel('Proportion significant') title('Verification of power computation') legend('Null hypothesis','Alternative hypothesis','Location','East')
不明な標準偏差を使用した正規平均の両側検定
次に、標準偏差がわからない場合に、両側検定、つまり、サンプル平均が高すぎるか低すぎるかによって帰無仮説を棄却する検定を実施するとします。
検定の統計値は t 統計値、つまり、サンプル平均と検定対象の平均の差を平均の標準誤差で割ったものです。帰無仮説のもとでは、これは N-1 の自由度をもつスチューデントの t 分布を含みます。対立仮説のもとでは、分布は、真の平均と検定対象の平均の正規化された差に等しい非心度パラメーターをもつ非心 t 分布になります。
この両側検定に対しては、帰無仮説のもとで、両方の裾に等しい 5% の誤差確率を割り当て、検定統計量がどちらかの方向において極端すぎる場合に棄却する必要があります。また、どの対立仮説のもとでも両方の裾を考慮する必要があります。
N = 16; df = N-1; alpha = 0.05; conf = 1-alpha; cutoff1 = tinv(alpha/2,df); cutoff2 = tinv(1-alpha/2,df); x = [linspace(-5,cutoff1), linspace(cutoff1,cutoff2),linspace(cutoff2,5)]; y = tpdf(x,df); h1 = plot(x,y); xlo = [x(x<=cutoff1),cutoff1]; ylo = [y(x<=cutoff1),0]; xhi = [cutoff2,x(x>=cutoff2)]; yhi = [0, y(x>=cutoff2)]; patch(xlo,ylo,'b'); patch(xhi,yhi,'b'); title('Distribution of t statistic, N=16'); xlabel('t'); ylabel('Density'); text(2.5,.05,sprintf('Reject if t>%.4g\nProb = 0.025',cutoff2),'Color','b'); text(-4.5,.05,sprintf('Reject if t<%.4g\nProb = 0.025',cutoff1),'Color','b');
帰無仮説のみのもとでの検出力関数と、mu の単一の対立仮説値を調べる代わりに、それを mu の関数として見ることができます。検出力は、mu がどちらかの方向において帰無仮説から遠ざかるにつれて増加します。検出力を計算するには、関数 sampsizepwr
を使用できます。検出力の計算には、標準偏差の値を指定する必要があります。これは約 20 になることが予想されます。標本サイズ N=16 の検出力関数の図は次のようになります。
N = 16; x = linspace(90,127); power = sampsizepwr('t',[100 20],x,[],N); plot(x,power); xlabel('True mean') ylabel('Power') title('Power function for N=16')
平均が 110 のときに検出力が 80% になる必要があります。このグラフによると、N=16 という標本サイズでは検出力が 50% 未満になります。どの標本サイズを使用すれば目的の検出力が得られるのでしょうか。
N = sampsizepwr('t',[100 20],110,0.8)
N = 34
約 34 という標本サイズが必要です。前の例と比較すると、不明な真の標準偏差を補って両側検定を行うには、観測値を 9 個追加する必要があります。
検出力関数のプロットは、さまざまな N の値に対して作成できます。
Nvec = 2:40; power = sampsizepwr('t',[100 20],110,[],Nvec); plot(Nvec,power,'bo-',[0 40],[DesiredPower DesiredPower],'k:'); xlabel('N = sample size') ylabel('Power') title('Power function for the alternative: \mu = 110')
そして、前回と同じようなシミュレーションを実行し、必要な検出力が得られるかどうかを確認できます。
nsamples = 400; samplenum = 1:nsamples; N = 34; h0 = zeros(1,nsamples); h1 = h0; for j = 1:nsamples Z0 = normrnd(mu0,sig,N,1); h0(j) = ttest(Z0,mu0,alpha); Z1 = normrnd(mu1,sig,N,1); h1(j) = ttest(Z1,mu0,alpha); end p0 = cumsum(h0) ./ samplenum; p1 = cumsum(h1) ./ samplenum; plot(samplenum,p0,'b-',samplenum,p1,'r-') xlabel('Sample number') ylabel('Proportion significant') title('Verification of power computation') legend('Null hypothesis','Alternative hypothesis','Location','East')
50 という標本サイズを用意する余裕があるとします。対立仮説値 mu=110 を検出するための検出力は 80% より大きくなることが予想されます。80% の検出力を維持する場合は、どのような対立仮説を検出できるでしょうか。
mu1 = sampsizepwr('t',[100 20],[],.8,50)
mu1 = 108.0837
比率の検定
次に、2 つの比率を区別するために必要な標本サイズを判断する問題を考えましょう。約 30% の人々が候補者を支持している人口をサンプリングする場合において、この値を 33% とは区別できるだけの十分な数の人々をサンプリングするとします。
ここでの観念は、前回と同じです。ここでは、サンプル数を検定統計量として使用できます。この数には二項分布が含まれています。任意の標本サイズ N に対して、帰無仮説 P=0.30 を棄却するためのカットオフを計算できます。たとえば N=100 の場合は、サンプル数が以下のように計算されたカットオフ値よりも大きい場合に帰無仮説を棄却します。
N = 100; alpha = 0.05; p0 = 0.30; p1 = 0.33; cutoff = binoinv(1-alpha, N, p0)
cutoff = 38
二項分布は離散分布であるため、複雑さが伴います。カットオフ値を超過する確率は、正確に 5% ではありません。
1 - binocdf(cutoff, N, p0)
ans = 0.0340
ここでも、標本サイズの範囲における対立仮説値 P=0.33 に対する検出力を計算しましょう。ここでは 30% より大きい対立仮説値のみに注目する必要があるため、片側 (右裾) 検定を使用します。
Nvec = 50:50:2000; power = sampsizepwr('p',p0,p1,[],Nvec,'tail','right'); plot(Nvec,power,'bo-',[0 2000],[DesiredPower DesiredPower],'k:'); xlabel('N = sample size') ylabel('Power') title('Power function for the alternative: p = 0.33')
関数 sampsizepwr
を使用して、80% の検出力に必要な標本サイズを要求します。
approxN = sampsizepwr('p',p0,p1,0.80,[],'tail','right')
Warning: Values N>200 are approximate. Plotting the power as a function of N may reveal lower N values that have the required power. approxN = 1500
警告メッセージにより、回答は近似であることが伝えられます。異なる標本サイズに対して検出力関数を見てみると、二項分布は離散であるため、関数は一般的に増加しますが、不規則であることがわかります。1470 ~ 1480 の標本サイズ範囲において、p=0.30 および p=0.33 の両方に対して帰無仮説を棄却する確率を見てみましょう。
subplot(3,1,1); Nvec = 1470:1480; power = sampsizepwr('p',p0,p1,[],Nvec,'tail','right'); plot(Nvec,power,'ro-',[min(Nvec),max(Nvec)],[DesiredPower DesiredPower],'k:'); ylabel(sprintf('Prob[T>cutoff]\nif p=0.33')) h_gca = gca; h_gca.XTickLabel = ''; ylim([.78, .82]); subplot(3,1,2); alf = sampsizepwr('p',p0,p0,[],Nvec,'tail','right'); plot(Nvec,alf,'bo-',[min(Nvec),max(Nvec)],[alpha alpha],'k:'); ylabel(sprintf('Prob[T>cutoff]\nif p=0.30')) h_gca = gca; h_gca.XTickLabel = ''; ylim([.04, .06]); subplot(3,1,3); cutoff = binoinv(1-alpha, Nvec, p0); plot(Nvec,cutoff,'go-'); xlabel('N = sample size') ylabel('Cutoff')
このプロットは、検出力関数の曲線 (上部のプロット) が不規則なだけでなく、一部の標本サイズで減少していることを示しています。これらの標本サイズでは、5% 以下の有意水準 (中央のプロット) を保持するためにカットオフ値 (下部のプロット) を増加する必要があります。この範囲内で、目的の 80% の検出力を提供するより小さい標本サイズを見つけることができます。
min(Nvec(power>=0.80))
ans = 1478
相関の検定
これまでに検討した例では、特定の有意水準を達成するための検定統計量のカットオフを求め、対立仮説のもとでそのカットオフを超える確率を計算することができました。最後の例では、少し難しい問題を考えてみましょう。
2 つの変数 X および Y からのサンプルを使用する場合に、これらが相関していないか、あるいは相関が最高で 0.4 にも及ぶことがあるかを検査するために必要な標本サイズを知りたいとします。サンプルの相関の分布を t 分布に変換することで解を得ることは可能ですが、検定統計量の分布に対処できない問題においても使用できる方法を採用してみましょう。
任意の標本サイズに対してモンテ カルロ シミュレーションを使用すると、相関の検定に適したカットオフ値を判断できます。大きなシミュレーションを実行してこの値を正確に求めてみましょう。まず、25 の標本サイズから始めます。
nsamples = 10000; N = 25; alpha = 0.05; conf = 1-alpha; r = zeros(1,nsamples); for j = 1:nsamples xy = normrnd(0,1,N,2); r(j) = corr(xy(:,1),xy(:,2)); end cutoff = quantile(r,conf)
cutoff = 0.3372
次に、対立仮説に基づいてサンプルを生成し、検定の検出力を推定します。
nsamples = 1000; mu = [0; 0]; sig = [1 0.4; 0.4 1]; r = zeros(1,nsamples); for j = 1:nsamples xy = mvnrnd(mu,sig,N); r(j) = corr(xy(:,1),xy(:,2)); end [power,powerci] = binofit(sum(r>cutoff),nsamples)
power = 0.6470 powerci = 0.6165 0.6767
検出力は 65% であると推定し、95% の信頼をもって真の値が 62% ~ 68% の範囲内にあるとします。80% の検出力を得るには、より大きな標本サイズが必要です。N を 50 に増加し、この標本サイズのカットオフ値を推定して、検出力のシミュレーションを繰り返すことができます。
nsamples = 10000; N = 50; alpha = 0.05; conf = 1-alpha; r = zeros(1,nsamples); for j = 1:nsamples xy = normrnd(0,1,N,2); r(j) = corr(xy(:,1),xy(:,2)); end cutoff = quantile(r,conf) nsamples = 1000; mu = [0; 0]; sig = [1 0.4; 0.4 1]; r = zeros(1,nsamples); for j = 1:nsamples xy = mvnrnd(mu,sig,N); r(j) = corr(xy(:,1),xy(:,2)); end [power,powerci] = binofit(sum(r>cutoff),nsamples)
cutoff = 0.2315 power = 0.8990 powerci = 0.8786 0.9170
この標本サイズでは、80% の目標を上回る検出力が得られます。このようにして試行錯誤を繰り返すと、ここでの要件を満たす 50 未満の標本サイズを見つけることができます。
まとめ
Statistics and Machine Learning Toolbox の確率関数を使用すると、仮説検定で要求される検出力のレベルを実現するために必要な標本サイズを判断できます。問題によっては、標本サイズを直接計算できる場合と、正しい値が見つかるまで標本サイズの範囲を探し求める必要が生じる場合があります。乱数発生器は、目的の検出力が満たされているかどうかを確認するのに役立ちます。また、対立条件下で特定の検定のパワーを調べるためにも使用できます。