【読後メモ】データ解析のための統計モデリング入門（その1）

よくある統計学の入門書では、確率の基礎や確率分布の紹介が多くの紙面を割いていることがあります。しかし本書では、一般化線形モデルに焦点を絞ることで、確率・統計の基礎についてはあまり触れられてはいません。これだけ書くと難しいように感じるかもしれませんが、逆に焦点を絞ることで本書を通して一貫した話題になっており、一般化線形モデルについて理解しやすいのではないかと思います。また、例題をベースにしているのも理解が進みやすいポイントかなと思います。

入門書と謳っておきながら、中身はだいぶハードという書籍もありますが、こちらは本当に入門書として良いのかなと思います。

本記事は、この書籍の各章（一部飛ばしていますが）でどのような手法が書かれているのか、どのようなことが書かれているのかといったメモです。長くなったので、数回に分けてまとめていきます。「その1」として一般化線形モデルまで（書籍の1,2,3,6章）をまとめます。

なお、私の別のブログで具体的事例（チョコボールのエンゼル推定）を題材にして統計モデリングの基礎についてまとめた記事があります。この書籍の内容はこの記事でもある程度カバーしていますので、そちらも参照いただけたらと思います。

chocolate-ball.hatenablog.com

【トップに戻る】

第1章データを理解するために統計モデルを作る

第1章はIntroductionとして、統計モデルとは何か、なぜモデリングするのかについて述べられています。

統計モデル

1.1節では統計モデルというものについて概説しています。

私はこの章（1章全体）はすごく重要だなと思っています。企業活動に於いて、データを集めて統計量を算出するみたいなことは日常的に行われる活動と思いますが、データを集めて何がしたいのか、統計量（平均など）を出して何がしたいのかを考えずにとりあえずデータを集めるのが重要と思っているケースを私自身何度か見てきました。データを集めて解析することで、その背景の仕組み*1を理解して活用することが目的なはずです。データ背景の仕組みを理解するために、確率分布などの統計学の手法を利用することで、活用がしやすくなります。

このあたりの考えは、当ブログの以下の記事にも想いを書いています。 learning-with-machine.hatenablog.com

ブラックボックスな統計解析の悪夢

1.2節では、統計モデリングの中身を理解しないまま手元のデータに当てはめてしまうことへの警鐘を鳴らしています。 PythonやRには有名な統計解析ライブラリがあり、何らかの結果が欲しければライブラリを使えばとりあえずの結果は出せてしまいます。

しかしデータ分析の結果を「活用」する場面においては、中身を理解しないまま統計モデリングのライブラリの出す結果がわかったところで、そこから先に進まないのではないかなと思います。この書籍では「線形」モデルを題材にしているので、全てのデータ解析シーンに十分な訳ではありませんが、まずはこの書籍の内容を抑えることがスタート地点なのかなと思います。

ライブラリに頼る解析は脱して、正しく統計モデリングできるように学んでいきたいですね。自戒を込めて。

本書籍では、「線形モデル」→「一般化線形モデル」→「一般化線形混合モデル」→「階層ベイズモデル」の順に徐々にモデルを複雑にして現実の様々な問題に適用できるように説明されています。データ解析の基本は、「シンプルなモデルから考える」です。この記事でも、その流れで見ていこうと思います。

【トップに戻る】

第2章確率分布と統計モデルの最尤推定

第2章では、植物の種子数を題材にして、統計モデリングの最も基礎となる単一パラメータの推論を解説しています。

データを解析するには、以下のような流れを踏むのが主流と思います。

データを確認
データを当てはめる確率分布を仮定
確率分布のパラメータを推定
当てはまりを確認

そして、改良の余地があれば最初に戻って繰り返します。ということで、最もシンプルなモデルから初めて見ます。

確率分布

