正しいモデル選択をするためにモデル評価が大事になってきます。いくつか方法がありますが、代表的なものを列挙します。
- Hold Out法
- Cross Validation法
- Leave One Out法
- Bootstrap法
- WAIC
- WBIC
leave one outはcross validationにおいて、データ数とグループ数が等しい場合を指します。まずはこの方法を使って、性能評価を行い、正しくモデル選択ができるか検証してみます。
テストデータの作成
下の因果グラフに基づき、テストデータを生成します。誤差項はすべて、正規分布に従うとします。1回のサンプリング数は10とします。
n= 10
mu = c(0,0)
sd = c(1,1)
cor =.2
cov = sd %*% t(sd) * cor
diag(cov) = sd^2
generate_data = function(){
mvrnorm(n,mu,cov) %>% as.tibble() %>% mutate(V3 = rnorm(n,0,1),y=3 * V1 + rnorm(n,0,.5))
}
作成されたデータの概観をみてみましょう。
モデルの候補
1. Y ~ V1 (正しいモデル)
1. Y ~ V1+V2
1. Y ~ V1+V2+V3
1. Y ~ (V1+V2)^2 (相互作用も考慮したモデル)
1. Y ~ (V1+V2+V3)^2
Leave one out法
では、Leave one outに基づいてモデルを選択してみます。
models = list(y~V1, y~V1+V2, y~V1+V2+V3)
result = models %>% map(~{
model = .x
1:nrow(df) %>% map(
~(lm(model, data=df[-.x,]) %>% predict(df[.x,]) - df$y[.x])^2
)
})
result %>% simplify_all() %>% set_names(str_c("model_",1:3)) %>% as.tibble %>% summarise_all(mean) %>% which.min
# >>model_1
上記のデータセットでは正しくモデル1が選ばれたようです。
次に上記の流れ(データ生成+モデル選択)を1000回繰り返します。果たして正しいモデルを選択できたのでしょうか?
models = list(y~V1, y~V1+V2, y~V1+V2+V3, y~(V1+V2)^2, y~(V1+V2+V3)^2)
1:100 %>% map(~{
df = generate_data()
result = models %>% map(~{
model = .x
1:nrow(df) %>% map(
~(lm(model, data=df[-.x,]) %>% predict(df[.x,]) - df$y[.x])^2
)
})
result %>% simplify_all() %>% set_names(1:length(models)) %>% as.tibble %>% summarise_all(mean) %>% which.min
})
result %>% unlist %>% as_tibble %>%
ggplot(aes(value)) + geom_bar() + labs(x="モデル",y="カウント数")
さて、気になる結果です。一番選択されやすいモデルは正しいモデル①であったものの、それ以外のモデルも少なからず選択されています。
| Model | Count |
|---|---|
| 1 | 639 |
| 2 | 130 |
| 3 | 101 |
| 4 | 96 |
| 5 | 34 |
次にサンプルサイズを100に増やして、同じように実験してみます。
結果はこちら。
| Model | Count |
|---|---|
| 1 | 701 |
| 2 | 95 |
| 3 | 89 |
| 4 | 76 |
| 5 | 39 |
想像どおり、サンプル数を増やすことによって、より正しくモデル選択ができるようになりました。
しかし、それでも正答率は70%ほど。データのみから正しいモデルを選択することはかなり困難な印象です。Leave one out以外にWAICも汎化誤差の近似として使いますが、性能は似たようなものらしいです(うろ覚え)。
今回は正規分布かつ線形モデルでテストデータを作成しました。
しかし、現実は背景因子の確率分布も複雑であり、因果関係も線形構造ではないことが多いでしょう。
やはり、その道の専門家と手を取り合い、少しでも多くのドメイン知識に基づいて、モデルを作成していく&候補を絞っていく作業が重要になってくるのではないでしょうか。
Bootstrap法(追記)
1つの標本からリサンプリングを繰り返して生成する標本をブートストラップ標本と呼びます。
手元のデータから最尤推定で求めたパラメタでMSEを計算しても、真のデータ分布を使ったMSEはそれより大きくなる。それをブートストラップ標本を使って修正する方法です。詳しくは「はじパタ第2章」参照してください。
result = 1:100 %>% map(~{
df = generate_data()
models %>% map(~{
model=.x
expected_bias = 1:100 %>% map_dbl(~{
sample = df %>% sample_frac(replace=T)
fit = lm(model,data=sample)
fitted_MSE = fit$residuals^2 %>% mean
predicted_MSE = (predict(fit, df) - df$y)^2 %>% mean
bias = predicted_MSE-fitted_MSE
bias
}) %>% mean
lm(model,df)$residuals %>% .^2 %>% mean + expected_bias
}) %>% which.min
})
気になる結果は。。。
| Model | Count |
|---|---|
| 1 | 622 |
| 2 | 164 |
| 3 | 134 |
| 4 | 72 |
| 5 | 8 |
今回のテストデータだけからは結論付けれませんが、LOOVと比較すると性能が悪いような気がします。