BUSINESS DATA SCIENCEという本を最近読んでいるので内容を自分なりにまとめる

Business Data Science: Combining Machine Learning and Economics to Optimize, Automate, and Accelerate Business Decisions

作者:Taddy, Matt
発売日: 2019/08/21
メディア: ハードカバー

最近「ビジネスデータサイエンスの教科書」という邦題で日本語版も出たみたいです

ビジネスデータサイエンスの教科書

作者:マット・タディ
発売日: 2020/07/22
メディア: 単行本

この章はなにか

ビジネスをはじめとした実世界はノイズを多く含んでいるのでノイズを考慮したモデルを作成する必要がある。
この章では、ノイズに対して統計的手法で不確実性を扱うことによって対処する方法を提供している。

データなどは作者のgitにある。

Uncertaintyに対するブートストラップ法に関しては過去記事

knknkn.hatenablog.com

この節はなにか

推定した結果がどの程度信用していいか判断するために、検定をおこなう。その検定をおこなう際の注意点を書いている。

前回の記事と同様のデータセットを使う。

# A tibble: 10,000 x 7
      id anychildren broadband hispanic race  region spend
   <dbl>       <dbl>     <dbl>    <dbl> <chr> <chr>  <dbl>
 1     1           0         1        0 white MW       424
 2     2           1         1        0 white MW      2335
 3     3           1         1        0 white MW       279
 4     4           0         1        0 white MW       829
 5     5           0         1        0 white S        221
 6     6           0         1        0 white MW      2305
 7     7           1         1        0 white S         18
 8     8           1         1        0 asian MW      5609
 9     9           0         1        0 white W       2313
10    10           1         1        0 white NE       185
# … with 9,990 more rows

検定

今回は、以下の重回帰をおこなった結果の回帰係数について考える。

spendy = glm(log(spend) ~ .-id, data = data)
round(summary(spendy)$coef, 2) 
# Estimate Std. Error t value Pr(>|t|)
# (Intercept)     5.86       0.16   36.34     0.00
# anychildren     0.09       0.03    2.54     0.01
# broadband       0.52       0.04   11.93     0.00
# hispanic       -0.18       0.04   -4.30     0.00
# raceblack      -0.25       0.18   -1.41     0.16
# raceother      -0.41       0.31   -1.32     0.19
# racewhite      -0.21       0.15   -1.36     0.17
# regionNE        0.26       0.05    4.98     0.00
# regionS         0.01       0.04    0.13     0.90
# regionW         0.18       0.05    3.47     0.00

この結果の Pr(>|t|)が各係数が起こりうる確率を表している。

より厳密には、各係数の t value を求め、その t value が起こりうる確率をPr(>|t|)で表している。

これは、以下のような、x軸が t value y軸が確立となる確率密度関数であるt分布を考えている。

f:id:chito_ng:20191231155617p:plain

t value は 期待値の差/sqrt(分散/n) で求めることができる。ここでは係数=0が期待値となる場合と、その係数が期待値となる場合で有意に異なるかを t value で表し、それを確率で表しているのが Pr(>|t|) となっている。

そして、このPr(>|t|) がある一定の値以下(一般的には0.05)のときその係数を有意としている。

軽い説明になりましたが、詳細な説明は以下がわかりやすかったです。

logics-of-blue.com

なお、書籍での説明はt分布ではなくz分布で説明をしている。おそらくデータ量が多いので係数の期待値の分布を中心極限定理より正規分布に仮定しても問題なさそうなためか？

bellcurve.jp

多重検定

結果に対して同時に複数の検定をおこなうのは一般的に注意が必要となる。