第2章では、この流れに従って種子数がポアソン分布でモデル化できるとして、そのポアソン分布のパラメータをデータに基づいて推定しています（最尤推定）。ここで、「ポアソン分布」が唐突に出てくるのですが、2章の最後に確率分布の選び方の指針が示されています。曰く、次の3つの点に注目します。

- 説明したい量は離散か連続か？
- 説明したい量の範囲は？
- 説明したい量の標本平均と標本分散の関係は？

これらの点に従ってこの書籍では以下の確率分布が登場します。

- ポアソン分布：データが離散値、ゼロ以上の範囲、上限特になし、平均≒分散
- 二項分布：データが離散値、ゼロ以上で有限の範囲、分散は平均の関数
- 正規分布：データが連続値、範囲が[-∞, +∞]、分散と平均は無関係
- ガンマ分布：データが連続値、範囲が[0, +∞]、分散は平均の関数

これらの分布はメジャーでよく使われます。まずはこの4つの確率分布に従うとしてモデルを組み立てるのが良いのかなと思います。

最尤推定

上記の指針に基づいてモデルを設定し、そのパラメータを推定してみます。第2章では「最尤推定」という最も基本的な手法を使って推定してみます。

最尤推定では、まず尤度 $L(\theta)$ という「データの当てはまりの良さ」を定義し、その尤度 $L(\theta)$ が最大になるパラメータ $\theta$ を推定値 $\hat{\theta}$ とします。

まず、種子数のモデルとしてポアソン分布を設定します。

${ \displaystyle p(y|\lambda) = \frac{\lambda^{y}\exp{(-\lambda)}}{y!} }$

求めたいパラメータは $\lambda$ です。また、yはデータです( $y={y_1, y_2, \cdots y_n}$ )。この $\lambda$ の関数である尤度 $L(\lambda)$ は以下のように定義されます。

${ \displaystyle L(\lambda) = p(y_1 | \lambda)p(y_2 | \lambda)\cdots p(y_n | \lambda) \\ = \prod_{i}{p(y_i | \lambda)} = \prod{\frac{\lambda^{y_i}\exp{-\lambda}}{y_i!}} }$

これを対数変換し*2、パラメータ $\lambda$ で偏微分して0になる $\lambda$ が推定値 $\hat{\theta}$ です。

計算は省略しますが、結果は標本平均になります(簡単なので、手で計算してみましょう)。

${ \displaystyle \hat{\lambda} = \frac{1}{n}\sum_{i}{y_i} }$

ということで、長々書いてきましたが、結局のところ平均を求めれば良いのです。ポアソン分布、二項分布、正規分布は最尤推定量が標本平均になります。これは上記の計算手順に従えば容易に求められます。他の分布はどうなのか計算してないのでわからないのですが、基本を押さえておくのは重要と思いますので。

書籍では、このようにして得られた推定値 $\hat{\theta}$ を使ってシミュレーションすることで、データとモデルの当てはまりを確認できるということが書かれています。また、次に得られるデータの予測も推定値を使うことで可能になります。

【トップに戻る】

第3章一般化線形モデル（GLM）

第3章では一般化線形モデルが説明されています。この章で扱う例題は、第2章と同様に種子数を予測する問題です。第2章ではどの個体も共通でモデルを仮定しましたが、ここでは種子数は個体の特徴に依って変わるであろうという考えをモデルに組み込みます。一段階モデルが複雑になりましたね。

線形モデル

「"一般化"線形モデル」の前に、線形モデルから見ていきます。

線形モデルについては、この書籍ではあまり触れられていないのですが、以下の式のように期待値が説明変数( $\bf{x}$ )と係数( $\beta$ )の線型結合で表現されるモデルです。

${ \displaystyle y_i = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \epsilon_i = \bf{\beta}\bf{x} + \epsilon_i }$

ここで、 $\beta=[\beta_0, \beta_1, \beta_2, \cdots$ ]、 $x=[0, x_1, x_2, \cdots$ ]、 $\epsilon$ は正規ノイズ $(N(0, \sigma^{2}))$ です。

