非線形モデルにおける検出力の事前計算
この例では、6つの連続尺度の因子が、部品検査の合否におよぼす主効果を調べたいとしましょう。応答は二項分布に従い、実験は60回行えるものとします。
この例は、次の手順で進めます。
1. 実験のカスタム計画を作成します。計画の作成を参照してください。
注: カスタム計画は、非線形な状況に最適とはいえませんが、この例では、単純化するために、「非線形計画」プラットフォームではなく「カスタム計画」プラットフォームを使用します。「非線形計画」プラットフォームで作成した計画のほうが直交計画より優れていることを示す例については、『実験計画(DOE)』の非線形計画の例を参照してください。
2. 一般化線形モデルを使って、ロジスティックモデルをあてはめます。一般化線形モデルのあてはめを参照してください。
3. 尤度比検定のp値をシミュレートして、線形予測子によって異なる確率の差の検出力を調べます。検出力の考察を参照してください。
この例での計画
一般化線形モデルで、リンク関数としてロジット関数を使用して、成功確率をモデル化します。リンク関数としてロジット関数を使用すると、次のロジスティックモデルがあてはめられます。
ここで、p(X)は、X = (X1, X2, ..., X6)という因子設定のもと、部品が合格する確率を表します。
線形予測子L(X)は、次の式で表されます。
線形予測子の係数が次のような場合だったときの検出力を求めてみましょう。
|
係数 |
値 |
|---|---|
|
b0 |
0 |
|
b1 |
1 |
|
b2 |
0.9 |
|
b3 |
0.8 |
|
b4 |
0.7 |
|
b5 |
0.6 |
|
b6 |
0.5 |
線形予測子の切片は0なので、全因子が0に設定されている場合、部品が合格する確率は50%に等しくなります。i番目の因子における各水準に関連する確率は、その他の因子がすべて0である場合、次のようになります。
|
因子 |
Xi = 1の場合の合格率 |
Xi = -1の場合の合格率 |
差 |
|---|---|---|---|
|
X1 |
73.11% |
26.89% |
46.2% |
|
X2 |
71.09% |
28.91% |
42.2% |
|
X3 |
69.00% |
31.00% |
38.0% |
|
X4 |
66.82% |
33.18% |
33.6% |
|
X5 |
64.56% |
35.43% |
29.1% |
|
X6 |
62.25% |
37.75% |
24.5% |
たとえば、X1以外の全因子が0に設定されている場、検出したい差は46.2%になります。検出したい差が最小になるのは、X6以外の全因子が0に設定された場合で、その差は24.5%です。
計画の作成
注: この節の手順を省略するには、 ［ヘルプ］>［サンプルデータライブラリ］を選択し、「Design Experiment」フォルダの「Binomial Experiment.jmp」を開きます。「DOE シミュレート」 スクリプトの横にある緑色の三角形をクリックし、［応答をシミュレート］の設定に進んでください。
1. ［実験計画(DOE)］>［カスタム計画］を選択します。
2. 「因子」アウトラインで、「N個の因子を追加」の横のボックスに「6」と入力します。
3. ［因子の追加］>［連続変数］を選択します。
4. ［続行］をクリックします。
主効果の計画を作成したいため、「モデル」アウトラインには変更を加えません。
5. 「実験の回数」の「ユーザ定義」テキストボックスに「60」と入力します。
6. 「カスタム計画」の赤い三角ボタンをクリックし、メニューから［応答のシミュレート］を選択します。
この操作を行っておくと、計画テーブルを作成するために［テーブルの作成］をクリックしたときに、「応答をシミュレート」ウィンドウが開きます。
注: 乱数シード値（ステップ7）と開始点の数（ステップ8）を設定すると、以下の数値例と同じ実験設定が得られます。同じ実験設定でなくても良い場合は、これらの手順は不要です。
7. （オプション）「カスタム計画」の赤い三角ボタンをクリックして、［乱数シード値の設定］を選択します。「12345」と入力して［OK］をクリックします。
8. （オプション）「カスタム計画」の赤い三角ボタンをクリックして、［開始点の数］を選択します。「1」と入力して［OK］をクリックします。
9. ［計画の作成］をクリックします。
10. ［テーブルの作成］をクリックします。
注: 乱数を用いているため、「Y」および「Yのシミュレーション」の値は、図とは異なったものになります。
計画テーブル（一部）
「応答をシミュレート」ウィンドウ
計画テーブルと「応答をシミュレート」ウィンドウが表示されます。計画テーブルに次の2つの列が追加されます。
– 「Y」列：「応答をシミュレート」ウィンドウの設定に従ってシミュレートされたデータ値。
– 「Yのシミュレーション」列：計算式が含まれており、その計算結果が表示されている。「応答をシミュレート」ウインドウでモデルを設定すると、そのモデルに従った乱数が生成されます。計算式を表示するには、「列」パネルで列名の右側にある＋記号をクリックしてください。
次の節では、二項分布に従う応答をシミュレートして、一般化線形モデルをあてはめます。
［応答をシミュレート］の設定
ここでは、二項分布に従う応答データをシミュレートします。成功率はロジスティックモデルに従うものとします。［応答をシミュレート］の詳細については、『実験計画(DOE)』の応答のシミュレートを参照してください。
注: この節の手順を省略するには、「二項分布のシミュレート」スクリプトの横にある緑の三角ボタンをクリックします。この操作が済んだら、一般化線形モデルのあてはめに進んでください。
1. 「応答をシミュレート」ウィンドウ（図）で、「Y」の下に以下の値を入力します。
– 「切片」に「0」と入力します。
– 「X1」は、デフォルト値の「1」をそのまま使います。
– 「X2」に「0.9」と入力します。
– 「X3」に「0.8」と入力します。
– 「X4」に「0.7」と入力します。
– 「X5」に「0.6」と入力します。
– 「X6」に「0.5」と入力します。
2. 「分布」アウトラインから［二項］を選択します。
「N」の値は「1」のままにしておきます。これは、それぞれの試行回数が1回だけであることを示します。
設定後の「応答をシミュレート」ウィンドウ
3. ［適用］をクリックします。
データテーブルの「Yのシミュレーション」列が、二項分布に従う乱数を生成する計算式に置き換えられます。「Y試行回数」の列には、各実験における試行回数が示されます。
4. （オプション）「列」パネルで、「Yのシミュレーション」の右側にある＋記号をクリックします。
「Yのシミュレーション」に挿入される、二項分布に従う乱数の計算式
5. ［キャンセル］をクリックします。
一般化線形モデルのあてはめ
1. データテーブルの「モデル」スクリプトの横の緑の三角ボタンをクリックします。
2. ［Y］ボタンの横の変数「Y」をクリックし、［削除］ボタンをクリックします。
3. 「Yのシミュレーション」をクリックし、［Y］ボタンをクリックします。
応答変数Yとして、二項分布に従う乱数の列が設定されます。
4. 「手法」リストから［一般化線形モデル］を選択します。
5. 「分布」リストから［二項］を選択します。
［リンク関数］メニューに［ロジット］が表示されます。
6. ［実行］をクリックします。
二項分布のデータ1組に対して、モデルがあてはめられます。
検出力の考察
この例での計画では、線形予測子および係数と、係数によって決められる確率を示しました。ここでは、それらの係数に対する検出力を求めます。
1. 「効果の検定」アウトラインで、「p値(Prob>ChiSq)」列を右クリックして［シミュレーション］を選択します。
「シミュレーション」ウィンドウ
「切り替え元の列」リストで［Yのシミュレーション］列が選択されていることを確認します。この列は、モデルのあてはめに使用されたデータ値を含んでいます。「切り替え先の列」でも「Yのシミュレーション」を選択すると、切り替え元の「Yのシミュレーション」データ値を、「Yのシミュレーション」の計算式で指定されている乱数に置き換えて、シミュレーションが実行されます。
レポートで選択した「p値(Prob>ChiSq)」列は、「係数がゼロである」という帰無仮説に対する尤度比検定のp値です。「効果の検定」表にある、各効果に対するこの「p値(Prob>ChiSq)」をシミュレートしましょう。
2. 「標本数」に「500」と入力します。
3. ［OK］をクリックします。
「一般化線形モデル シミュレーション結果」データテーブルが表示されます。
注: シミュレーションは乱数によって行われるため、結果は図とは異なる可能性があります。
「シミュレーション結果」テーブル（一部）
テーブルの最初の行は、観測されたデータから得られた「p値(Prob>ChiSq)」であり、除外の状態になっています。2行目以降の500行は、シミュレーションしたデータから計算された結果です。
4. 「検出力の分析」スクリプトを実行します。
注: 応答値はシミュレートされた乱数であるため、検出力のシミュレーションの結果は、図とは異なる可能性があります。
一変量の分布（最初の3つの効果）
各主効果のヒストグラムには、シミュレーションによって得られた500個の「p値(Prob>ChiSq)」が描かれています。また、「検出力のシミュレーション」アウトラインには、500回のシミュレーションのうち有意性検定で棄却されたものの割合が表示されます。
次に、見やすくするために、レポートを縦に積み重ねて、プロットを非表示にしてみましょう。
5. 「一変量の分布」の赤い三角ボタンをクリックし、［積み重ねて表示］を選択します。
6. Ctrlキーを押しながら、「X1」の赤い三角ボタンをクリックし、［外れ値の箱ひげ図］の選択を解除します。
7. Ctrlキーを押しながら「X1」の赤い三角ボタンをクリックし、［ヒストグラムオプション］を選択し、［ヒストグラム］の選択を解除します。
注: 応答値はシミュレートされた乱数であるため、検出力のシミュレーションの結果は、図とは異なる可能性があります。
検出力（最初の3つの効果）
「検出力のシミュレーション」アウトラインの「棄却された割合」は、シミュレーションのうち、p値がαより小さくなったものの割合を示します。たとえば、係数が0.8、確率の差が38%である「X3」で、棄却された割合を見てみると、有意水準0.05の検定で379/500 = 0.758となっています。表は、シミュレーションで計算された、有意水準0.05の検定に対する検出力を、効果ごとに示したものです。「検出したい差」が小さくなるにつれ、検出力も小さくなっています。たとえば、「検出したい差」が24.5%（X6）における検出力は、約0.37にとどまっています。
注: 応答値はシミュレートされた乱数であるため、検出力のシミュレーションの結果は、表とは異なる可能性があります。
因子 | Xi = 1の場合の | Xi = -1の場合の | 検出したい差 | シミュレーションで得られた検出力（α=0.05で棄却された割合） |
|---|---|---|---|---|
X1 | 73.11% | 26.89% | 46.2% | 0.852 |
X2 | 71.09% | 28.91% | 42.2% | 0.828 |
X3 | 69.00% | 31.00% | 38.0% | 0.758 |
X4 | 66.82% | 33.18% | 33.6% | 0.654 |
X5 | 64.56% | 35.43% | 29.1% | 0.488 |
X6 | 62.25% | 37.75% | 24.5% | 0.372 |