非線形モデルにおける検出力の事前計算

この例では、6つの連続尺度の因子が、部品検査の合否におよぼす主効果を調べたいとしましょう。応答は二項分布に従い、実験は60回行えるものとします。

メモ: カスタム計画は、非線形な状況に最適とはいえませんが、この例では、単純化するために、「非線形計画」プラットフォームではなく「カスタム計画」プラットフォームを使用します。「非線形計画」プラットフォームで作成した計画のほうが直交計画より優れていることを示す例については、『実験計画(DOE)』の「非線形計画」章を参照してください。

一般化線形モデルを使って、ロジスティックモデルをあてはめます。第 “一般化線形モデルのあてはめ”を参照してください。

尤度比検定のp値をシミュレートして、線形予測子によって異なる確率の差の検出力を調べます。第 “検出力の考察”を参照してください。

この例での計画

一般化線形モデルで、リンク関数としてロジット関数を使用して、成功確率をモデル化します。リンク関数としてロジット関数を使用すると、次のロジスティックモデルがあてはめられます。

ここで、π(X)は、X = (X1, X2, ..., X6)という因子設定のもと、部品が合格する確率を表します。

線形予測子L(X)は、次の式で表されます。

線形予測子の係数が次のような場合だったときの検出力を求めてみましょう。

係数	値
β0	0
β1	1
β2	.9
β3	.8
β4	.7
β5	.6
β6	.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 シミュレート」スクリプトの横にある緑色の三角形をクリックし、第 “［応答をシミュレート］の設定”に進んでください。

［実験計画(DOE)］>［カスタム計画］を選択します。

「因子」アウトラインで、「N個の因子を追加」の横のボックスに「6」と入力します。

［因子の追加］>［連続変数］を選択します。

［続行］をクリックします。

主効果の計画を作成したいため、「モデル」アウトラインには変更を加えません。

「実験の回数」の「ユーザ定義」テキストボックスに「60」と入力します。

「カスタム計画」の赤い三角ボタンをクリックし、メニューから［応答のシミュレート］を選択します。

この操作を行っておくと、計画テーブルを作成するために［テーブルの作成］をクリックしたときに、「応答をシミュレート」ウィンドウが開きます。

メモ: 乱数シード値（第 7 步）と開始点の数（第 8 步）を設定すると、以下の数値例と同じ実験設定が得られます。同じ実験設定でなくても良い場合は、これらの手順は不要です。

（オプション）「カスタム計画」の赤い三角ボタンをクリックして、［乱数シード値の設定］を選択します。「12345」と入力して［OK］をクリックします。

（オプション）「カスタム計画」の赤い三角ボタンをクリックして、［開始点の数］を選択します。「1」と入力して［OK］をクリックします。

［計画の作成］をクリックします。

10.

［テーブルの作成］をクリックします。

メモ: 乱数を用いているため、「Y」および「Yのシミュレーション」の値は、計画テーブル（一部）とは異なったものになります。

図10.17 計画テーブル（一部）

図10.18 「応答をシミュレート」ウィンドウ

計画テーブルと「応答をシミュレート」ウィンドウが表示されます。計画テーブルに次の2つの列が追加されます。

–

「Y」列：「応答をシミュレート」ウィンドウの設定に従ってシミュレートされたデータ値。

–

「Yのシミュレーション」列：計算式が含まれており、その計算結果が表示されている。「応答をシミュレート」ウインドウでモデルを設定すると、そのモデルに従った乱数が生成されます。計算式を表示するには、「列」パネルで列名の右側にある＋記号をクリックしてください。

次の節では、二項分布に従う応答をシミュレートして、一般化線形モデルをあてはめます。

［応答をシミュレート］の設定

ここでは、二項分布に従う応答データをシミュレートします。成功率はロジスティックモデルに従うものとします。［応答をシミュレート］の詳細については、『実験計画(DOE)』の「カスタム計画」章を参照してください。

メモ: この節の手順を省略するには、「二項分布のシミュレート」スクリプトの横にある緑の三角ボタンをクリックします。この操作が済んだら、第 “一般化線形モデルのあてはめ”に進んでください。

「応答をシミュレート」ウィンドウ（「応答をシミュレート」ウィンドウ）で、「Y」の下に以下の値を入力します。

–

「切片」に「0」と入力します。

–

「X1」は、デフォルト値の「1」をそのまま使います。

–

「X2」に「0.9」と入力します。

–

「X3」に「0.8」と入力します。

–

「X4」に「0.7」と入力します。

–

「X5」に「0.6」と入力します。

–

「X6」に「0.5」と入力します。

「分布」アウトラインから［二項］を選択します。

「N」の値は「1」のままにしておきます。これは、それぞれの試行回数が1回だけであることを示します。

