一般化線形モデル(GLM)は統計解析のフレームワークとしてとにかく便利。
Rでもビルトインの関数から拡張までさまざまなライブラリから提供されている機能だが、
さまざまなライブラリがありすぎてどれを使えばいいのかわかりにくいのと、
さらに一般化線形モデル(GLM)自体にもいろいろな亜種があるため、
どの手法をどのライブラリの関数で実装すればいいかわからなくなる。
そこで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
は行列
- “gaussian”→
alpha = 0
でRidge回帰、alpha = 1
でLasso回帰、Elastic Netはその間で(alpha = 0.5
でhalf and half)。- 目的変数が負の2項分布で正則化回帰を行うには
glmregNB
{mpath}を使う
説明変数間に相関がある場合(PCR、PLSとGLM)
そもそも回帰分析、線形モデルは説明変数が相互に独立であることを前提としている。
しかし説明変数の数が多くなると、説明変数間に相関の見られるものが出てくるようになる。
そのままだと不安定な、おかしな分析結果になる。
一つの解決方法は相関係数に基づいて説明変数を除外する(絞り込む)ことだが、どの変数を使えばいいのか問題にもなる。
その代わりに変数を集約してGLMに適用する方法が
- 主成分回帰(PCR: Principal Component Regression)
- 部分最小二乗回帰(PLS: Partial Least Square)
である。
PCR(主成分回帰)
目的変数とは無関係に、説明変数自体を縮約した結果の変数を回帰分析する手法。
主成分分析を行い、上位の(負荷量の大きい)主成分を説明変数とする。
Rでは
- 説明変数を
prcomp()
{stats}で縮約(主成分分析) - 縮約した変数を説明変数として
glm()
{stats}
PLS(部分的最小二乗回帰)
説明変数間の相関と、目的変数と説明変数の相関を両方考慮して生み出した潜在変数を使った回帰分析。
RではplsRglm()
{plsRglm}
データの加工や分析で使うRの使い方 の記事一覧