こんにちは，kur0ckyです．最近は辛い料理を作るのにハマっているのですが，そのせいで腸内環境が大変なことになってます．

さて，エフロン＆ヘイスティの「大規模計算時代の統計推論」に載っている解析をすべてRでやってみる，ということをやっていたのでした．前回記事はこちらです．

つづく第2章は「頻度派的な推論」というタイトルで，そもそも頻度主義とはなにか，というところから解説されています．本章は解析が少なく内容が軽くなってしまったので，私の知らなかった概念や，少しややこしい概念について補足的なポエムも書いてみました．誤りを見つけたときは石を投げてください．それではやっていきます．

はじめに
2.1 実践における頻度主義
- 不動統計量について
- 表2.1
2.2 頻度派的な最適性
- 統計的仮説検定について
- どのような検定が望ましいか
2.3 補足と詳細
おわりに

例のごとくtidyverseを読み込んでおきます．ggplot2の設定もしておきます．

library(tidyverse)
theme_set(theme_classic(base_family = "HiraKakuPro-W3"))

はじめに

頻度派主義の巨人について

ピアソン，フィッシャー，ネイマンが頻度主義の大家であることは知っていましたが，本書ではホテリングの名前もその一人として挙がっていました．ホテリングは古典的な異常検知で有名なホテリング理論ぐらいしか知らないのですが（無知），時間があればその辺りも調査してみたいです．言われてみればたしかに，データの背後に多次元正規分布を仮定し，分布から遠い（非常にまれな事象）を異常値とみなす，というのはなんだか頻度主義っぽいかもしれません．

腎臓病のデータについて

本章ではまずはじめに，腎臓病検査のgfrという数値を例にあげ，ある統計量の確率的性質をみることの有用性を示しています．

さっそく公式HPからデータをダウンロードしましょう．変数は1つしかなく，レコードが改行で区切られているtxtファイルになっています． readr::read_csv()などでも読むことができますが，ここではbase::scan()を使ってみましょう． scan()はファイルからベクトルを読み込む関数です． read_lines()でも読むことができる（ただしcharacter型になる）ので，そちらの例も合わせてしめしておきます．ちなみに，readr::read_lines_raw()を用いると，エンコーディング前の未加工の状態でテキストを読み込むことができます．エンコーディングが予めわからない場合などに用いると有用でしょう．

gfr <- scan("https://web.stanford.edu/~hastie/CASI_files/DATA/gfr.txt",
            what = double())
# read_lines("https://web.stanford.edu/~hastie/CASI_files/DATA/gfr.txt") %>% 
#   as.integer()

続いて，平均値とその標準誤差を算出します．平均値の標準誤差は，各要素が同一の確率分布から独立に得られていると考えると，以下のように算出することができます．平均値は54.2654，標準誤差は0.9445844となり，は本書の値とも一致しますね．本書ではこれらの数値をもとに，小数点以下の数値が信用に足らないことを指摘しています．

簡単なのでライブラリなどは使うまでもないですが，強いて使うならば，rapportools::se.mean()などがあります．

# 平均値の算出
mean(gfr)
# 標準誤差の算出
sqrt(var(gfr) / length(gfr))
# パッケージを使う例
rapportools::se.mean(gfr)

分布を確認するためにヒストグラムを描きましょう． ggplot2で書くために一旦tibble（data.frame）に直します．本書の図にできるだけ近づけるため，binsも色々ためしてみましたが，結局完璧には合わせられず断念しました笑．

df_gfr <- tibble(gfr)
df_gfr %>% 
  ggplot(aes(gfr))+
  geom_histogram(bins = 30, colour = "black", fill = "green")+
  labs(y = "度数",
       x = "gfr")

f:id:kur0cky:20200823003751p:plain

頻度主義的推論の手続きについて

頻度主義でよく行われる推論について，感覚的な紹介がされています．まず，以下のような数学的対象が導入されます．

$F$ ：データが従う確率分布
$X$ ：理論的標本（確率変数）． $F \rightarrow {\bf X}, {\bf X} = (X _ 1, \dots, X _ n)$
$x$ ：標本の実現値． ${\bf x} = (x _ 1, \dots, x _ n)$
$\theta$ ：確率分布 $F$ の母数（もしくは性質）
$\hat\Theta = t({\bf X})$ ： $\theta$ の理論的な標本．いわゆる推定量．確率的な性質をもつ
$\hat\theta = t(\bf x)$ ：推定量に実際に観測された標本を入れたもの．推定値．単一の値

