はじめよう実験計画

実験を早く終わらせるための技術

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))

 

上記コードでは、まずpre_designに事前の計画を格納し(ここでは、"Helicopter_before_aug.txt”というファイルから読み取る)、拡張の候補点と合わせた計画XをoptFederovに渡しています。

D最適計画を選択する際の、モデル式は二次モデルで、試行回数は20回としています。全く新しくD最適計画を作成するときと異なるのは、augment=Tとして、事前の計画の行番号をrowsで指定することです。こうすることで、事前の計画の実験点が拡張後も含まれるようになります。

このようにして、拡張した計画augment_Doptを表2に示しました。

f:id:Sturgeon:20200515210729p:plain

表2. 拡張後のD最適計画

最後に

本記事はRを使ったD最適計画の作成方法を紹介しました。特に事前の計画から拡張するという手順は、どんな実験にも応用が効く方法なので、非常に便利だと思います。

理論をしっかりと抑えたいという方には以下の参考書がおすすめです。

 

sturgeon.hatenablog.com