DSDにおけるモデル作成
DSDの実験を行ったあと、どのようにしてモデル式を作成するかということは非常に重要な問題です。本記事では、いくつかのモデル作成方法について説明します。
本記事では、DSDの実験で得られる応答Yを仮定して、DSDの実験結果とします。その実験結果に基づき、いくつかの方法でモデルを作成し、出来たモデルがもともとの式を再現しているか、ということを確認します。
計画の作成
DSDの計画作成と応答Yの計算を以下のコードで行います。
#計画Xの作成
library(daewr)#DSDのためのパッケージ X0 <- DefScreen(m=8) X <- X0[1:6]#1~6を列を使用(7,8列目はfake factor)
ここでは、6要因のDSDを作成するにあたり、8要因のDSDを作り、最初の6列のみを使用します。使われなかった2要因はfake factorと呼ばれるもので、後述します。
作成された計画Xは次のような17回の計画になります。
この計画で得られる実験結果Yを以下の式(1)でシミュレートします。
Y=3+2A+4B-C+3D-2AA-2AB+CC+N(0,s=0.3) ・・・(1)
#応答Yのシミュレート #Y=3+2A+4B-C+3D-2AA-2AB+CC+N(0,s=0.3)・・・式(1) y <- function(A,B,C,D){ return(3+2*A+4*B-C+3*D-2*A**2-2*A*B+C**2) } set.seed(1) Y <- y(X$A,X$B,X$C,X$D)+rnorm(numeric(nrow(X)), mean = 0, sd = 0.1)
このシミュレートした応答Yと実験計画Xを一緒に示したものが下の表2です。
シミュレーション実験の結果が準備できました。それでは、これからモデルの作成手法について見ていきましょう。
単純な線形回帰
まずは、「DSDにおいて主効果はほぼ確実である」という特徴を利用して、すべての主効果による1次式を作成してみたいと思います。
> result1 <- lm(Y~A+B+C+D+E+F, data = XY) > summary(result1) Call: lm(formula = Y ~ A + B + C + D + E + F, data = XY) Residuals: Min 1Q Median 3Q Max -3.2512 -2.1299 0.8134 1.8206 1.9378 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.18501 0.62404 3.501 0.005713 ** A 2.00605 0.68766 2.917 0.015374 * B 4.00479 0.68766 5.824 0.000167 *** C -0.97220 0.68766 -1.414 0.187797 D 2.96773 0.68766 4.316 0.001523 ** E 0.02260 0.68766 0.033 0.974433 F -0.03642 0.68766 -0.053 0.958800 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 2.573 on 10 degrees of freedom Multiple R-squared: 0.8631, Adjusted R-squared: 0.781 F-statistic: 10.51 on 6 and 10 DF, p-value: 0.0007925
上記result1のEstimateが各項の係数を表しており、その式は
Y=2.18+2.0A+4.0B-0.97C+2.97D+0.02E-0.03F
となります。もとの式(1)とA,B,C,Dの係数はよく一致していることが分かります。しかし、2次効果を含まないため、モデルはあまりよく当てはまっていません。
そこで、このモデルに2次効果を追加していくことが考えられますが、もとの式(1)を知らない状態では、どの2次効果が有効であるか分かりません。もっともシンプルな方法は、その実験における事前の知識を活用して、有効そうな2次効果をいくつかピックアップして、モデルを作成することです。
例えば、要因Aに対しては、最小値が実験範囲内にあることが分かっているとすれば、下に凸であることを表すAA項をモデルに加えることが考えられます。
ちなみに、全ての2次効果を含んだモデル式を作成することは、線形回帰の特性上できないことに注意が必要です（項の自由度が足りない）。試しに、全ての2次効果を含んだ2次モデル作ろうとすると、以下のようになり、係数は途中までしか計算されません。
Coefficients: (13 not defined because of singularities) Estimate Std. Error t value Pr(>|t|) (Intercept) 2.99838 0.09357 32.044 0.000972 *** A 2.00605 0.02501 80.217 0.000155 *** B 4.00479 0.02501 160.142 3.9e-05 *** C -0.97220 0.02501 -38.876 0.000661 *** D 2.96773 0.02501 118.672 7.1e-05 *** E 0.02260 0.02501 0.904 0.461585 F -0.03642 0.02501 -1.457 0.282547 AB -1.92013 0.04570 -42.017 0.000566 *** AC -0.38199 0.05131 -7.445 0.017569 * AD -0.15196 0.03951 -3.846 0.061427 . AE 0.60977 0.07443 8.193 0.014573 * AF -0.86707 0.07443 -11.650 0.007287 ** BC -0.76973 0.06512 -11.820 0.007082 ** BD -0.20975 0.06512 -3.221 0.084373 . BE NA NA NA NA BF NA NA NA NA CD NA NA NA NA CE NA NA NA NA CF NA NA NA NA DE NA NA NA NA DF NA NA NA NA EF NA NA NA NA AA -0.98767 0.11023 -8.960 0.012228 * BB NA NA NA NA CC NA NA NA NA DD NA NA NA NA EE NA NA NA NA FF NA NA NA NA
情報量基準を用いた解析手法
情報量基準にはAIC, BIC, MICなどがありますが、DSDの初期の論文では、AICによるモデル作成が紹介されています。実際にはAICではなく、サンプル数が少ないことによる影響を補正したAICc（有限補正付き赤池情報量基準）を用いるのが良いと私は考えています。AICというのは次の式で表されるモデルの"良さ"を評価するための指標です。
AICを最小にするモデルが得られた点に対して当てはまりがよく、かつ、オーバーフィッティングでないという意味で"良い"モデルになります。
そこで、可能性のあるモデルに対して1つ1つAICを計算し、AICが最小値になったモデルを採用します。
このプロセスを以下のプログラムに示しました。
#AICcによるモデル選択 > stepwise(XY, y="Y", selection = "forward", select="SL", Choose="AICc", sle=0.2) #実行結果 $process Step EffectEntered EffectNumber Select Choose 1 0 intercept 1 1.000000 78.77538 2 1 B 2 -2.512721 71.15448 3 2 D 3 -2.434319 63.65892 4 3 A 4 -1.948462 58.66953 5 4 AB 5 -3.719537 42.81293 6 5 C 6 -2.447871 34.66278 7 6 AA 7 -3.746816 17.05902 8 7 CC 8 -6.668402 -24.33868 9 8 BE 9 -1.173878 -19.58194 $variate [1] "intercept" "B" "D" "A" "AB" "C" "AA" [8] "CC"
Chooseの列がAICcの計算結果です。Stepが増えるごとにAICcは小さくなっており、最後のStep8で増加し、これ以上の項の追加にモデル性能の向上を見込めないと判断されています。Variateが有効な項であり、上では
intercept (切片), B, D, A, AB, C, AA, CC
の順に効果が大きいことが分かります。これらの項を使用して2次モデルの回帰を行ったのが以下の結果です。
Call: lm(formula = Y ~ B + D + A + AB + C + AA + CC, data = XY) Residuals: Min 1Q Median 3Q Max -0.17073 -0.03612 0.01161 0.07472 0.13738 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.96120 0.08357 35.44 5.61e-11 *** B 4.00479 0.03010 133.03 3.89e-16 *** D 2.96773 0.03010 98.58 5.77e-15 *** A 2.00605 0.03010 66.64 1.95e-13 *** AB -1.98939 0.03506 -56.74 8.26e-13 *** C -0.97220 0.03010 -32.29 1.29e-10 *** AA -1.97777 0.07321 -27.01 6.32e-10 *** CC 1.03525 0.07872 13.15 3.51e-07 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.1126 on 9 degrees of freedom Multiple R-squared: 0.9998, Adjusted R-squared: 0.9996 F-statistic: 5444 on 7 and 9 DF, p-value: 1.383e-15
上の結果より、AICcで作成されたモデル式は
Y=2.96+2.01A+4.00B-0.97C+2.97D-1.99AB-1.97AA+1.03CC ・・・(2)
となり、もとの式(1)とよく一致しています。
効果的モデル選択
前述のように、AICcによるモデル選択はシンプルで、多くの場合良いモデルを与えてくれます。しかし、AICcよりも有効な効果を見落としにくい効果的モデル選択という法が提案されています。
効果的モデル選択はDSDの特徴をフルに活用したモデルの作成方法です。この手法は2つのStepに分かれています。
Step1: 有効な主効果の特定
Step2: 2次効果の推定
Step1
Step1では、fake factorと呼ばれる偽の要因を用いて誤差を推定し、その誤差に対して大きな主効果を有効とします。fake factorは表1に示したように、実験において要因が割り当てられなかった"偽"要因のこと。表1のG列とH列がfake factorです。
fake factorを使うと、実験の推定誤差sfを次式で計算できます。
ここで、X'はXの転置行列、XFは計画行列におけるfake factorの列(表1の列G,列H)、mfはfake factorの個数で今回は2です。
詳細はこちらの論文 https://doi.org/10.1080/00401706.2016.1234979
にかかれていますが、ざっくり言うと、fake factorの実際の効果はゼロなので、fake factorの変動平方和は誤差の推定に使えるということです。
上式の計算をRで書くと、以下のようになります。
> #Fake factorによる誤差の推定 > mf <- 2 > XF <- cbind(X0$G, X0$H) > PF <- XF%*%solve(t(XF)%*%XF)%*%t(XF) > sf <- sqrt(t(Y)%*%PF%*%Y/mf) > sf [,1] [1,] 0.09357043
sfは確かに式(1)で設定した標準偏差sd=0.1とよく一致しています。
sfが求まったら、すべての主効果(1次効果)から有効な主効果を選別します。各効果の係数=0という帰無仮説の下で、
（項の係数/sf）が自由度mfのt分布に従う
という性質を利用し、|項の係数|/sf>t分布の両側20％点・・・★であれば有効な主効果と見なします。そしてそれ以外の主効果は除外します。
以下に実際のプログラムを示します。
result1 <- lm(Y~A+B+C+D+E+F+0, data = XY) barplot(abs(result1$coefficients)/sf[1,1])#図1 abline(h=qt(0.9, 2), col="red")
図1より、★の条件を満たす主効果はA,B,C,Dのみとなります。これは仮定したもとの式(1)と同じであり、正しく有効な主効果を選別できていることが分かります。
Step2
有効な主効果を決めたら、次はStep2「2次効果の推定」です。ふつうは有効な主効果のみで構成される2次項を追加します。今回の例で言えば、A,B,C,Dだけで構成される2次効果AA, AB, AC, ,...,, BB, CC, DDが候補であり、有効でない主効果から構成される2次効果(例えばAFやEF)は候補から除外します。
2次効果を追加する際は、有効な主効果のみのモデルの残差に対して、2次効果当てはめます。2次効果の推定は変数増加法で行っていき、何らかの基準で選べば良いです。以下ではAICcで選んでみました。
> main_effects <- c("A", "B", "C", "D")#有効な主効果 > result_first <- lm(Y~A+B+C+D, data = XY) > X2 <- add_quadratic(X[main_effects], only_additional = T)#main_effects由来の2次効果 > X2$Y2 <- result_first$residuals > stepwise(X2, y="Y2", selection = "forward", select="SL", Choose="AICc", sle=0.2) $process Step EffectEntered EffectNumber Select Choose 1 0 intercept 1 1.000000 44.975546 2 1 AB 2 -6.801765 16.064490 3 2 AA 3 -5.104659 -5.607646 4 3 CC 4 -9.364250 -52.598421 5 4 BD 5 -1.489529 -54.786697 6 5 AD 6 -1.069988 -53.870566 $variate [1] "intercept" "AB" "AA" "CC" "BD"
以上より、有効な2次効果として、AB, AA, CC, BDが選ばれました。intereptも2次効果と同時に評価します。ここで、Step 4, 5, 6においてChoose(AICc)は約1しか変化していないことに注意します。Step2では追加する2次効果の数をできるだけ少なく抑えることが重要です。なぜなら、2次効果が多すぎるとオーバーフィッティングになる危険があるからです。そのため、ここでは、AB, AA, CCのみをモデルに含め、BD, ADはモデルに含めないという判断が出来ます。そうすると、有効な効果は
A, B, C, D, AB, AA, CC
となり、前節の式(2)と一致しますね。
ところで、先の結果からAD, BDを含めることも考えられますが、実際にこれらの効果を含めて回帰分析を行うと...
> result_second <- lm(Y~A+B+C+D+AB+AA+CC+BD+AD, data = XY) > summary(result_second) Call: lm(formula = Y ~ A + B + C + D + AB + AA + CC + BD + AD, data = XY) Residuals: Min 1Q Median 3Q Max -0.09399 -0.03104 -0.01326 0.03787 0.10724 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.01164 0.08151 36.948 2.77e-09 *** A 2.00605 0.02380 84.289 8.71e-12 *** B 4.00479 0.02380 168.271 6.91e-14 *** C -0.97220 0.02380 -40.849 1.37e-09 *** D 2.96773 0.02380 124.696 5.63e-13 *** AB -1.99422 0.03159 -63.128 6.57e-11 *** AA -2.02932 0.06131 -33.099 5.95e-09 *** CC 1.02555 0.07769 13.200 3.35e-06 *** BD 0.06596 0.03119 2.115 0.0723 . AD -0.04668 0.02973 -1.570 0.1604 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.08905 on 7 degrees of freedom Multiple R-squared: 0.9999, Adjusted R-squared: 0.9997 F-statistic: 6776 on 9 and 7 DF, p-value: 6.032e-13
BDとADの係数は小さくP値も大きいため、あまり効果のない項であると考えられます。
DSD用の無料アプリ
こちらで公開しているDSD専用のWebアプリでは、本ページで解説した、AICcを利用したモデル選択と、効果的モデル選択の両方がプログラミングすることなしに実行できます。
https://my-first-dsd.shinyapps.io/DSDApp_ver2/
デモンストレーション動画もそのうち用意する予定です。