例えば、100個の変数を用いた重回帰分析を考える。このとき、真の結果として100個の変数のうち5つが目的変数と関係があるとする。つまり、真のモデルでは5つの変数の係数は0以外で、他の95個は0(目的変数に無関係）となる。

5%水準での検定は言い換えると、5%の確率で「本当は有意ではない」のに「有意」だと判断してしまう。
そのため、目的変数に無関係の95個の変数は95 * 0.05 = 4.75となることから、95個中5つが誤って有意と判定されてしまう。

100個の変数を入れることはあまりないかもしれないが、10個の変数の場合を考えても10 * 0.05 = 0.5となり、変数のうち1つ誤って有意とする可能性が50%となることから、50%の確率で誤ったモデルを作成することになる(「全ての検定結果があっている」可能性が50%)。
そして、当たり前だが誤った変数を入れたモデルの推定は誤っているのでこのことは非常に重要な問題となる。

Benjamini-Hochberg(BH)法による FDRコントロール

このとき、誤って有意としてしまう（真の帰無仮説を棄却してしまう）割合を一定に抑えたい。つまり、 $\frac{誤って有意と判定された数(FP) } {有意と判定された数(TP+FP)}$ で定義されるFalse Discovery Proportion(FDP)が一定になるように新たな検定水準を考える。

例えば、先程同様に100変数の重回帰分析のうち、真に有意な5変数と誤って有意とした5変数の計10変数を考える。このとき、FDR = 0.2とすることは誤って有意とする変数の数を10 * 0.2 = 2個未満になるように有意な変数を再考することとなる。

しかし、FDPの分子の「誤って有意と判定された数」は何が真かわからないので知ることはできない。
一方で、FDPの期待値はBenjamini-Hochberg(BH)法を用いることで利用することができる。つまり、BH法を用いることで「有意判定のうち誤って有意となる割合(FDP)の期待値をq未満に抑える」といったことをおこなうことができる。なお、FDPの期待値をFalse Discovery Rate(FDR)と呼ぶ。

BH法の定義は以下。
1. 変数の数をN、各係数のp値を昇順に [tex: p\{1}, p\{2}, p\{3}...p\{N}] とし、FDRをqとする
2. kを1から順に増加させたときの [tex: p\{k} \leq q \frac{k}{N}] の最大となる [tex: p\k] をカットオフ値pとする
3. $p\_k$ 以下となる全p値を全て有意とし、それ以外を全て有意ではないとして再定義をおこなう(新たな有意水準をpとする)

なお、このときのFDRはq以下となるし、言い換えるとFDPは平均的にはq以下となる。

BH法の実例

先ほどと同様のデータで考える。

spendy = glm(log(spend) ~ .-id, data = data)
round(summary(spendy)$coef, 2) 
# Estimate Std. Error t value Pr(>|t|)
# (Intercept)     5.86       0.16   36.34     0.00
# anychildren     0.09       0.03    2.54     0.01
# broadband       0.52       0.04   11.93     0.00
# hispanic       -0.18       0.04   -4.30     0.00
# raceblack      -0.25       0.18   -1.41     0.16
# raceother      -0.41       0.31   -1.32     0.19
# racewhite      -0.21       0.15   -1.36     0.17
# regionNE        0.26       0.05    4.98     0.00
# regionS         0.01       0.04    0.13     0.90
# regionW         0.18       0.05    3.47     0.00

p値のみ抜き取ると以下。

pval = summary(spendy)$coef[-1, 'Pr(>|t|)']
# anychildren    broadband     hispanic    raceblack    raceother    racewhite     regionNE 
# 1.122431e-02 1.323047e-32 1.696801e-05 1.594329e-01 1.878053e-01 1.740374e-01 6.440516e-07 
# regionS      regionW 
# 8.977425e-01 5.322114e-04

これを昇順に並べると以下。

# A tibble: 9 x 1
      pval
     <dbl>
1 1.32e-32
2 6.44e- 7
3 1.70e- 5
4 5.32e- 4
5 1.12e- 2
6 1.59e- 1
7 1.74e- 1
8 1.88e- 1
9 8.98e- 1

FDRを0.1としたときの $p_{k} \leq q \frac{k}{N}$ を考えると、
1. $p_1 = 1.32e-32 \lt 0.1 \frac{1}{9}$
2. $p_2 = 6.44e- 7 \lt 0.1 \frac{2}{9}$
3. $p_3 = 1.70e- 5 \lt 0.1 \frac{3}{9}$
4. $p_4 = 5.32e- 4 \lt 0.1 \frac{4}{9}$
5. $p_5 = 1.12e- 2 \lt 0.1 \frac{5}{9}$
6. $p_6 = 1.59e- 1 \gt 0.1 \frac{6}{9}$

以上より、k = 5が最大となり、新たな有意水準は $p_{5} = 1.12e- 2$ となった。つまり、9変数のうち、5変数が有意となり4変数が有意とはいえない。このとき、誤って有意としている割合の期待値(FDR)は0.1以下となっている。

なお、これはp値と昇順rankの散布図に $\frac{q}{N}$ となる傾きの直線を引いたときその直線以下を有意としている、とも捉えることができる。

pval = summary(spendy)$coef[-1, 'Pr(>|t|)']

tibble(pval) %>%
  arrange(pval) %>% 
  mutate(rank = row_number()) %>% 
  ggplot(aes(rank,pval)) +
  geom_point() + 
  geom_abline(slope = 0.1/9)

f:id:chito_ng:20200101034600p:plain

BH法が使えない場合

BH法はp値が一様分布から独立に生成されているという仮説のもとなりたつ。そして、例えば変数をひとつ外すだけで多くの場合でp値が変わるように、この仮定は満たされないことが多い。

この欠点を克服、つまりp値同士が独立していなくても問題がない手法として、カットオフ値を求める際にkに応じたweightを定めるBenjamini-Yekutieli法(BY法)が存在する*1。

stats.biopapyrus.jp

また、ベイズ的なアプローチを取りMCMCによって事後分布を推定するという方法もある。ただし、この手法は事前分布などに仮定を多く置いているので誤った推定をおこなうことが往々にして起こる。そのため、ブートストラップ法のような仮定が必要ないノンパラメトリックな手法が使えない場合(超多変量、n数が少ないなど)を除くとノンパラメトリックな手法を用いる法がロバストなので推奨されている。