図10.19 設定後の「応答をシミュレート」ウィンドウ

［適用］をクリックします。

データテーブルの「Yのシミュレーション」列が、二項分布に従う乱数を生成する計算式に置き換えられます。「Y試行回数」の列には、各実験における試行回数が示されます。

（オプション）「列」パネルで、「Yのシミュレーション」の右側にある＋記号をクリックします。

図10.20 「Yのシミュレーション」に挿入される、二項分布に従う乱数の計算式

［キャンセル］をクリックします。

一般化線形モデルのあてはめ

データテーブルの「モデル」スクリプトの横の緑の三角ボタンをクリックします。

［Y］ボタンの横の変数「Y」をクリックし、［削除］ボタンをクリックします。

「Yのシミュレーション」をクリックし、［Y］ボタンをクリックします。

応答変数Yとして、二項分布に従う乱数の列が設定されます。

「手法」のメニューから［一般化線形モデル］を選択します。

［分布］メニューから［二項］を選択します。

［リンク関数］メニューに［ロジット］が表示されます。

［実行］をクリックします。

二項分布のデータ1組に対して、モデルがあてはめられます。

検出力の考察

第 “この例での計画”では、線形予測子および係数と、係数によって決められる確率を示しました。ここでは、それらの係数に対する検出力を求めます。

「効果の検定」アウトラインで、「p値(Prob>ChiSq)」列を右クリックして［シミュレーション］を選択します。

図10.21 「シミュレーション」ウィンドウ

「切り替え元の列」の「Yのシミュレーション」は、モデルの推定に使用されたデータ値を含んでいます。「切り替え先の列」でも「Yのシミュレーション」を選択すると、切り替え元の「Yのシミュレーション」データ値を、「Yのシミュレーション」の計算式で指定されている乱数に置き換えて、シミュレーションが実行されます。

レポートで選択した「p値(Prob>ChiSq)」列は、「係数がゼロである」という帰無仮説に対する尤度比検定のp値です。「効果の検定」表にある、各効果に対するこの「p値(Prob>ChiSq)」をシミュレートしましょう。

「標本数」に「500」と入力します。

（オプション）「乱数シード値」に「123」と入力し、テキストボックスの外をクリックします。

［OK］をクリックします。

「一般化線形モデルシミュレーション結果」データテーブルが表示されます。

メモ: シミュレーションは乱数によって行われるため、結果は「シミュレーション結果」テーブル（一部）とは異なります。

図10.22 「シミュレーション結果」テーブル（一部）

テーブルの最初の行は、観測されたデータから得られた「p値(Prob>ChiSq)」であり、除外の状態になっています。2行目以降の500行は、シミュレーションしたデータから計算された結果です。

「検出力の分析」スクリプトを実行します。

メモ: 応答値はシミュレートされた乱数であるため、検出力のシミュレーションの結果は、一変量の分布（最初の3つの効果）とは異なります。

図10.23 一変量の分布（最初の3つの効果）

各主効果のヒストグラムには、シミュレーションによって得られた500個の「p値(Prob>ChiSq)」が描かれています。また、「検出力のシミュレーション」アウトラインには、500回のシミュレーションのうち有意性検定で棄却されたものの割合が表示されます。

次に、見やすくするために、レポートを縦に積み重ねて、プロットを非表示にしてみましょう。

「一変量の分布」の赤い三角ボタンをクリックし、［積み重ねて表示］を選択します。

Ctrlキーを押しながら、「X1」の赤い三角ボタンをクリックし、［外れ値の箱ひげ図］の選択を解除します。

Ctrlキーを押しながら「X1」の赤い三角ボタンをクリックし、［ヒストグラムオプション］を選択し、［ヒストグラム］の選択を解除します。

メモ: 応答値はシミュレートされた乱数であるため、検出力のシミュレーションの結果は、検出力（最初の3つの効果）とは異なります。

図10.24 検出力（最初の3つの効果）

「検出力のシミュレーション」アウトラインの「棄却された割合」は、シミュレーションのうち、p値がαより小さくなったものの割合を示します。たとえば、係数が0.8、確率の差が38%である「X3」で、棄却された割合を見てみると、有意水準0.05の検定で379/500 = 0.758となっています。検出力のシミュレーション（有意水準0.05）シミュレーションで計算された、有意水準0.05の検定に対する検出力を、効果ごとに示したものです。「検出したい差」が小さくなるにつれ、検出力も小さくなっています。たとえば、「検出したい差」が24.5%（X6）における検出力は、約0.37にとどまっています。

メモ: 応答値はシミュレートされた乱数であるため、検出力のシミュレーションの結果は、検出力のシミュレーション（有意水準0.05）とは異なります。

表10.1 検出力のシミュレーション（有意水準0.05）
因子	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