はじめよう実験計画

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

モデル項の選択のためのPartial F test

 

はじめに

本記事ではPartial F testの話をしたいと思います。

回帰モデルの分散分析では、t検定またはF検定を用いてモデルの係数が0かどうかの検定を行います。Partial F testでは、ここからもう少し踏み込んで、モデルにある項を加える意味があるかを判断できます。つまり、図1のようにModel Aにx3という項を加えてModel Bにすること意味があるかどうか?という問題を、Partial F testで議論することができます。

図1

 

回帰モデルの分散分析

Partial F testの説明に入る前に、基本となる回帰モデルの分散分析について、簡単にまとめてみました。

図2 回帰モデルの分散分析

平方和分解→平方和を自由度で割って平均平方→F検定という流れは、身体に沁みついていてほしいので、まだ沁みてないよ~という方は分散分析の記事をご覧いただければと思います。

図2のF値を計算したら、自由度(k, n-k-1)のF分布の水準αのパーセント点と比べて、計算したF値の方が大きければモデルは有意、という結論になります。

 

www.doe-get-started.com

 

 

 

Partial F test

それではPartial F testの説明に入ります。

説明変数をp個含むモデルAに対し、新たに項を加え、説明変数q個を含むモデルBを作る状況を考えます(図3)。

図3

この時、それぞれモデルの平方和を図的に表すと、図4のようになります。

図4 モデルAとモデルBの平方和分解
SeAとSeBの自由度はそれぞれn-p-1, n-q-1なので、SeA-SeBの自由度はq-p

 

図4の斜線部分は、モデルAに新たな項を加えたことによる残差の減少分であり、大きさはSeA-SeB、自由度はq-pになります。この残差の減少分とSeBを比較して、

 \displaystyle \frac{S_{eA}-S_{eB}/(q-p)}{S_{eB}/(n-q-1)}\approx F_{q-p, n-q-1} \tag{1}

を統計量としたF検定を行うことで、残差の変化に意味があるか、すなわち、モデルAと項数の多いモデルBに有意な差があるかを判定するできます。

 

Partial F testをRで実装

丁寧にモデルごとに残差平和和を計算する方法

式(1)で表されるPartial F testを以下の流れで実装してみます。

  1. データを適当に準備
  2. 項x1のみのモデルAを作ってanova
  3. 項x1とx2のモデルBを作ってanova
  4. anova(モデル)の実行結果のResidualsのSum Sqに残差平方和Seが出力されるので、それを式(1)に代入してF値を計算
  5. 4で求めたF値とF(α=0.05, 自由度1, 自由度1)を比較し、4のF値が大きれば項x2が有意
#データの準備
x1 <- c(-1,-1,1,1)
x2 <- c(-1,1,-1,1)
y <- c(-4,2,-2,5)

###### Partial F test (自分で) #################
modelA <- lm(y~x1)#modelA: y=b0+b1x1
anova(modelA)#実行結果A
SeA <- 42.5#実行結果AのResidualsのSum Sq
modelB <- lm(y~x1+x2)#modelB: y=b0+b1x1+x2
anova(modelB)#実行結果B
SeB <- 0.25#実行結果BのResidualsのSum Sq
F_value <- ((SeA-SeB)/(2-1))/(SeB/1)#式(1)
F_value#F値
qf(0.95, 1,1)#F(1,1)の確率α点

#-------------------------------------------------------#
#実行結果A
> anova(modelA)
Analysis of Variance Table
Response: y
          Df Sum Sq Mean Sq F value Pr(>F)
x1         1   6.25    6.25  0.2941 0.6419
Residuals  2  42.50   21.25
#実行結果B
> anova(modelB)
Analysis of Variance Table
Response: y
          Df Sum Sq Mean Sq F value  Pr(>F)  
x1         1   6.25    6.25      25 0.12567  
x2         1  42.25   42.25     169 0.04887 *
Residuals  1   0.25    0.25
#
> F_value
[1] 169
> qf(0.95, 1,1)
[1] 161.4476
#F_valueを>qf(0.95,1,1)のため、modelAとmodelBの残差の差が有意
#すなわち、modelBで加えた項x2が有効である.

 

上でF_value=169 は qf(0.95,1,1)=161 より大きいため、Modelにx2を項に加えることが意味があると考えられます。

anova (modelA, modelB)で簡単に済ます方法

以上は、計算の過程を示すために自分で実装しましたが、Rでanova(modelA, modelB)を実施することでParital F testは実行できます。

###### Partial F test (簡単に) #################
anova(modelA, modelB)#上の内容
#-------------------------------------------------------#
> anova(modelA, modelB)
Analysis of Variance Table

Model 1: y ~ x1
Model 2: y ~ x1 + x2
  Res.Df   RSS Df Sum of Sq   F  Pr(>F)  
1      2 42.50                           
2      1  0.25  1     42.25 169 0.04887 *
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
> 

 

 

Partial F testを利用したStepwise法の実装

