Rを使ったD最適計画の選択 - はじめよう実験計画

はじめに

D最適計画については、まず以下の記事をご覧ください。こちらの記事でD最適性の意味を説明しています。今回は、Rを使用してD最適計画を作成する手順について説明します。

sturgeon.hatenablog.com

Rを用いたD最適計画の作り方

D最適計画を作成の流れは、まず候補となる実験点を作成し、その中から指定した回数、指定したモデルに対してD最適性がベストになるような計画を選択する、というもの。

ここでは、要因をA、B、Cの3つとし、それぞれ（-1, 0, 1）の3水準としたばいいのD最適計画を作成しましょう。

以下のコードでは、まず候補点のデータ（data）を作成しています。3水準3要因の全通りなので、27個の候補点があります。ここで、水準に制約がある場合などは、実施可能な実験点のみをdataに格納すればOKです。

> A <- c(-1,0,1)
> B <- c(-1,0,1)
> C <- c(-1,0,1)
> data <- expand.grid(A,B,C)
> data
   Var1 Var2 Var3
1    -1   -1   -1
2     0   -1   -1
3     1   -1   -1
4    -1    0   -1
5     0    0   -1
6     1    0   -1
7    -1    1   -1
8     0    1   -1
9     1    1   -1
10   -1   -1    0
11    0   -1    0
12    1   -1    0
13   -1    0    0
14    0    0    0
15    1    0    0
16   -1    1    0
17    0    1    0
18    1    1    0
19   -1   -1    1
20    0   -1    1
21    1   -1    1
22   -1    0    1
23    0    0    1
24    1    0    1
25   -1    1    1
26    0    1    1
27    1    1    1

次に、Rの"AlgDesign"というパッケージを使って、候補となるdataの中から、1次モデル（~.で表示）に対する試行回数12回のD最適計画を選択します。optFederovという関数は最適計画を選ぶ関数ですが、criterionとしてA最適性やI最適性なども指定できます。

最後の出力結果がD最適計画を示しています。上のコード中の候補点の中から、12点を選んだものだということがわかります。

>library(AlgDesign) > desD <- optFederov(~., data, nTrials=12, criterion = "D") > desD$design Var1 Var2 Var3 1 -1 -1 -1 2 0 -1 -1 3 1 -1 -1 7 -1 1 -1 8 0 1 -1 9 1 1 -1 19 -1 -1 1 21 1 -1 1 22 -1 0 1 24 1 0 1 25 -1 1 1 27 1 1 1

D最適性を利用した計画の拡張

今までは、新しくD最適計画を作成する場合でしたが、追加の実験（拡張計画の作成）を行うときにも、D最適計画の考え方を使うことが出来ます。

どういうことかと言うと、図1のように、事前に行った実験に対し、追加の実験をしたいと思ったとき、全体の計画がD最適になるような追加の実験点を選ぶ、ということです。

f:id:Sturgeon:20200515205655p:plain — 図1. 事前に与えられた計画に対するD最適な拡張

これを行うには、先ほどのoptFederovという関数に、拡張の候補点ともともとの実験点を合わせて渡し、もともとの実験点が拡張後の計画に必ず含まれるように指定します。具体的にはRで以下のように計算します。

pre_design <- read.table("Helicopter_before_aug.txt", header = T)#事前の計画
b <- seq(-1.75,-0.25,by=0.25)
f <- seq(-2,-0.55,by=0.25)
r <- seq(0.36,2.36,by=0.25)
candidate <- expand.grid(b,f,r)#拡張の候補点
colnames(candidate) <- c("B", "F", "R")
X <- rbind(pre_design, candidate)#事前の計画＋拡張の候補点
augment_Dopt <- optFederov(~quad(.), X, nTrials=20, criterion = "D", augment=T, rows = 1:nrow(pre_design))