D最適な計画 - はじめよう実験計画

はじめに
最小二乗法の係数
分散共分散行列の例
D最適とは
まとめ
D最適計画による計画の拡張（実例）
Rを使用したD最適計画の選択
D最適計画について詳しく理解するための参考書

はじめに

D最適計画はある試行回数・指定したモデルについて、もっとも効率的に効果の推定ができる計画です。

この性質を利用すると、図1のように、最初の実験の後に追加で実験をする際「指定した回数で想定したモデルの項を推定するのに、もっとも効率的な実験点」を追加することが出来ます。実験の拡張だけではなく、最初からD最適計画を用いて実験計画を作成することも可能です。

本記事では、基本となる最小二乗法と分散共分散行列について簡単にまとめてから、D最適計画の理論について説明したいと思います。
また、Ｒを用いた具体的なＤ最適計画の作成方法についてはこちらの記事をご覧ください。
sturgeon.hatenablog.com

最小二乗法の係数

はじめに最小二乗法の行列表式を書いておきます。ちょっとややこしいかもしれませんが、D最適計画を立てる際に必要となる式なので、一応頭の片隅に置いといてください。
まずは回帰モデルを

$y_i=\beta_0+\beta_1 x_{i1}+\beta_2 x_{i2}+\beta_3 x_{i3}+\cdots+\beta_k x_{ik}+\varepsilon_i, i=1,...,n$

とおきます。y_iが実験値、x_iが変数（要因の値）です。nは観測の数、kはモデルに含まれる項の数。この式を行列で表現すると

$\boldsymbol{y}=\boldsymbol{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}$

と書けます。ここで、

$\boldsymbol{y}=\begin{bmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{n} \end{bmatrix},\ \boldsymbol{X}=\begin{bmatrix} 1 & x_{11} & x_{12} &\cdots & x_{1k} \\ 1 & x_{21} & x_{22} &\cdots & x_{1k} \\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_{n1} & x_{n2} &\cdots & x_{1k} \end{bmatrix}, \ \boldsymbol{\beta}=\begin{bmatrix} \beta_{1} \\ \beta_{2} \\ \vdots \\ \beta_{n} \end{bmatrix},\ \boldsymbol{\varepsilon}=\begin{bmatrix} \varepsilon_{1} \\ \varepsilon_{2} \\ \vdots \\ \varepsilon_{n} \end{bmatrix}$

で、Xは計画行列とかデザインマトリックスとか呼ばれる行列です（後でまた出てきます）。
上式中の係数βが最小二乗法で推定するべきもので、

$\hat{\boldsymbol{\beta}}=(\boldsymbol{X}'\boldsymbol{X})^{-1}\boldsymbol{X}'\boldsymbol{y}'$

で計算することができます。
さらに、天下り的ではありますがβの分散共分散行列という行列がありまして、

$Cov(\hat{\boldsymbol{\beta}})=\sigma^2(\boldsymbol{X}'\boldsymbol{X})^{-1}$

で計算されます。

分散共分散行列の例

分散共分散行列の対角成分はモデルの項β_iの分散であり、非対角成分はβ_iとβ_jの共分散を表しています。ここでは、具体的なイメージを持ってもらうために、下の3要因配置計画に対して、分散共分散行列を実際に計算してみましょう。

実験で得られる値yを、変数をxとして、モデル式

$\boldsymbol{y}=\boldsymbol{\beta_0}+\boldsymbol{\beta_1 x_1}+\boldsymbol{\beta_2 x_2}+\boldsymbol{\beta_3 x_3}$

にフィットさせることを考えると、はじめに説明したデザインマトリックスは

$\boldsymbol{X}=\begin{bmatrix} 1 & -1 & -1 & -1 \\ 1 & 1 & -1 & -1 \\ 1 & -1 & 1 & -1 \\ 1 & 1 & 1 & -1 \\ 1 & -1 & -1 & 1 \\ 1 & 1 & -1 &1\\ 1 & -1 & 1 & 1 \\ 1& 1 &1 &1 \end{bmatrix}$
となります。したがって、分散共分散行列は

$\sigma^2(\boldsymbol{X'X})^{-1}=\sigma^2\begin{bmatrix} 1/8 & 0 & 0 & 0 \\ 0 & 1/8 & 0 & 0 \\ 0 & 0 & 1/8 & 0 \\ 0 & 0 & 0 & 1/8 \end{bmatrix}$
と計算できます。この分散共分散行列は対角成分のみに値があり、非対角成分はすべてゼロになっています。このことはβ₁とβ₂が完全に独立に推定できることを意味しています。これは要因配置計画では常に成り立つ性質です。

D最適とは

前章で見たように、分散共分散行列はモデルの項が独立に推定可能かどうか（交絡がない）を示してくれます。先の例は主効果のみのモデル式でしたが、交互作用項、2乗項が加わっても同じように分散共分散行列を計算すれば各項が独立かどうかわかります。
したがって、分散共分散行列の成分ができるだけゼロに近いほど、交絡が少ない計画ということになります。このような考えから、

■D最適計画■

$det((\boldsymbol{X}'\boldsymbol{X})^{-1})$ を最小にする計画XをD最適計画と呼び、もっともモデル項の交絡が少ない効率的な計画である。

と言えます。「Xの行列式を最小にする」というのが「D」最適の意味で、限られた試行回数の中で、モデルの項の推定がもっとも正確になります。D最適以外にも、行列式以外を基準としたA最適やI最適など、さまざまな最適計画が存在しますが、D最適がもっともメジャーな計画です。

まとめ

以上でD最適計画の理論的な話はおしまいです。結構分かりにくい概念だと思うので、もう一度イメージを確認しておきますね（図1再掲）。

はじめに図1の左の計画表があります。この計画に従って実験をして、追加の実験が必要になったとします。「その追加の実験点を含めた20回の計画Xは、要因B,F,Rの1次モデルに対して、 $det((\boldsymbol{X}'\boldsymbol{X})^{-1})$ が最も小さくなるとき、モデルの推定精度が最も高くなる」というわけです。

このようにして作成した計画を、3次元(B,F,R)の実験空間にプロットしたものが図2です。最初の点が青、D最適計画によって拡張された点が赤で示されています。

D最適計画による計画の拡張（実例）

こちらの記事ではD最適計画によりPlackett-Burmann計画の拡張を行いました。
sturgeon.hatenablog.com

Rを使用したD最適計画の選択

Rを使ってＤ最適計画が出来ます。詳しくはこちらをご覧ください。
sturgeon.hatenablog.com

D最適計画について詳しく理解するための参考書

D最適計画については『実験計画法-方法編-』（山田　秀著）でも説明されています。ただ、例題がないのと、実際の作成手順については記述されてないです。
一方、"Design and Analysis of Experiment 10th edition"は例題も豊富なのでおすすめです。拡張された点が実験空間上でどのように分布するかなど、図を使った説明も多く、英語ではありますがかなり理解が深まると思います。

実験計画法方法編―基盤的方法から応答曲面法、タグチメソッド、最適計画まで

amazon

Design and Analysis of Experiments, Wiley; 第10版 (2020/6/23)

amazon