さて、実験で数多くの因子を扱う場合、どの因子をモデルに含めるべきか迷うときがあります。そんな時、Partial F testを使って、項をモデルに加える意味があるかどうか判断できます。この判断を自動で行うのがStepwise法です。

Stepwise法は、項を追加していくForward, 除外していくBackward, 追加と除外を組み合わせた方法がありますが、以下では、ForwardのStepwiseについて説明します。

Partial F testを使用したStepwise法(Foward)は次のフローで実行できます。

  1. .初期モデルの作成 Model0: y~x0
  2. 候補となる項xk, (k=1,…,n) を一つずつ追加してModel: y’ ~ yi + xk を計算
  3. Model y’とyiでPartial F testを行い、p値pk, (k=1,…,n) を計算
  4. 閾値pthより小さなp値の中で最小値に対応する項xkを加え、modelをyi → yi+1に更新
  5. Model yi+1の修正つき赤池情報量基準(AICc)を計算

AICcはモデルの良し悪しの指標です。AICc小さいほど良いモデルと覚えておけばいいと思います。

Stepwise法の実装

以上の内容を実装したのが下の関数selectX_allです。

selectX_all <- function(X_resevoir, Y, enter=0.2){
  selected <- list()
  M <- length(Y)-1-1-1
  df <- data.frame("Y"=Y)
  X_model <- c()
  
  for (j in 1:M) {
    model <- lm(Y~., df)#model
    selected[["params"]][j] <- list(colnames(X_model))
    selected[["AICc"]][j] <- list(AICc(model))
    p <- c()
    for (i in 1:ncol(X_resevoir)) {
      if (j==1) {
        X_add <- data.frame(X_resevoir[i])
      }else{
        X_add <- data.frame(cbind(X_model, X_resevoir[i]))
      }
      df_add <- cbind(X_add, "Y"=Y)
      added.model <- lm(Y~., df_add)#added model
      a <- anova(model, added.model)
      p <- append(p,a$`Pr(>F)`[2])
    }
    p_min <- min(na.omit(p))
    if (p_min<enter) {
      added_term_index <- which(p==p_min)#最小のp値を持つ項のindex
      if (j==1) {
        X_model <- data.frame(X_resevoir[added_term_index])
      }else{
        X_model <- cbind(X_model, X_resevoir[added_term_index])
      }
      X_resevoir <- X_resevoir[-added_term_index]#最小のp値を持つ項を除く
      df <- cbind(X_model, "Y"=Y)
    }else{
      break 
    }
  }
  return(selected) 
}

 

Plackett-Burman計画に対するStepwise法の適用

関数selectX_allを実際に、12回のPlackett-Burman計画で得られたデータに対して使ってみます。

まず12回のPlackett-Burmann計画Xを生成し(要因A,B,C,D,...,L), 12回の試行に対し、仮想的なデータyiを次式で生成します。

 y_i=10+4A+5B+6D+\varepsilon_i

  \varepsilon_i\approx N(0,2^2), i=1,...12

以下のコードでは、データyiから有効な主効果であるA, B, Dを抽出できるか、という実験をしました。

 

library(FrF2)
X <- pb(12, randomize = F)
X_n <- data.frame(apply(X, 2, as.numeric))
colnames(X_n) <- colnames(X)
set.seed(2)
#データをy=10+4A+5B+6D+εで生成
y <- 10+4*X_n$A+5*X_n$B+6*X_n$D+rnorm(12, sd=2)

#Forward-Stepwise
selectX_all(X_n, y)

##############以下、実行結果##############
> X_n
    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
> 
> y
 [1] 15.206171 15.369698 16.175691 22.739249  2.839496 15.264841 -3.584091 -5.479396 10.968948  4.722426 13.835302
[12] 18.963506

#Xとyのデータから、もともとのy=10+4A+5B+6Dを再構築できるか?
> selectX_all(X_n, y)
$params
$params[[1]]
NULL
$params[[2]]
[1] "D"
$params[[3]]
[1] "D" "B"
$params[[4]]
[1] "D" "B" "A" #最小のAIC
$params[[5]]
[1] "D" "B" "A" "E"

$AICc
$AICc[[1]]
[1] 90.91758
$AICc[[2]]
[1] 87.09049
$AICc[[3]]
[1] 81.94675
$AICc[[4]]
[1] 66.9364 #最小
$AICc[[5]]
[1] 70.24708

 

実行結果にピックアップされた主効果とStepごとのモデルのAICcを出力するようにしています。最小のAICcとなるときの主効果はD, B, Aの順になっており、有効な主効果を正しく抽出できることが確認できました。

 

まとめ

モデルに新しい項を追加する意味があるかどうかを判断するPartial F testについて説明しました。Partial F-testを一言でまとめると、Model AとModel Bの残差平方和の差に対するF検定だと言えます。

また、本記事後半ではPartia F testを用いたForward-Stepwise法の紹介をしました。実験では、数多くの因子のうち、どの因子(項)をモデルに加えるかを判断する必要がありますが、この判断をStepwise法によって行うことができるので、その基本となるPartial F testはとっても大事。