ページタイトルのスペル間違えました、すみません。
誤：Burmann → 正: Burman
Plackett-Burman計画（プラケット・バーマン計画）とは
まずは、Plackett-Burman計画（以下、PB計画）がどのようなものか見てみましょう。以下コードにN=12のPB計画を示します。
> library(FrF2) > pb(12)#pb計画の作成 A B C D E F G H J K L 1 -1 1 1 1 -1 -1 -1 1 -1 1 1 2 1 1 -1 1 1 1 -1 -1 -1 1 -1 3 1 1 1 -1 -1 -1 1 -1 1 1 -1 4 -1 -1 -1 1 -1 1 1 -1 1 1 1 5 1 -1 1 1 -1 1 1 1 -1 -1 -1 6 1 -1 1 1 1 -1 -1 -1 1 -1 1 7 -1 1 1 -1 1 1 1 -1 -1 -1 1 8 1 1 -1 -1 -1 1 -1 1 1 -1 1 9 -1 1 -1 1 1 -1 1 1 1 -1 -1 10 -1 -1 1 -1 1 1 -1 1 1 1 -1 11 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 12 1 -1 -1 -1 1 -1 1 1 -1 1 1 class=design, type= pb
N=12のPB計画は、2水準（±1）でA~Lまで最大11個の要因を含みます。
12回という少ない試行数にも関わらず、多くの要因を調べることができるため、因子のスクリーニングで使用することができます。
しかしながら、因子間の交絡関係（相関関係）が複雑であるというデメリットがあるため、解析には注意が必要です。
本記事は、PB計画を扱うときに理解しておきたい特徴や、解析における注意点についてまとめました。
PB計画の交絡関係
PB計画では、主効果同士はすべて直交しているため（図1）、独立に推定することができます。したがって、主効果のみに注目するときは、特別注意が必要なことはありません。
交互作用を考える際は注意が必要です。例えば、ある主効果Aに対して、Aを含む交互作用AB, AC, AD,...は主効果Aと直交します。しかし、主効果Aを含まない交互作用（BC, DEなど）はAと一部交絡します。
以上のことを、相関係数のカラーマップ（図2）と合わせてまとめますと、
- 主効果同士は直交する
- 主効果AとAを含む交互作用AB,ACなどの関係は直交する。
- 主効果AとＡを含まない交互作用BC, BDなどは相関係数がゼロではない（0.333)
PB計画と他の計画との比較（メリット・デメリット）
さらに12回のPB計画の特徴を理解するため、図3に一部実施要因配置計画と比較してみました。
図3左に示す211-7の一部実施要因配置計画と、N=12のPB計画を比較すると、PB計画の大きなメリットは、主効果と交互作用が完全には交絡しないということです。
何が言いたいかというと、一部実施要因配置計画では主効果の推定をする際、主効果と”完全に”交絡する交互作用が存在するため、主効果と交互作用の分離は基本的にはできないんですね。一方、PB計画では主効果と交互作用が完全には交絡しないため、主効果も交互作用もどちらも検出できるメリットがあるというわけです。
ただ、完全には交絡しないといっても直交しているわけではないので、主効果と交互作用の両方を常に精度よく推定できるわけではありません。この点については具体例を用いて後述します。
PB計画の使いどころ
以上のことを総合すると、PB計画（N=12）の使用時の条件としては
- 主効果は11個まで
- 有効な主効果を含む交互作用（主効果Ａが有効であれば、交互作用AB, ACなど）は検出できる可能性が高い
- 主効果と交互作用は一部交絡が存在するため、あまり多くの効果は検出できない
言葉だけではいまいちピンとこないですよね？次から具体例を用いてPB計画の使用について考えてみようと思います。
PB計画の解析における注意点
本章では、表1示すPB計画を扱います。
この12回のPB計画の実験結果yA, yB, yCはそれぞれ、
で生成したデータです。
以下では、yA, yB, yCそれぞれのケースで、生成された12データから逆に元の式を再構築できるかどうか？を考えていきます。
ケースA（主効果のみが存在する場合）
式yAは主効果A,B,Dのみからなるケースです。
主効果だけを推定したい場合、通常の回帰分析をすればよいですが、要因がA~Lまでの11個あるため、誤差の自由度がなく、分散分析はできません。このような状況で役に立つのがHalf Normal Plotです。
図4は主効果の回帰分析の結果をもとに、Half Normal Plotを作図したもので、D, A, Bが直線から大きく外れているため、有効な主効果であると判断できます。
なお、上図は以下のコードで作図しています。
#PB計画の作成
library(FrF2)
X <- pb(12)
#Xはfactor型なのでnumericに変換
X_n <- data.frame(apply(X, 2, as.numeric))
#行列形式に変換 Xm <- as.matrix(X_n)
#主効果のコントラストを求める contrast_main <- abs(solve(t(Xm)%*%Xm)%*%t(Xm)%*%as.matrix(y1)*2)
#Half Normal Plotの作成 halfnormal(contrast_main, labs = rownames(contrast_main), alpha = 1)
Half Normal Plotについてはこちらの記事をご覧ください。
ケースB（主効果より小さな交互作用が存在する場合）
式yBは主効果A,B,Dに加えて、係数が小さな交互作用AB, ADが存在します。このようなケースでは先述のHalf Normal Plotに加えて、Forward Stepwise法を用いた有効因子の選択ができます。
X2 <- add_interaction(X_n) > selectX_all(X2, y2) $params $params[[1]] NULL $params[[2]] [1] "D" $params[[3]] [1] "D" "B" $params[[4]] [1] "D" "B" "A" $params[[5]] [1] "D" "B" "A" "AD" $params[[6]] [1] "D" "B" "A" "AD" "AB" $params[[7]] [1] "D" "B" "A" "AD" "AB" "AF" #AICc最小値のとき $params[[8]] [1] "D" "B" "A" "AD" "AB" "AF" "DL" $AICc $AICc[[1]] [1] 92.31326 $AICc[[2]] [1] 87.56237 $AICc[[3]] [1] 84.18044 $AICc[[4]] [1] 81.07561 $AICc[[5]] [1] 80.87868 $AICc[[6]] [1] 58.8406 $AICc[[7]] [1] 51.81352 #AICc最小値 $AICc[[8]] [1] 76.18309
上のコードで、selectX_allはPartial F-testに基づいたForward Stepwise法を実行する関数で、こちらの記事で解説しています。
AICcが最も小さいのは上コードの7番目で、D, B, A, AD, AB, AFの順に有効であると判断されています。この効果で回帰分析を実施すると図5のようになります。
ここで、注意が必要なのは交互作用AFです。AFは大きな主効果B、Dと交絡があるため、B, Dに引きずられて検出されている可能性があります。そのため、AFが本当に有効かどうかは、追加の実験で見極める必要があります。とはいえAFはかなり小さい効果なので、yBとほぼ同一のモデルが再現できたと考えて良さそうです。
ケースC（大きな交互作用が存在する場合）
式yCは、主効果A,B,Dと大きさが同等の交互作用が存在する場合で、PB計画において最も注意が必要なケースです。
ケースBと同様に、Forward Stepwiseを実施すると以下のようになります。
X2 <- add_interaction(X_n) > selectX2_all(X2, y3) $params $params[[1]] NULL $params[[2]] [1] "EH" $params[[3]] [1] "EH" "FJ" $params[[4]] [1] "EH" "FJ" "D" $AICc $AICc[[1]] [1] 96.12203 $AICc[[2]] [1] 89.00958 $AICc[[3]] [1] 83.28067 $AICc[[4]] [1] 76.10567 $AICc[[5]] [1] 73.41559
本来存在しないEH, FJが入ってしまっていて、もとの式yCを全く再現できていないことが分かります。
これは、図2で説明したように、主効果と関係のない交互作用がその主効果と交絡するために起きる現象です。つまり、本来有意でない交互作用項の変動が、実際には有効な主効果の変動を"消費"してしまうので、Forward Stepwiseで一度意味のない効果が取り込まれてしまうと、間違ったモデルが選択されてしまうのです。
では、どうすればよいか？
もとのyCを再現するには、実験結果に対して以下の前提が成り立つと考えます。
- 交互作用より主効果の方が大きく、優先的に評価するべきである
- 有効な交互作用は有効な主効果からなる（Heredityという）
2に関して、例えば主効果Aが有効ならAB, ACなどは有効な交互作用の候補となりうるが、EHやFJなどは候補とならないと仮定する、ということです。
この前提が妥当であるとすれば、まずは、ケースAのようにHalf Normal Plotで有効な主効果を特定し、次に有効な主効果を含む交互作用の中から項を選ぶという順序でモデル選択を行います。
図6においても、やはりD,A,Bが有効な因子であることが分かります。そして、K,C,J,EとL,F,G,Hが2グループに分かれていることが見て取れますが、このように2グループに分かれることは、交互作用項が存在することを示唆しています。*1
次に、有効な主効果D,A,Bと、関連する交互作用AB, AD, BDで回帰分析を行うと図7のようになり、AB, ADをモデルに加えればよいことが分かり、式yCが再現できました。
交互作用項がたくさんあって、分散分析ができない（誤差の自由度がゼロ）場合には、候補となる交互作用だけを用いて、ケースBのように、Forward-Stepwiseしてもよいです。しかしながら、前述のように、Foward Stepwiseでは意味のない交互作用が一度含まれると、本来有効な効果の変動を消費してしまって、有効な効果を検出できない場合があります。そのため、できればFoward Stepwiseではなくすべてのモデルの組合せを総当たりで調べる方が好ましいとされています。*1
まとめ
長くなりました、本記事のまとめです。
・有効な主効果は高い確率で検出できる。
・交互作用項は主効果と一部交絡するため、通常のForward Stepwiseだと関係ない交互作用が誤って検出される恐れがある。
・このリスクを下げるため、まずHalf Normal Plotで有効な主効果を判定し、
・その主効果を含む交互作用項の候補からモデルを選択する。
*1:Analysis of the 12 run Plackett-Burman design, John Tyssedal and Oddgeir Samset, 1998, https://www.researchgate.net/publication/2284923_Analysis_of_the_12_run_Plackett-Burman_design