似たような記号がたくさん出てくるため，腑に落ちるまで少し時間がかかるかもしれません．また， $\hat\Theta$ や $X$ を考える意味も分かりにくいかもしれません．

ただ，上で太字で示したように， $\hat\theta$ や ${\bf x}$ はあくまでも実現した1つの値です．そして，「なんらかの計算アルゴリズム $t$ 」に対し，確率的な性質を考えることこそが重要であるため（確率のことばで議論するため），確率変数や理論的標本を考えるのだと私は整理しています．

2.1 実践における頻度主義

本節では，頻度主義でしばしば行われるいくつかの工夫を通して，「頻度主義とは何か」がフワっと述べられています．

上で触れたように，頻度主義では，ある計算アルゴリズム（たとえば平均など）に対して，その確率的性質を導き出します．このとき，データが従う確率分布 $F$ を知りもしないのに， $F \rightarrow {\bf X}$ を用いた $\hat\Theta = t({\bf X})$ の性質を議論しなくてはなりません．そこで，本書では次の5つの工夫を紹介しています

差し込み法 (plug-in principle)
テイラー級数近似（デルタ法; delta method）
パラメトリック分布族と最尤理論
模擬実験とブートストラップ法
不動統計量 (pivotal statistic)

不動統計量について

上の4つは有名ですし，後の章で詳しく触れられるのでここでは割愛します．しかし，これは私が無知なだけですが，不動統計量というのは初めて聞きました．本書では，不動統計量について以下のように説明し，その例としてスチューデントの2標本 $t$ 検定を紹介しています．

不動統計量 $\hat{\theta} = t(x)$ とは，背後にある確率分布 $F$ に依存しない分布を持つ統計量のことである．

$t$ 検定は標本が正規分布に従うことを仮定するので，背後にある確率分布 $F$ に依存しているじゃないか！！と思ったのですが，読み進めて次のように解釈しました．

2つの標本を $\bf x_1, x_2$ としたとき，検定統計量は次のようになります．

$\displaystyle \begin{align} t = \frac{\bar x_2 - \bar x_1}{\hat{\rm sd}}, \hat{\rm sd} = \hat{\sigma}\left(\frac{1}{n_1}+\frac{1}{n_2}\right)^{1/2} \end{align}$

この統計量は仮定した正規分布の母数 $\sigma$ に依存せず，自由度 $n_1 + n_2$ の $t$ 分布に従います．以上から推察するに，不動統計量とは，「特定の分布（族）を仮定するかさておき，その母数に依存しないように設計された統計量」ということでしょうか．

どなたか詳しい方がいたら教えて下さい．

表2.1

ブートストラップ法の使用例として，gfrデータに対する3つの位置推定値と，その標準誤差を推定するタスクを紹介しています．「位置」推定値には平均値，中央値，ウィンザー化平均値を使用しています．平均であればその標準誤差は陽に求まりますが，他2つに関しては難しいためブートストラップ法で推定することになります．

まずは点推定値を求めます．ウィンザー化平均の条件分岐にはみんな大好きcase_when()を使います．

# 予め第1，第3四分位点を求めておく
q1 <- quantile(gfr, 0.25)
q3 <- quantile(gfr, 0.75)
estimated_value <- df_gfr %>% 
  mutate(winsorized_gfr = case_when(gfr < q1 ~ q1,
                                    gfr > q3 ~ q3,
                                    TRUE ~ gfr)) %>% 
  summarise(mean = mean(gfr),
            median = median(gfr),
            winsorized_mean = mean(winsorized_gfr)) %>% 
  gather(var, estimated_value)

print(estimated_value)

# A tibble: 3 x 2
  var             estimated_value
  <chr>                     <dbl>
1 mean                       54.3
2 median                     52  
3 winsorized_mean            52.8

続いて，中央値とウィンザー化平均の標準誤差をブートストラップ法で求めます．ブートストラップ標準誤差の値が本書と若干ずれてますが，ブートストラップの反復回数が少ないということにしておきましょう．

