一般化線形モデル（GLM）は統計解析のフレームワークとしてとにかく便利。
Rでもビルトインの関数から拡張までさまざまなライブラリから提供されている機能だが、
さまざまなライブラリがありすぎてどれを使えばいいのかわかりにくいのと、
さらに一般化線形モデル（GLM）自体にもいろいろな亜種があるため、
どの手法をどのライブラリの関数で実装すればいいかわからなくなる。

そこでRに実装されている代表的なGLM系の関数と特徴についてまとめてみた。

Contents

1 一般化線形モデルのおさらい
2 RのGLM用関数

一般化線形モデルのおさらい

一般化線形モデルとは

$$
y = g^{-1}(\alpha + \beta_1 x_1 + \beta_2 x_2 + … + \beta_i x_i) + \epsilon
$$

で指定されるモデル。

一般化線形モデルを決定するのは

誤差項 $\epsilon$ の分布→目的変数の分布から決まる
線形予測子 $\alpha + \beta_1 x_1 + \beta_2 x_2 + … + \beta_i x_i$ →最初に考えるモデル式（説明変数の組み合わせ）
リンク関数 $g()$（線形予測子と推定値の関係式）→誤差項の分布から大体決まる

の3つ。

手順は

モデル式（の候補）を洗い出す→パラメータを推定する→モデルを評価する

モデルの評価は

anova(fit)でモデル自体が有意かどうかを見る
AICなどで当てはまりを見る
par(mfrow=c(2,2));plot(fit)で残差の分布を見る
- Residuals vs Fitted: 残差全体の挙動
- Normal Q-Q: 残差が正規分布でないGLMでは見る必要ない
- Scale-Location: 残差の大きさ
- Residuals vs Leverage: 外れ値の確認
検証用データで予測精度を見る

RのGLM用関数

glm(): デフォルトで組み込まれている一般化線形モデルの関数
glmnet(): パッケージ{glmnet}に含まれるElastic Net（Ridge+Lasso）を実装した一般化線形モデル

`glm` {stats}

デフォルトで組み込まれている関数

glm(モデル式, family = 目的変数の分布, data = データフレーム)
# リンク関数も指定する場合
glm(モデル式, family = 目的変数の分布(link = リンク関数), data = データフレーム)

familyは誤差項ではなく目的変数の確率分布を指す。
対応しているfamilyとリンク関数

binomial(link = "logit")
- 目的変数が2値変数
- 対応しているリンク関数は
  - ‘logit’（ロジスティック回帰／ロジットモデル）
  - ‘probit’（プロビットモデル）
  - ‘cauchit’ ( = Cauchy)
  - ‘log’
  - ‘cloglog’ (complementary log-log)
- リンク関数の違い: http://www.karlin.mff.cuni.cz/~kulich/vyuka/pokreg/R/glm_binary_links.html
poisson(link = "log")
- 目的変数が0以上の離散変数、分散がそこまで大きくない
- 対応しているリンク関数は
  - ‘log’（対数線形モデル）
  - ‘identity’
  - ‘sqrt’
gaussian(link = "identity")
- 目的変数が正規分布に従う
- 対応しているリンク関数は
  - ‘identity’（線形予測子＝推定値）
  - ‘log’（対数正規）
  - ‘inverse’
Gamma(link = "inverse")
- 目的変数が0以上の連続変量
- 対応しているリンク関数は
  - ‘inverse’
  - ‘identity’
  - ‘log’
inverse.gaussian(link = "1/mu^2")
- 対応しているリンク関数は
  - ‘1/mu^2’
  - ‘inverse’
  - ‘identity’
  - ‘log’
quasi(link = "identity", variance = "constant")
quasibinomial(link = "logit")
quasipoisson(link = "log")

glm() {stats}で対応していないタイプの主要なモデル

目的変数が負の2項分布（ポアソンよりゼロが多い）→glm.nb() {MASS}
多項ロジット→mlogit() {mlogit}
Coxの比例ハザードモデル→coxph() {survival}

ポアソン回帰と負の2項回帰

カウントデータの回帰ではポアソン分布が代表的に用いられるが、
期待値と分散が一致していることが前提の分布なので、分散が大きくなる場合はポアソン回帰が当てはまりにくくなる。
特にゼロが多いデータでは過分散の状態になる。

→そこでdispersion test（過分散の検定）を行う。

帰無仮説：過分散でない＝ポアソン回帰が適当

が棄却されたら負の2項分布を使うのが適切になる。

過分散の検定を行う関数

dispersiontest(ポアソンglmオブジェクト) {AER}
odTest(glm.nbオブジェクト) {pscl}

これらは引数にとるオブジェクトの種類が異なる。

参考
http://tjo.hatenablog.com/entry/2013/09/23/232814

カウントデータの回帰

カウントデータでは上記のように

分散≒平均→ポアソン回帰glm(..., family = poisson(link = "log"), ...)
分散が大きい→負の2項回帰glm.nb() {MASS}

がよく用いられるが、ゼロが多いデータでは上記のほかにGLMではないが

ゼロ切断モデル（Hurdle model、一つのカウントデータの分布で、ハードルを超えなければ0として扱う）→hurdle() {pscl}
ゼロ過剰モデル（Zero-inflated model、ゼロになるかどうかの2項回帰＋ゼロにならなければ通常のカウントデータの回帰）→zeroinfl() {pscl}

がある。ゼロの発生する原理に基づいて使い分けるといいかも。

参考
http://www.slideshare.net/logics-of-blue/2-6-25621115

`glmnet` {glmnet}

回帰分析において過学習を防ぐため、モデルに正則化項を付けたGLMをあてはめる。
代表的な正則化にはLasso（L1正則化）とRidge（L2正則化）があり、それらの混合をElastic Netという。
glmnet {glmnet}はこのElastic Net扱う関数（ライブラリ）。

glmnet(説明変数, 目的変数, family = 目的変数の分布, alpha = 1)
# データフレームはこうやって指定できる
glmnet(as.matrix(data.train[,-cv]), data.train[,cv], family = 目的変数の分布, alpha = 1)

glmnetの特徴

Elastic Net（Ridge+Lasso）を実装
多項ロジットを扱える（glm() {stats}にはなかったためmlogit() {mlogit}などが必要だった）
データフレームではなくMatrix形式（sparse matrixなど）を使う
familyの指定の仕方がglm() {stats}と異なる。
- “gaussian”→yは連続量のベクトル
- “binomial”→yは2値変量のベクトル
- “poisson”→yは離散量（≧0）のベクトル
- “multinomial”: 多項ロジット→yは因子ベクトル
- “cox”: Coxの比例ハザードモデル→yは(時間,ステータス)のn×2行列（Surv() {survival}で生成できる）
- “mgaussian”: 多変量正規分布→yは行列
alpha = 0でRidge回帰、alpha = 1でLasso回帰、Elastic Netはその間で（alpha = 0.5でhalf and half）。
目的変数が負の2項分布で正則化回帰を行うにはglmregNB {mpath}を使う