つまり、ノイズは正規分布に従うと仮定し、線型結合（直線のようなもの）を平均としたモデルとなります。なお、説明変数( $x$ )は、 $x^{2}, x^{3}$ のように多項式を使っても線型モデルになります（多項式回帰）。

このモデルを図で示すと以下のようになります（手書きでごめんなさい）。

f:id:hippy-hikky:20190219233312p:plain:w350

書籍でも述べられているように、直線と正規ノイズを仮定しているため、種子数のようなカウントデータのモデルには合わないことがわかります。所得のモデルなども同様ですね。

一般化線形モデル

線形モデルには限界が容易に見えることは上述の通りです。特に、本書で扱う種子数をモデリングする場合には負の種子数など存在しないので、線形モデルで予測するには無理がありそうだということがわかります。

そこで、線形モデルを一般化した「一般化線形モデル」が登場します。一般化線形モデルをすごく簡単に説明すると、目的変数(y)と説明変数(x)の関係に直線以外を利用できるようにし、ノイズも正規分布以外の分布を利用できるようにしたものと言えます。

書籍では、ポアソン回帰を使って、種子数のモデリングをしています。ポアソン回帰では、データはポアソン分布に従うと仮定します。先の線形モデルでは、正規ノイズであるとしていたところを、ポアソン分布に従って揺らぐと仮定しているということです。また、そのパラメータ $\lambda_i$ は、

${ \displaystyle \lambda_i = \exp{(\beta_0 + \beta_1 x_1)} }$

つまり、

${ \displaystyle \log{\lambda_i} = \beta_0 + \beta_1 x_1 }$

と仮定します。この $\log{\lambda_i}$ は対数リンク関数と呼ばれ、パラメータ $\lambda_i$ を線型結合で表現するための関数です。先の線形モデルでは、パラメータ $\mu$ には関数がかかっておらず、そのような場合「恒等リンク関数」と呼びます。

このモデルを図で示すと以下のようになります。

f:id:hippy-hikky:20190219233550p:plain:w350

以上のように、一般化線形モデルは、データの揺らぎとリンク関数を自由に選択することで、データの生成過程をよりよくモデリングしようという手法になります。線形モデルは、一般化線形モデルの特殊形であるということも容易にわかると思います。

次に係数( $\beta$ )の推定ですが、書籍ではR等のライブラリを利用して計算しています。線形モデルであれば、最小二乗法を使うことで容易に計算できますが、一般化線形モデルでは勾配法などを利用して数値的に解く必要があります。これらの計算については、書籍でも触れていないので割愛します。

「ブラックボックスな統計解析の悪夢」という節で安易なライブラリの使用は避けるように推奨していました。なのに、ここではライブラリ使うのかよと思うかもしれませんが、モデリングの肝心な部分はどのようなモデルを採用するかを考えることです。実際、データを当てはめて係数を推定する部分までスクラッチで書いていたら時間が勿体無いですし、何よりライブラリを使うことで計算速度が速くバグの混入も少ないでしょうから。計算は知っておくに越したことはないので、いずれまとめたいと思います。

【トップに戻る】

第6章 GLMの応用範囲を広げる

本章では、一般化線形モデルの他の応用例として、確率分布、リンク関数の別の例を紹介しています。

ロジスティック回帰

ロジスティック回帰は、私の感覚では線形回帰と並んでよく使われる手法だと思います *3。

ロジスティック回帰では、確率分布に二項分布

${ \displaystyle p(y | N, q) = \binom{N}{y}q^{y}(1-q)^{N-y} }$

、リンク関数に「ロジットリンク関数」を指定します。「ロジットリンク関数」とは何かということですが、まず、二項分布のパラメータは事象の生じる確率 $q$ です。確率なので、0~1の範囲に収まる必要があります。そこで、ロジスティック関数を利用して以下のように表現します。

${ \displaystyle q_i = logistic(z_i) = \frac{1}{1+\exp{(-z_i)}} }$

