一般化線形モデル（GLM）は統計解析のフレームワークとしてとにかく便利。Rでもビルトインの関数から拡張までさまざまなライブラリから提供されている機能だが、さまざまなライブラリがありすぎてどれを使えばいいのかわかりにくいのと、さらに一般化線形モデル（GLM）自体にもいろいろな亜種があるため、どの手法をどのライブラリの関数で実装すればいいかわからなくなる。
そこでRに実装されている代表的なGLM系の関数と特徴についてまとめてみた。
目次
一般化線形モデルのおさらい
一般化線形モデルとは
$$
y = g^{-1}(\alpha + \beta_1 x_1 + \beta_2 x_2 + … + \beta_i x_i) + \epsilon
$$
で指定されるモデル。
一般化線形モデルを決定するのは
- 誤差構造：目的変数の分布
- 線形予測子 $\alpha + \beta_1 x_1 + \beta_2 x_2 + … + \beta_i x_i$ →説明変数の組み合わせ
- リンク関数 $g()$（線形予測子と推定値の関係式）→目的変数の分布から大体決まる
の3つ。
手順は
- 説明変数（の候補）を洗い出す
- モデルを構築する（パラメータを推定する）
- モデルを評価する
モデルの評価は、modelを構築したモデルのオブジェクトとすると
anova(model)でモデル自体が有意かどうかを見る- AICなどで当てはまりを見る
- 検証用データで予測精度を見る
- 線型モデル（
lm()）のみpar(mfrow=c(2,2));plot(model)で残差の分布を見る- Residuals vs Fitted: 残差全体の挙動
- Normal Q-Q: 残差が正規分布に従うかどうか
- Scale-Location: 残差の大きさ
- Residuals vs Leverage: 外れ値の確認
RのGLM用関数
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(..., family = poisson(link = "log"), ...)が代表的に用いられるが、期待値と分散が一致していることが前提の分布なので、分散が大きくなる場合はポアソン回帰が当てはまりにくくなる。特にゼロが多いデータでは過分散の状態になる。
そこでdispersion test（過分散の検定）を行う。「帰無仮説：過分散でない＝ポアソン回帰が適当」が棄却されたら違うモデルを使うのが適切になる。
過分散の検定を行う関数
dispersiontest(ポアソンglmオブジェクト){AER}odTest(glm.nbオブジェクト){pscl}
これらは引数にとるオブジェクトの種類が異なる。
参考
http://tjo.hatenablog.com/entry/2013/09/23/232814
ゼロが多い（過分散の）場合、
- 負の2項分布を想定したGLM
glm.nb(){MASS} - ゼロ切断モデル（Hurdle model、一つのカウントデータの分布で、ハードルを超えなければ0として扱う）→
hurdle(){pscl} - ゼロ過剰モデル（Zero-inflated model、ゼロになるかどうかの2項回帰＋ゼロにならなければ通常のカウントデータの回帰）→
zeroinfl(){pscl}
が適切になる。ゼロの発生する原理に基づいて使い分けるといいかも。
参考
http://www.slideshare.net/logics-of-blue/2-6-25621115
説明変数が多い回帰の問題対策
そもそも回帰分析、線形モデルは説明変数が相互に独立であることを前提としている。しかし説明変数の数が多くなると、説明変数間に相関の見られるものが出てくるようになる。そのままだと構築したモデルの回帰係数が不安定な解釈できないものになる（多重共線性）。また説明変数が多くなると、学習データに対する当てはまりはいいが、実のデータに対して予測をすると不正確になることがある（過学習）。
正則化回帰
glmnet() {glmnet}
正則化回帰とはモデルの説明変数に過度に依存しないように、回帰係数が大きくなった場合にペナルティがかかるように調整した回帰モデル。代表的な正則化にはLasso（L1正則化）とRidge（L2正則化）があり、それらの混合をElastic Netという。
- Lasso回帰では多重共線性のある説明変数の係数をゼロにする。変数を減らす目的で使われる。
- Ridge回帰では多重共線性のある変数の係数を残す（係数を小さくする）。過学習対策。
glmnet(説明変数, 目的変数, family = 目的変数の分布, alpha = 1)
# データフレームはこうやって指定できる
glmnet(as.matrix(data.train[,-label]), data.train[,label], 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で半々）。- 目的変数が負の2項分布で正則化回帰を行うには
glmregNB{mpath}を使う
手動で説明変数を絞り込む
回帰係数ごとにVIFを計算する。vif() {car}で算出できる。値が大きいと係数が不安定になっている。10以上でNGとすることが多い。
- VIF=3で偏相関係数0.8
- VIF=5で偏相関係数0.9
- VIF=10で偏相関係数0.95
特にVIFの大きな説明変数を除外する手がある。
主成分（相関のない合成変数）を使った回帰
相関関係のある説明変数群から無相関の合成変数を作り、その合成変数を説明変数とした回帰を行う。合成変数の生成方法によって以下の２種類の方法がある。
- 主成分回帰（PCR: Principal Component Regression）
- 部分最小二乗回帰（PLS: Partial Least Square）
PCR（主成分回帰）
目的変数とは無関係に、説明変数自体を主成分分析によって無相関化、縮約した合成変数を作る。説明変数に対して主成分分析を行い、負荷量上位の主成分を説明変数とする。
Rでは
- 説明変数を
prcomp(){stats}で縮約（主成分分析） - 縮約した変数を説明変数として
glm(){stats}
問題点として、主成分の分散が大きいからといって目的変数との関連が大きいとは限らない。下位の負荷量の主成分でも目的変数との関連が大きいことはあるし、目的変数の変動の様子と合成変数の作り方には関係がない。
PLS（部分的最小二乗回帰）
plsRglm() {plsRglm}
目的変数の変動を吸収するために最適化された合成変数の作り方を採用する。具体的には目的変数と説明変数の共分散を最大化する主成分を生成する（PCRでは主成分の分散を最大化する主成分を生成）。
特殊な回帰
多項ロジット
mlogit() {mlogit}
複数の選択肢から1個だけ選択するケース（通常のロジットモデルは2個から1個を選択する）
Coxの比例ハザードモデル
coxph() {survival}
イベントの発生を予測する
1人につきイベントが1回しか発生しない
coxph(Surv(time, event) ~ x1 + x2 + x3, data=mydata)
1人につきイベントが複数回発生する（recurrent event model）
coxph(Surv(start, stop, event) ~ x1 + x2 + x3, id=customer_id, data=mydata)
折れ線回帰
説明変数の区間によって目的変数との関係性（回帰係数）が異なる場合がある。たとえば訪問回数が特定の回数まではポジティブな影響があるが、それを超えると横ばいになる、逆にネガティブになるなどのケースである。この場合、説明変数を区間に区切って別の変数として扱う方法がある。たとえば「訪問回数」という説明変数を
- 訪問回数10回までの訪問回数
- 訪問回数11回以上の訪問回数
という2個の変数に分けるということである。これを折れ線回帰（piecewise regressionまたはsegmented regression）という。
説明変数の閾値を手動で指定して折れ線回帰を実行する
x1を10を閾値に2分割、x2を30と60を閾値に３分割して折れ線回帰を実行する例
# 分割用の関数を作成
piecewise <- function(x, lower, upper=Inf){
# lowerは区間の下限で、upperは区間の上限。
if (missing(lower)) {
return(pmin(x, upper))
} else {
return(pmin(x-lower,upper-lower)*(x>lower))
}
}
# 分割した変数を生成
mydata$x1_rng1 <- piecewise(mydata$x1, upper=10)
mydata$x1_rng2 <- piecewise(mydata$x1, lower=10)
mydata$x2_rng1 <- piecewise(mydata$x2, upper=30)
mydata$x2_rng2 <- piecewise(mydata$x2, lower=30, upper=60)
mydata$x2_rng3 <- piecewise(mydata$x2, lower=60)
# 分割した変数を説明変数として回帰を実行
lm_model <- lm(y ~ x1_rng1 + x1_rng2 + x2_rng1 + x2_rng2 + x2_rng3, data=mydata)
説明変数の閾値を自動計算して折れ線回帰を実行する
segmented() {segmented}
- あらかじめ当該変数を使った
lm(),glm(),Arima()のいずれかのモデルを生成 - そのモデルと区間分割する変数と分割点の初期値を手動で与えて
segmented()関数を実行
するとそこから最適な閾値を探索して折れ線回帰を実行する。
lm_model <- lm(y ~ x1 + x2, data=mydata)
o <- segmented(lm_model, seg.Z = ~ x1 + x2, psi = list(x1=c(10), x2=c(30,60)))
slope(o)
データの加工や分析で使うRの使い方 の記事一覧
