はじめに
本記事ではPartial F testの話をしたいと思います。
回帰モデルの分散分析では、t検定またはF検定を用いてモデルの係数が0かどうかの検定を行います。Partial F testでは、ここからもう少し踏み込んで、モデルにある項を加える意味があるかを判断できます。つまり、図1のようにModel Aにx3という項を加えてModel Bにすること意味があるかどうか？という問題を、Partial F testで議論することができます。
回帰モデルの分散分析
Partial F testの説明に入る前に、基本となる回帰モデルの分散分析について、簡単にまとめてみました。
平方和分解→平方和を自由度で割って平均平方→F検定という流れは、身体に沁みついていてほしいので、まだ沁みてないよ～という方は分散分析の記事をご覧いただければと思います。
図2のF値を計算したら、自由度(k, n-k-1)のF分布の水準αのパーセント点と比べて、計算したF値の方が大きければモデルは有意、という結論になります。
Partial F test
それではPartial F testの説明に入ります。
説明変数をp個含むモデルAに対し、新たに項を加え、説明変数q個を含むモデルＢを作る状況を考えます（図3）。
この時、それぞれモデルの平方和を図的に表すと、図4のようになります。
SeAとSeBの自由度はそれぞれn-p-1, n-q-1なので、SeA-SeBの自由度はq-p
図4の斜線部分は、モデルAに新たな項を加えたことによる残差の減少分であり、大きさはSeA-SeB、自由度はq-pになります。この残差の減少分とSeBを比較して、
を統計量としたF検定を行うことで、残差の変化に意味があるか、すなわち、モデルＡと項数の多いモデルBに有意な差があるかを判定するできます。
Partial F testをRで実装
丁寧にモデルごとに残差平和和を計算する方法
式(1)で表されるPartial F testを以下の流れで実装してみます。
- データを適当に準備
- 項x1のみのモデルAを作ってanova
- 項x1とx2のモデルBを作ってanova
- anova(モデル)の実行結果のResidualsのSum Sqに残差平方和Seが出力されるので、それを式(1)に代入してF値を計算
- 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）は次のフローで実行できます。
- .初期モデルの作成 Model0: y~x0
- 候補となる項xk, (k=1,…,n) を一つずつ追加してModel： y’ ~ yi + xk を計算
- Model y’とyiでPartial F testを行い、p値pk, (k=1,…,n) を計算
- 閾値pthより小さなp値の中で最小値に対応する項xkを加え、modelをyi → yi+1に更新
- 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を次式で生成します。
以下のコードでは、データ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はとっても大事。
