AIC_PMSE.Rmd

---
title: "いつでもAICを使えばよいというものではない"
author: "伊東宏樹"
date: "2022-06-07"
output:
   revealjs::revealjs_presentation:
    self_contatined: false
    theme: black
    transition: slide
    css: style.css
    reveal_options:
      slideNumber: true
      previewLinks: false
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(MASS)
library(AICcmodavg)
library(parallel)
# Set number of CPU cores
options(cl.cores = 8)
```

# 生態学の統計解析で<br />わりとあるパターン {#intro}

::: {style="margin-top: 2em; font-size: 125%;"}
1.  AIC（赤池情報量基準）でモデル選択
2.  「ベストモデル」で、係数の有意性検定
:::

## しかし {#intro2}

[ そもそも ]{style="font-size: 150%;"}

[ **「AICは正しいモデルを選ぶためのものではない」** ]{style="font-size: 200%;"}

::: {style="text-align: right;"}
粕谷 (2015)
:::

## AICは予測性能が高いモデルを選ぶ {#intro3}

-   「AICの理論は平均的な予測がよいモデルを相対的に選ぶもの」
-   「AIC最小のモデルが真のモデルとどれだけ確実に一致するか、あるいは他のモデルが否定されるかについて述べているわけではない」

::: {style="text-align: right;"}
粕谷 (2015)
:::

## モデルの当てはまりの良さでもない {#intro4}

::: {style="text-align: left; padding: 1em;"}
「モデルの相対的な当てはまりの良さ」を「AICにより比較した」という論文もあったりするが... <!-- 日林誌 vol.103 p.274-->
:::

::: {style="text-align: left;"}
-   「AICは統計モデルのあてはまりの良さ (goodness of fit) ではなく、予測の良さ (goodness of prediction) を重視するモデル選択基準です」(久保 2012)
-   「モデルの当てはまりと予測性能は違う」(伊庭 2017)
:::

# モデル選択後に検定? {#post-selection-test}

::: {style="text-align: left;"}
-   有意性検定: 多くのばあい、係数が0であるかないかを検定
-   モデル選択で はずれた説明変数は、係数=0としていることと等価
    -   残った変数で、あらためて係数=0か検定する意義とは?
:::

## 回帰分析の目的 {#regression}

::: {style="text-align: left;"}
1.  目的変数と説明変数の関係の記述
2.  説明変数による目的変数の予測
3.  目的変数に対する説明変数の介入効果の推定
:::

::: {style="text-align: right;"}
竹内ほか (2022)
:::

## 予測と因果はちがう {#predition-causual}

::: {style="font-size: 125%; text-align: left;"}
「予測を目的とした回帰モデルの偏回帰係数に因果関係的な解釈を持ち込もうとする(期待する)のはもちろんのこと，個々の偏回帰係数が生物学的・生態学的にどのような意味を持つのかを考察することも，本来の解析目的からは逸脱した，明確な理論的根拠に欠ける行為である」
:::

::: {style="text-align: right;"}
竹内ほか (2022)
:::

## モデル選択後の検定・推定 {#pms}

::: {style="text-align: left;"}
-   AICに限らず、選択されたモデルという条件付きでは
    -   推定値の分布がゆがむことがある
    -   帰無仮説が棄却される確率が 有意水準どおりにならないことがある
-   正しく推定できる手法は、統計学者が研究している (Kuchibhotla et al. 2022, Zhang et al. 2022, etc.)
:::

# 実例1 {#example1}

```{r example1}
# Number of replications
R <- 10^4

# Sample size
N <- 1000