boot_size = 1000
size = length(gfr)
set.seed(1)
se <- df_gfr %>% 
  sample_n(size = n()*boot_size,
           replace = T) %>% 
  mutate(boot_index = rep(1:boot_size, each = size)) %>% 
  group_by(boot_index) %>% 
  mutate(q1 = quantile(gfr, 0.25),
         q3 = quantile(gfr, 0.75),
        winsorized_gfr = case_when(gfr < q1 ~ q1,
                                    gfr > q3 ~ q3,
                                    TRUE ~ gfr)) %>% 
  summarise(median = median(gfr),
            winsorized_mean = mean(winsorized_gfr)) %>% 
  summarise(median = sd(median),
            winsorized_mean = sd(winsorized_mean)) %>% 
  mutate(mean = sqrt(var(gfr) / size)) %>% 
  gather(var, se)

print(se)

# A tibble: 3 x 2
  var                se
  <chr>           <dbl>
1 median          0.874
2 winsorized_mean 0.905
3 mean            0.945

最後に，点推定値と合わせておしまいです．表示順をslice()で調整したり，mutate_if()とround()を組み合わせて数値型のカラムだけ小数点を丸めたりしています．こういうときに***_all()系や***_at()，***_if()は便利ですね．

se %>% 
  right_join(estimated_value, ., by = "var") %>% 
  slice(c(3, 2, 1)) %>% 
  mutate_if(is.numeric, round, 2) %>% 
  as.data.frame() %>%
  column_to_rownames("var") %>% 
  rename(推定値 = estimated_value,
            標準誤差 = se) %>% 
  knitr::kable()

	推定値	標準誤差
mean	54.27	0.94
winsorized_mean	52.81	0.91
median	52.00	0.87

綺麗な表ができあがりました．それにしても，ウィンザー化平均というものは初めて知りました（無知）．調べてみるとロバスト統計ではしばしば用いられているようです．知らないことだらけです．

2.2 頻度派的な最適性

本節では，頻度主義のメリットを「考えるべきは確率モデル $F$ とアルゴリズム $t$ の選択だけで済むこと」としています．また，より良い $t$ の選択手法の例としてフィッシャー情報量限界 (Fisher information bound) とネイマン=ピアソンの補題 (Neyman-Pearson lemma) を取り上げています．

統計的仮説検定について

これは私の感想ですが，ネイマン=ピアソンの補題をこのページ数で取り上げるのは少し無理がある気がしました．本書は，仮説検定の枠組みをよく知っていることが前提にされている気がします．そこで，理解のために必要だと思われるいくつかの概念について，触りだけ紹介しようと思います．より詳しく知りたい方は宿久先生の特集論文などを読むと良いと思います．

www.jstage.jst.go.jp

ここでは次の概念を紹介します．

仮説検定
第一種の過誤
第二種の過誤
ネイマン=ピアソンの基準
単純仮説
複合仮説
検出力
最強力検定
ネイマン=ピアソンの補題

まず仮説検定についてですが，その手続きをざっくりと書くと，以下のようになります．

観測したデータが帰無仮説のモデルのもとで起こる確率を計算する．
設定した有意水準未満でしか起こらない事象を観測したとき，
帰無仮説の仮定のもとでは，起こっている事象はレアすぎると考える．
つまり，帰無仮説が成り立っていないのでは，と判断する．

もちろん，帰無仮説を否定することはできても，支持はできないことには注意が必要です．ここで，有意水準を閾値として判断を行っているため，判断には一定の誤りが伴うことになります．この誤りを次の2つに分類します．

第一種の過誤：帰無仮説が正しいのに，帰無仮説を棄却してしまう誤り
第二種の過誤：逆に，対立仮説が正しいのに，帰無仮説を棄却しない誤り

一般に，第一種の過誤と第二種の過誤はトレードオフの関係にあります．どのように検定を設計するのがよいでしょうか．

どのような検定が望ましいか

ネイマン=ピアソン基準では，どこまで第一種の過誤を許容できるか，つまり有意水準 $\alpha$ を先に定め，そのもとで第二種の過誤 $\beta$ を最小化します．また，対立仮説が正しいときに，帰無仮説を棄却する確率を検出力といいます．つまり，検出力は $1-\beta$ です．まとめると，ネイマン=ピアソン基準を用いるということは，ある $\alpha$ のもとで $1-\beta$ がより大きい検定を「良い検定」とみなすことです．