p 値関数のアニメーション - ジョンとヨーコのイマジン日記

前置き

この文書は統計的仮説検定とかを一度は学んだことがある人向けに書いている. 小ネタです.

オリジナリティのない教科書がそうであるように, 1標本の t 検定を例にとることをお許しください.

プロット; 信頼区間

t 検定では, 統計モデルとして正規分布を用い, 正規分布の平均パラメータが $\mu=\mu_0$ であることを帰無仮説とする.

正規分布を用いるのになぜ t 検定という名前になっているかというと, 帰無仮説のもとで検定統計量が従う分布が正規分布でなく t 分布だからである.

（t 検定の解説はいろんなところにあると思うのでここでは省略する.）

統計ソフト R では t.test 関数に t 検定の検定統計量とかその他もろもろが実装されている.

R の例題用データの一つに sleep というのがあって, これはなんか睡眠薬かなんかを飲んで睡眠時間の伸びを測ったデータである.

こんなふうだ.

sleep データのプロット

ヘルプによると "The 'group' variable name may be misleading about the data: They represent measurements on 10 persons, not in groups." です.

このデータに対して, 帰無仮説を $\mu = 0$ , 対立仮説を $\mu \neq 0$ とした t 検定（両側検定）を行うには例えばこんなふうにする.

sleep2 <- unstack(sleep, extra~group)
x_d <- sleep2[,2] - sleep2[,1]
res0 <- t.test(x_d)

単にデフォルトが $\mu = 0$ になっているだけで, 帰無仮説ではパラメータが 0 でなきゃいけないみたいなルールがあるわけではないので, ためしに帰無仮説 $\mu = \mu_0$ の検定を $\mu_0$ をいろいろ変えてかたっぱしからやってみる.

len <- 100
mu0 <- seq(xbar - 1.5, xbar + 1.5, length.out = len)
res <- lapply(mu0, function(mu)t.test(x_d, mu = mu))

有意水準を5%とし, 帰無仮説が棄却された $\mu_0$ を青, されなかった $\mu_0$ を赤で示す.

信頼区間のアニメーション

縦の点線の内側の範囲では帰無仮説が棄却されないことがわかる.

縦の点線はなにかというと平均 $\mu$ の95%信頼区間で, なんのことはないという感じですね.

プロット; p 値関数

棄却するしないの2値の情報だけだとアニメーションが単調なので, 縦軸を p 値（検定と検定統計量を与えたとき帰無仮説を棄却できる最小の有意水準）にしてみる.

$\mu_0$ をいろいろ変えて p 値を表示.

p 値のアニメーション

反対に p 値を動かして指定した p 値が得られる $\mu_0$ を表示.

p 値のアニメーション（逆）

2次元の情報しかないので実はアニメーションにする必要はなくて, gifアニメ上の点の軌跡は次の図で示すとおり.

p 値関数

これを称して p 値関数という.

95%信頼区間が有意水準5%で棄却されないパラメータの範囲を示すことを思いだすと, p 値関数を求めることは ( $1-\alpha$ ) 信頼区間を任意の $\alpha$ で求めることに相当する.

R のコード

まとめて貼るぞ.

library(gganimate)

prop_t <- function(mu, xbar = 0, se = 1, df = 1, 
                   alternative = "two.sided"){
  if (!alternative %in% c("two.sided", "less", "greater")) {
    stop("alternative must be either \"two.sided\", \"less\", or \"greater\".")
  }
  v <- (xbar-mu)/se
  p0 <- 0
  if(alternative == "two.sided"){
    p0 <- 2*pt(abs(v), df, lower.tail = FALSE)
  }else if(alternative == "less"){
    p0 <- pt(v, df)
  }else if(alternative == "greater"){
    p0 <- pt(v, df, lower.tail = FALSE)
  }
  return(p0)
}

ggplot(sleep, aes(x = group, y = extra, group = ID)) +
  geom_line() +
  geom_point() +
  theme_bw(24)

ggsave("./Desktop/sleep.png")

sleep2 <- unstack(sleep, extra~group)
x_d <- sleep2[,2] - sleep2[,1]
res0 <- t.test(x_d)
print(res0$statistic)

xbar <- mean(x_d)
se <- sd(x_d)/sqrt(nrow(sleep2))

len <- 100
mu0 <- seq(xbar - 1.5, xbar + 1.5, length.out = len)
res <- lapply(mu0, function(mu)t.test(x_d, mu = mu))
alpha <- seq(0, 0.5, length.out = len)

resdf <- data.frame(
  time_dummy = 1:len,
  p_dummy = c(rev(1:(len/2)),1:(len/2)),
  mu0 = mu0,
  stat = sapply(res, function(x)x$statistic),
  pval = sapply(res, function(x)x$p.value)
)

ggplot(resdf, aes(x = mu0, y = 0, colour = pval < 0.05)) +
  geom_point(shape = 4, size = 4, stroke = 2) +
  geom_vline(xintercept = res0$conf.int, linetype = 2) +
  labs(x = expression(mu[0]), y = "", colour = "reject") +
  theme_bw(24) +
  theme(axis.ticks.y = element_blank(), axis.text.y = element_blank()) +
  transition_time(time = time_dummy)

anim_save(file = "./Desktop/reject.gif")

ggplot(resdf, aes(x = mu0, y = pval))+
  geom_hline(aes(yintercept = pval), linetype = 2) +
  geom_vline(aes(xintercept = mu0), linetype = 2) +
  geom_point(size = 4)+
  labs(x = expression(mu[0]), y = "p-value")+
  theme_bw(24) +
  transition_time(time = time_dummy)

anim_save(file = "./Desktop/muve_mu.gif")

ggplot(resdf, aes(x = mu0, y = pval)) +
  geom_point(size = 4) +
  geom_hline(aes(yintercept = pval), linetype = 2) +
  geom_vline(xintercept = res0$conf.int, linetype = 2)+
  labs(x = expression(mu[0]), y = "p-value") +
  theme_bw(24) +
  transition_time(time = p_dummy)

anim_save(file = "./Desktop/muve_p.gif")

ggplot(resdf, aes(x = mu0))+
  stat_function(fun = prop_t,
                args = list(xbar = xbar, se = se, df = 9),
                size = 1, n = 501) +
  labs(x = expression(mu[0]), y = "p-value")+
  theme_bw(24)

ggsave("./Desktop/p_fun.png")