set.seed(123)
sd1 <- 1
sd2 <- 3
```

::: {style="text-align: left;"}
-   サンプルサイズ: `r N`
-   説明変数: $x_1, x_2, x_3 \sim \mathrm{Normal}(0, `r sd1`)$
-   目的変数: $y \sim \mathrm{Normal}(0, `r sd2`^2)$

目的変数は説明変数とは無関係
:::

```{r ex1_data, cache=TRUE}
data <- lapply(1:R, function(i) {
  x <- mvrnorm(N, c(0, 0, 0),
               matrix(c(sd1, 0, 0,
                        0, sd1, 0,
                        0, 0, sd1),
                      3, 3))
  y <- rnorm(N, 0, sd2)
  data.frame(x1 = x[, 1], x2 = x[, 2], x3 = x[, 3], y = y)
})
```

::: {style="text-align: left; margin-top: 2em; font-size: 80%;"}
※この例は、Kuchibhotla et al. (2022) を参考にしました。
:::

## モデル {#ex1_model}

1.  `y ~ 1       ←「正しい」モデル`
2.  `y ~ x1`
3.  `y ~ x2`
4.  `y ~ x3`
5.  `y ~ x1 + x2`
6.  `y ~ x1 + x3`
7.  `y ~ x2 + x3`
8.  `y ~ x1 + x2 + x3`

::: {style="text-align: left;"}
-   AICによるモデル選択を 10\^`r log10(R)` 回くりかえし<br />
-   選択されたモデルの頻度分布を求める
:::

```{r ex1_model, cache=TRUE}
cl <- makeCluster(getOption("cl.cores", 2))
rlist <- clusterMap(cl, function(d) {
  m <- list()
  m[[1]] <- lm(y ~ 1, data = d)
  m[[2]] <- lm(y ~ x1, data = d)
  m[[3]] <- lm(y ~ x2, data = d)
  m[[4]] <- lm(y ~ x3, data = d)
  m[[5]] <- lm(y ~ x1 + x2, data = d)
  m[[6]] <- lm(y ~ x1 + x3, data = d)
  m[[7]] <- lm(y ~ x2 + x3, data = d)
  m[[8]] <- lm(y ~ x1 + x2 + x3, data = d)

  aic <- sapply(1:8, function(i) AIC(m[[i]]))
  aicc <- sapply(1:8, function(i) AICcmodavg::AICc(m[[i]]))
  bic <- sapply(1:8, function(i) BIC(m[[i]]))
  coef <- sapply(1:8, function(i) coef(m[[i]]))

  # coefficient of x1 and its p-value in the smallest-AIC model
  min_aic <- which.min(aic)
  if (min_aic %in% c(2, 5, 6, 8)) {
    s <- summary(m[[min_aic]])
    coef_x1 <- s$coefficients["x1", "Estimate"]
    pv_x1 <- s$coefficients["x1", "Pr(>|t|)"]
  } else {
    coef_x1 <- pv_x1 <- NA
  }
  
  # p-value for coefficient of x1 in the full model (model 8)
  pv_x1_m8 <- summary(m[[8]])$coefficients["x1", "Pr(>|t|)"]
  list(AIC = aic, AICc = aicc, BIC = bic, coef_x1 = coef_x1,
       pv_x1 = pv_x1, pv_x1_m8 = pv_x1_m8)
}, data)
stopCluster(cl)
```

## AIC {#ex1_AIC}

各モデル（1〜8）が選択された割合

```{r ex1_AIC}
# smallest-AIC model
min_aic <- sapply(1:R, function(i) which.min(rlist[[i]]$AIC))
pr <- summary(factor(min_aic)) / R
print(pr[1:4])
print(pr[5:8])
```

「正しい」モデル（model1）が選ばれた割合はおよそ`r round(summary(factor(min_aic))[1] / R * 100, 0)`%

## x1の係数の頻度分布 {#ex1_coef_x1}

x1を含むモデルが選択されたばあい

```{r ex1_coef_x1}
# coefficient of x1 in the smallest-AIC models which contain x1
coef_x1 <- sapply(1:R, function(i) {
         rlist[[i]]$coef_x1
  })
coef_x1 <- coef_x1[!is.na(coef_x1)]
hist(coef_x1, main = "In the \"best\" models which contain x1",
     xlab = "Coefficient of x1")
```

## x1の係数についてのp値の頻度分布 {#ex1_p-value_x1}

x1を含むモデルが選択されたばあい

```{r ex1_p-value_x1}
# p-value for x1 in the smallest-AIC models which contain x1
pv_x1 <- sapply(1:R, function(i) {
  rlist[[i]]$pv_x1
})
pv_x1 <- pv_x1[!is.na(pv_x1)]
hist(pv_x1, main = "In the \"best\" models which contain x1",
     xlab = "p-value for the coefficientnt of x1",
     xlim = c(0, 1), breaks = seq(0, 1, 0.0125))