ここで、 $z_i = \beta_0 + \beta_1 x_1 \cdots$ 。この式を変形すると次のようになります。

${ \displaystyle \log{\frac{p_i}{1-q_i}} = z_i }$

この左辺がロジット関数と呼ばれます。右辺は線型結合なので、これがロジットリンク関数ということです。

このモデルを図で表現すると次のようになります。

f:id:hippy-hikky:20190220004228p:plain:w350

このロジスティック回帰を使った例として、書籍では植物の種子の生存数をモデリングしています。観察する種子数Nに対して生存する種子の数yとし、それが体サイズ等に依存して変わるであろうというモデルです。

現実的な例では、バナー広告のクリック有無をロジスティック回帰でモデリングし、あるユーザがクリックするか否か予測するような問題で応用事例があります。この場合は、二項分布のNが1の場合（ベルヌーイ分布）を確率分布に想定していると思います。

連続量のGLM

書籍では、連続量のGLMとして正規分布の例とガンマ分布の例をあげています。

正規分布のGLMとしては、本記事4.1節で書いたように、確率分布に正規分布をとり、恒等リンク関数を使って表現したものと同等であるということが説明されています。

もう一つの例として、ガンマ分布を使った例を解説します。

例題は、葉の重量を説明変数として、目的変数に花重量としたものです。ガンマ分布(以下の式)は、0以上の連続量をとる確率分布なので、重量はちょうど良い題材ですね。

${ \displaystyle p(y | s, r) = Gamma(s, r) = \frac{r^{s}}{\Gamma{(s)}}y^{s-1}\exp{(-ry)} }$

平均 $\mu$ は $\mu=\frac{s}{r}$ 、分散は $\frac{s}{r^{2}}$ です。この平均 $\mu$ を次のように書き下します。

${ \displaystyle \mu_i = A{x^{b}}_{i} = \exp{(a)}{{x_i}^{b}} = \exp{(a+b \log{x_i})} }$

ここで、 $A=\exp{(a)}$ としています。対数を取ることで、

${ \displaystyle \log{\mu_i} = a+b\log{x_i} }$

となり、対数リンク関数で線形予測子が $a+b\log{x_i}$ です。

割り算値の統計モデリングはやめよう

私自身よくやっていたのですが、データとして得られる説明変数同士から割合のような変数を作り出すというようなことはやられることは多いのかなと思います。「特徴量エンジニアリング」として、いかに良い変数を作り出せるかが良い予測モデルを作るための鍵といった考えもあります。

書籍では、変数を変換しなければいけないような場面はほぼ無いと主張されています。具体的な手段として、人口密度のような変数に興味がある場合には、オフセット項を利用することで解決できるとして例が示されています。

この部分については、よく飲み込めていないので、別途記事にまとめていこうと思います。

【トップに戻る】

まとめ

ここまで、「データ解析のための統計モデリング入門」の一般化線形モデルについてを見てきました。ここまでの内容としては、リンク関数の説明が非常にわかりやすく、線形モデルと一般化線形モデルの差異について理解が深まるのでは無いかなと思います。また、書籍第2章で、データの性質と合う確率分布のまとめが載っています。今後出てくるベイズ推定でも、このまとめはとても有意義ですので是非覚えたいところです。

データ解析の初学者としては、ここまでの内容を理解することが第一歩なのかなと感じます。そういう意味で、本書は初学者にとってまず読むべき本と言えるのかなと思います。

【トップに戻る】

参考文献

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

良い本です。まだ読んでない方は読みましょう。

*1:書籍では「自然の仕組み」と書かれていますが、自然現象だけでなく人工的な現象を扱う場面もあると思います。

*2:対数変換した尤度関数を「対数尤度関数」と呼びます

*3:「回帰」ですが、分類でよく使われる印象です

機械と学習する

統計解析、機械学習について学習したことをまとめていきます

【読後メモ】データ解析のための統計モデリング入門（その1）

はじめに