- 尤度比検定
- 逸脱度の差に注目
- ネストしているモデルたちを比較できる
- パラメトリック
- ノンパラメトリック
5.1 統計学的な検定のわくぐみ¶
Neyman-Pearsonの検定のわくぐみ（AICによるモデル選択との比較，図5.1）
- 使用するデータを決める
- 適切な統計モデルの設定，パラメータの最尤推定
- 帰無仮説と対立仮設
- 帰無仮説棄却の危険率を評価
- 検定統計量を決める
- 帰無仮説が成り立つと仮定した場合の検定統計量の値が発生しても珍しくない範囲を定める（95%なら有意水準5%）
- 帰無仮説が棄却できるかを判定
5.2 尤度比検定の例題: 逸脱度の差を調べる¶
- 3章データ（種子数）
- 平均が $\lambda = \exp(\beta_1 + \beta_2 x_i)$ のポアソン分布のGLM
- 一定モデルとxモデルを比較，一定モデルが棄却できるかどうかを調べる
In [ ]:
if (interactive()) {
(function(x){
notinstalled <- x[!(x %in% installed.packages()[, "Package"])]
if (length(notinstalled) > 0) {
cat(sprintf("Install %d packages: %s \n", length(notinstalled), notinstalled))
install.packages(notinstalled)
} else {
cat("Already installed\n")
installed.packages()[, "Version"][x]
}
})(c("pipeR", "dplyr", "tidyr", "ggplot2", "readr", "knitr", "devtools"))
}
In [1]:
sapply(c("pipeR", "dplyr", "tidyr", "ggplot2", "readr"), require, character.only = TRUE)
In [2]:
d <- read_csv("data/chap03/data3a.csv")
str(d)
一定モデル
In [3]:
fit.null <- glm(y ~ 1, data = d, family = "poisson")
fit.null
In [4]:
logLik(fit.null)
xモデル
In [5]:
fit.x <- glm(y ~ x, data = d, family = "poisson")
fit.x
In [6]:
logLik(fit.x)
フルモデルの最大対数尤度，逸脱度，AIC
In [7]:
fit.full <- list()
fit.full$maxLogLik <- dpois(d$y, lambda = d$y) %>>% log() %>>% sum()
fit.full$dev <- -2 * fit.full$maxLogLik
fit.full$aic <- -2 * (fit.full$maxLogLik - length(d$y))
fit.full
In [8]:
options(repr.plot.width = 4, repr.plot.height = 4)
In [9]:
ggplot(aes(x = x, y = y), data = d) +
geom_point() +
stat_smooth(formula = y ~ 1, method="glm", method.args = list(family = "poisson"), se = FALSE, size = 0.5, linetype = "dashed") +
stat_smooth(formula = y ~ x, method="glm", method.args = list(family = "poisson"), se = FALSE, size = 0.5) +
xlab(expression("Body size"~italic(x)[i])) +
ylab(expression("Number of seeds"~italic(y)[i]))
In [10]:
data.frame(
"model" = c("Null", "x", "Full"),
"k" = c(1, 2, 100),
"max.log.L" = c(logLik(fit.null), logLik(fit.x), fit.full$maxLogLik),
"deviance" = c(fit.null$deviance, fit.x$deviance, 0) + fit.full$dev,
"resd.dev." = c(fit.null$deviance, fit.x$deviance, 0),
"aic" = c(fit.null$aic, fit.x$aic, fit.full$aic)
) %>>% mutate_each(funs(round(., 1)), max.log.L:aic) %>>% knitr::kable()
尤度比:
$$\frac{L_1^\ast}{L_2^\ast} = \frac{\text{Max. likelyhood of Null model}}{\text{Max. likelyhood of x model}}$$尤度比検定の検定統計量: 逸脱度の差（尤度比の対数 × -2）
$$\Delta D_{1, 2} = -2 \times (\log L_1^\ast - \log L_2^\ast)$$例題の場合だと，約4.5の改善
In [11]:
(fit.null$deviance - fit.x$deviance)
5.3 2種類の過誤と統計学的な検定の非対称性¶
統計学的な検定
- 帰無仮説: Null model
- 対立仮設: x model
2種類の過誤
- Type I error: 真のモデルは Null model なのに 観測データが Null model からではなく x model から生じたと判断する（帰無仮説を棄却）
- Type II error: 真のモデルはNull model ではないのに，観測データが Null model から生じたと判断する（帰無仮説を棄却できない）
Type I error を検討
- 帰無仮説が正しいとする（$\hat{\beta}_1 = 2.06$ を得る）
- $\hat{\beta}_1 = 2.06$ のモデルからいくつもデータを生成して，$\Delta D_{1,2}$ の分布を得る
- $\Delta D_{1,2} \ge 4.5$ なる確率 $P$ を計算する
- $P$ を検討
5.4 帰無仮説を棄却するための有意水準¶
有意水準$\alpha$を設定しておいて，
- $P \ge \alpha$: 帰無仮説を棄却できない
- $P \lt \alpha$: 帰無仮説を棄却
とする（$\alpha = 0.05$ が多い）．
以下では2種類の方法で$P$値を計算する
5.4.1 方法（1）汎用性のあるパラメトリックブートストラップ法¶
- 平均
mean(d$y)のポアソン乱数を生成- Null model で推定した平均種子数: $\exp(\beta_1)$
- 生成した乱数から $\Delta D_{1,2}$ を計算
In [12]:
readLines("Data/chap05/pbnew.R", encoding = "UTF-8") %>>% cat(sep="\n")
In [13]:
readLines("Data/chap05/pb.R", encoding = "UTF-8") %>>% cat(sep="\n")
In [14]:
source("Data/chap05/pb.R")
# or
# source("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/lrtest/pb.R")
pb
In [15]:
dd12 <- pb(d, n.bootstrap = 1000)
In [30]:
dd12.x <- pb(d, n.bootstrap = 10000)
In [16]:
summary(dd12)
In [17]:
options(repr.plot.width = 4, repr.plot.height = 4)
$10^3$回
In [18]:
data.frame(x = dd12) %>>%
ggplot(aes(x = x)) +
geom_histogram(bins = 100, colour = "black", fill = "white", size = 0.3) +
scale_x_continuous(limits = c(0, 20), breaks = seq(0, 20, 5)) +
theme_bw() +
geom_vline(xintercept = 4.5, size = 0.5, linetype = "dotted") +
xlab(expression(Delta~italic(D)[`1,2`]))
$10^4$回
In [33]:
data.frame(x = dd12.x) %>>%
ggplot(aes(x = x)) +
geom_histogram(bins = 100, colour = "black", fill = "white", size = 0.3) +
scale_x_continuous(limits = c(0, 20), breaks = seq(0, 20, 5)) +
theme_bw() +
geom_vline(xintercept = 4.5, size = 0.5, linetype = "dotted") +
xlab(expression(Delta~italic(D)[`1,2`]))
$\Delta D_{1,2} \ge 4.5$ となった個数
In [19]:
sum(dd12 >= 4.5)
In [31]:
sum(dd12.x >= 4.5)
In [20]:
as.numeric(c(TRUE, FALSE))
$P$ 値は
In [21]:
sum(dd12 >= 4.5) / 1000
$P = 0.05$ となる$\Delta D_{1,2}$は
In [22]:
quantile(dd12, 0.95)
$\Delta D_{1,2} = 4.5$ の $P$値は有意水準 0.05 より小さいので，帰無仮説を棄却
library(boot) を使うと，
In [23]:
library(boot)
In [24]:
b <- boot(
d,
function(data){
fit1 <- glm(y ~ 1, data = data, family = poisson)
fit2 <- glm(y ~ x, data = data, family = poisson)
fit1$deviance - fit2$deviance
},
R = 1000,
sim = "parametric",
ran.gen = function(data, mle){
out <- data
out$y <- rpois(nrow(out), lambda = mle)
out
},
mle = mean(d$y)
)
b
In [25]:
str(b)
In [26]:
summary(b$t[,1])
In [27]:
sum(b$t >= 4.5)
In [28]:
sum(b$t >= 4.5) / 1000
In [29]:
quantile(b$t, 0.95)
5.4.2 方法（2）$\chi^2$分布を使った近似計算法¶
$\Delta D_{1,2}$ の確率分布は自由度1の$\chi^2$分布で近似できる場合がある
- サンプルサイズが大きい場合に有効な近似計算
In [32]:
anova(fit.null, fit.x, test = "Chisq") %>>% print()
$P = 0.034$ で帰無仮説は棄却される
5.5 「帰無仮説を棄却できない」は「差がない」ではない¶
5.6 検定とモデル選択，そして推定された統計モデルの解釈¶
5.7 まとめ¶
- 尤度比検定: 逸脱度の差を調べる
- Neyman-Pearson のわくぐみ
- Type I error に注目
- 検出力（1 - Type II error の確率）を話題にすることもある
- 有意水準の決め方に根拠はない
- 帰無仮説が棄却できない場合は何も言えない
- Type I error に注目
- 検定やモデル選択の結果だけでなく，推定されたモデルが現象をどのように予測するのかも確認すべき
In [34]:
devtools::session_info()