abline(v = 0.05, col = "red", lty = 2)
text(0.6, 100,
     paste("Prop. p < 0.05:", round(sum(pv_x1 < 0.05) / length(pv_x1), 2)),
     cex = 2)
```

## モデル選択をしないとき {#ex1_p-value_full}

フルモデルでの、x1の係数についてのp値

```{r ex1_p-value_x1_in_model8}
# p-value for x1 in the full model
pv_x1_m8 <- sapply(1:R, function(i) {
  rlist[[i]]$pv_x1_m8
})
hist(pv_x1_m8, main = "In the full model",
     xlab = "p-value for the coefficientnt of x1",
     xlim = c(0, 1), breaks = seq(0, 1, 0.0125))
abline(v = 0.05, col = "red", lty = 2)
```

## BIC {#ex1_BIC}

各モデル（1〜8）が選択された割合

```{r ex1_BIC}
min_bic <- sapply(1:R, function(i) which.min(rlist[[i]]$BIC))
pr <- summary(factor(min_bic)) / R
pr2 <- rep(0, 8)
names(pr2) <- 1:8
for (i in names(pr)) {pr2[as.numeric(i)] <- pr[i]}
print(pr2[1:4])
print(pr2[5:8])
```

「正しい」モデル（model1）が選ばれた割合はおよそ`r round(summary(factor(min_bic))[1] / R * 100, 0)`%

BICは、候補に「正しい」モデル（データを生成したモデル）を含んでいて、 サンプルサイズがじゅうぶんに大きければ、「正しい」モデルを選択するが...

# 実例2 {#example2}

```{r ex2_settings}
R <- 10^4
N <- 50

set.seed(123)
var1 <- 1
var2 <- 9
cov <- 1
coef <- 0.3
sd1 <- 4
```

::: {style="text-align: left;"}
-   サンプルサイズ: `r N`  （小サンプル）
-   説明変数: $x_1, x_2$
-   目的変数: $y$

以下のような式でデータが生成される。
:::

$$
\begin{pmatrix}x_1 \\ x_2 \end{pmatrix} \sim
\mathrm{MultiNormal}\left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} `r var1` & `r cov` \\ `r cov` & `r var2` \end{pmatrix} \right)
$$

$$
y \sim \mathrm{Normal}(x_1 + `r coef` x_2, `r sd1`^2)
$$

::: {style="text-align: left; margin-top: 2em; font-size: 80%;"}
※この例は、大久保 (2022) を参考にしました。
:::

```{r ex2_data, cache=TRUE}
data <- lapply(1:R, function(i) {
  x <- mvrnorm(N, c(0, 0), matrix(c(var1, cov, cov, var2), 2, 2))
  y <- x[, 1] + coef * x[, 2] + rnorm(N, 0, sd1)
  data.frame(x1 = x[, 1], x2 = x[, 2], y = y)
})
```

## データの例 {#ex2_data}

```{r ex2_view_data}
pairs(data[[1]])
```

理論的な<abbr title="Variance Inflation Factor">VIF</abbr> = `r 1 / (1 - cov^2 / (var1 * var2))`

## モデル {#ex2_model}

1.  `y ~ x1`
2.  `y ~ x1 + x2   ←「正しい」モデル`

::: {style="text-align: left;"}
-   AICによるモデル選択を 10\^`r log10(R)` 回くりかえし<br />
-   選択されたモデルの頻度分布を求める
:::

```{r ex2_model, cache=TRUE}
cl <- makeCluster(getOption("cl.cores", 2))
rlist <- clusterMap(cl, function(d) {
  m1 <- lm(y ~ x1, data = d)
  m2 <- lm(y ~ x1 + x2, data = d)

  c(AIC1 =  AIC(m1), AIC2 = AIC(m2),
    AICc1 = AICcmodavg::AICc(m1), AICc2 = AICcmodavg::AICc(m2),
    BIC1 =BIC(m1), BIC2 = BIC(m2),
    coef1 = coef(m1)[-1], coef2 = coef(m2)[-1])
}, data)
stopCluster(cl)

vars <- names(rlist[[1]])
results <- t(matrix(unlist(rlist), nrow = length(vars)))
colnames(results) <- vars
```

## AIC {#ex2_AIC}

Model 2 (「正しい」モデル) が選ばれた割合

```{r ex2_AIC}
# AIC
sum(results[, "AIC1"] > results[, "AIC2"]) / R
```

およそ`r round(sum(results[, "AIC1"] > results[, "AIC2"]) / R * 100, 0)`%

## AICc {#ex2_AICc}

Model 2 (「正しい」モデル) が選ばれた割合

```{r ex2_AICc}
# AICc
sum(results[, "AICc1"] > results[, "AICc2"]) / R
```

およそ`r round(sum(results[, "AICc1"] > results[, "AICc2"]) / R * 100, 0)`%

## BIC {#ex2_BIC}

Model 2 (「正しい」モデル) が選ばれた割合

```{r ex2_BIC}
# BIC
sum(results[, "BIC1"] > results[, "BIC2"]) / R
```

およそ`r round(sum(results[, "BIC1"] > results[, "BIC2"]) / R * 100, 0)`%

## AICでmodel1が選択されたばあい {#ex2_mod1}

x1の係数（真値=1）: 大きいほうにかたよる。

```{r ex2_m1_x1}
m1_x1 <- results[results[, "AIC1"] < results[, "AIC2"],
                 "coef1.x1"]
hist(m1_x1, main = "In case model1 is selected",
     xlab = "Coefficient of x1")
abline(v = 1, col = "red", lty = 2)
```

## model2が選択されたばあい {#ex2_mod2}

x1の係数（真値=1）: 小さいほうにかたよる。

```{r ex2_m2_x1}
# in case model2 is selected
m2_x1 <- results[results[, "AIC1"] > results[, "AIC2"],
                 "coef2.x1"]
hist(m2_x1, main = "In case model2 is selected",
     xlab = "Coefficient of x1")
abline(v = 1, col = "red", lty = 2)
#sum(m2_x1 > 1) / length(m2_x1)
```

------------------------------------------------------------------------

x2の係数（真値=`r coef`）: 大きいほうにかたよる。

```{r ex2_m2_x2}
m2_x2 <- results[results[, "AIC1"] > results[, "AIC2"],
                 "coef2.x2"]
hist(m2_x2, main = "In case model2 is selected",
     xlab = "Coefficient of x2")
abline(v = coef, col = "red", lty = 2)
```

# まとめ {#summary}

-   AICは予測を最適化するもの
-   モデル選択後に係数の検定をするのはアブない
-   説明・予測・介入（因果）のなにを問題としているのか意識するようにする

# 参考文献 {#references}

::: {style="font-size: 66%; text-align: left; line-height: 100%;"}
-   伊庭幸人 (2017) モデル選択超速習 AICからスパースまで. 岩波データサイエンス Vol.5 p.6--18
-   粕谷英一 (2015) [生態学におけるAICの誤用: AICは正しいモデルを選ぶためのものではないので正しいモデルを選ばない.](https://doi.org/10.18960/seitai.65.2_179) 日本生態学会誌 65:179--185
-   久保拓弥 (2012) データ解析のための統計モデリング入門---一般化線形モデル・階層ベイズモデル・MCMC. 岩波書店
-   Kuchibhotla et al. (2022) [Post-Selection Inference.](https://doi.org/10.1146/annurev-statistics-100421-044639) Annual Review of Statistics and Its Application 9:505--527
-   大久保祐作 (2022) "生態学と因果推論1" <https://github.com/OhkuboYusaku/blog/tree/main/causal>
-   竹下ほか (2022) [哺乳類研究における非実験データセットに対する回帰分析の適用について.](https://doi.org/10.11238/mammalianscience.62.69) 哺乳類科学 62: 69--79
-   Zhang et al. (2022) [Post-model-selection inference in linear regression models: An integrated review.](https://doi.org/10.1214/22-SS135) Statistics Surveys 16: 86--136
:::

## ファイルなど {#files}

::: {style="text-align: left;"}
今回のスライドは

<https://ito4303.github.io/AIC_PMSE.html>

にあります。
:::

::: {style="margin-top: 1em; text-align: left;"}
例題のRコード（と、スライド作成のコード）は

<https://github.com/ito4303/AIC_PMSE>

においてあります。
:::