R言語

【R言語】一般化線形モデル 関数glmの使いかた

  1. HOME >
  2. R言語 >

【R言語】一般化線形モデル 関数glmの使いかた

スポンサーリンク

一般化線形モデルを行う方法を解説していきます。

R言語で回帰分析を行う方法については、以下の記事で紹介しています。

【R言語】関数lmによる線形回帰 単回帰分析・重回帰分析 回帰分析で予測する

R言語の回帰分析の実行方法を紹介していきます。回帰分析の実行方法だけでなく、回帰分析により作ったモデルでテストデータの予測も行っていきます。 実行方法では、回帰分析の結果の見方やその抽出方法を詳しく解 ...

続きを見る

一般化線形モデルを用いることで、母分布が正規分布でなくても線形モデルを構築することが可能となります。

例えば、目的変数が0と1のような二値である場合でも、回帰(予測)することが可能となります。

また、この記事で紹介しているプログラミングコードは以下のリンクからダウンロードすることができます。

RStudioで解析したい方へ
¥3,278 (2022/06/02 22:52時点 | Amazon調べ)
\楽天ポイント5倍セール!/
楽天市場
\ポイント5%還元!/
Yahooショッピング
Rstudioで統計解析がしたい方へ
¥3,058 (2022/06/01 17:32時点 | Amazon調べ)
\楽天ポイント5倍セール!/
楽天市場
\ポイント5%還元!/
Yahooショッピング

関数glm

R言語で一般化線形モデルを実行するには、関数glmを用います。

基本的な使いかたは回帰分析で用いたlmと変わりません。

lmとの違いとして、引数に誤差構造リンク関数を指定する点があります。

関数glmは以下の引数から成ります。

引数は以下の通りです。

関数glmの引数
formula回帰に用いるモデル(式)offsetオフセット項の指定
family誤差構造とモデルに用いられるリンク関数。control適合のプロセスを制御するパラメータのリスト
data回帰に用いるデータセットmodel返り値にモデルを含むかどうか
weightsprior weights. 事前に重みを与えたい場合数値ベクトルで与える。methodモデルを適合する際に用いるアルゴリズム
subset回帰する際に用いる観測値の部分集合x返り値に計画行列Xを含むかどうか
na.actiondataがNAを含む場合にとる行動。y返り値に目的変数yを含むかどうか
start回帰係数の初期値singular.ok計画行列Xが非正則の場合、回帰を行うか
etastart線形予測しの初期値contrastsオプションのリスト。詳細はcontrasts.arg。
mustart線形予測値の平均値の初期値  

上で紹介した引数のうち、familyについて詳しく見ていきます。

引数familyを変更することで、誤差構造とそれに対応するリンク関数を変更することが可能となります。

誤差構造リンク関数を以下にまとめました。

familyオブジェクト
familylink
binomial"logit", "probit", "cauchit", "log", "cloglog"
gaussian"identity", "log", "inverse"
Gamma"inverse", "identity", "log"
inverse.gaussian"1/mu^2", "inverse", "identity", "log"
poisson"log", "identity", "sqrt"
quasi"logit", "probit", "clolog", "identity", "inverse", "log","1/mu^2", "sqrt"
quasibinomial"logit", "probit", "cauchit", "log", "cloglog"
quasipoisson"log", "identity", "sqrt"

例として、目的変数がc(1, 3, 6, 3, 9, 2, 3, 4)のような離散値である場合、母集団にポアソン分布</spanを指定します。

また、c(1, 0, 1, 1, 0, 1, 1)のような二値である場合、二項分布を指定します。

この後、より具体的な実行方法をみていきます。

glmの実行例

誤差構造がポアソン分布

まず、glmの使用例として、ポアソン回帰の実行方法についてみていきます。

誤差構造がポアソン分布に従う場合の回帰分析の実行方法を解説していきます。

簡単に言うと、目的変数が「1, 2, 3, 4」といった離散値であるときの回帰分析の手順について説明します。

ここでは、以下の線形モデルから成るポアソン回帰の実行方法について触れます。

リンク関数にlogを使っていることに注意してください。

ポアソン回帰

ポアソン回帰は次の線形モデルで表されます。

\begin{align}\log \lambda_i = \log \mathrm{E}[Y_i | \boldsymbol{x}]= \beta_0 + \beta_1x_{1i} + \cdots + \beta_kx_{ki}, \ \ i = 1, \ldots, n,\end{align}

ここに\(\lambda_i\)は事象の発生回数の期待値であり、\(\boldsymbol{x}_i = (x_{1i}, \ldots, x_{ki})^T\)である。

では早速、ポアソン回帰の実行方法をみていきましょう。

データセットの用意

目的変数がポアソン分布に従う場合のglmの実行方法について解説します。

次のようなデータセットポアソン分布に従う列countsを持つデータを用います。

次のように、dataは離散変数のcountsと「1, 2, 3」を水準に持つ因子型の列outcomeとtreatmentから成ります。

また、パッケージggplot2を用いてcountsとoutcomeの経験分布及び箱ひげ図を描画すると 、outcomeの水準によってcountsが変化することが分かります。

上を実行することで、次のような経験分布と箱ひげ図がプロットされます。outcomeが1のcountsは2, 3のものより比較的大きいことが読み取れます。

経験分布と箱ひげ図1

また、treatmentとcountsの経験分布と箱ひげ図も以下のように引数を変更することでプロットすることができます。

経験分布と箱ひげ図2

treatmentが1のcountsは3のものと比べて、ばらつきが非常に小さいことが分かります。

回帰分析の実行

回帰分析を実行する前に、まずdataをトレインデータとテストデータに分けます。後にテストデータを用いた予測をするために、dataを学習用と検証用に分割します。

今回は簡単に前から9個のデータをトレインデータ、残りのデータをテストデータに分割します。

改めて、次の9行から成るデータに対し回帰分析を実行していきます。

ポアソン回帰を実行するには、次のように関数glmの引数を与えます。

1つ目の引数formulaと2つ目の引数dataに関しては、関数lmと全く同じです。formulaにはモデルを表す式を代入し、dataには回帰に用いるデータセットを代入します。

これらの引数に加えて、3つ目の引数familyにpoissonを代入することでポアソン回帰ができます。

regressionを参照するとモデルの回帰係数やAICがあるのが分かります。

次にこれらの結果をまとめる方法について解説していきます。

次のように、関数summayによって、これらの解析結果をリストの要素としてまとめることができます。

ポイント

summaryRegを参照すると、Coefficientsの欄の(intercept)とoutcome2の一番右に「*」がついているのが分かります。p値が0.05未満であり、統計的に有意である変数であることが言えます。

また、解析結果を得る方法についてみていきます。summaryReg$. . . でsummaryReg中の要素を参照することができます。

ポイント

上の例では、summaryReg$coefficientsで回帰係数やp値などが格納されたデータフレームに取得し、summaryReg$aicでモデルのaicを取得しています。

最後に、これらの解析結果をcsvファイルに出力しましょう。

関数write.csvを用いて、上で取り上げた解析結果をcsvファイルに出力します。

データフレームcoeffsとaicを含むデータフレームresultsを以下のように与えます。

resultsは以下のようになっています。

以下のように、write,csvでデータフレームresultsをcsvファイルへと出力します。

上を実行すると、作業ディレクトリ中に"ポアソン回帰.csv"が保存されているのが確認できます。

csvファイルを開くと次の画像のように適切に出力されているのが分かります。

ポアソン回帰

テストデータを用いた予測

次に、先ほど作ったテストデータをモデルにあてはめて、テストデータを予測していきます。

テストデータvalidationDataの値をモデルにあてはめるには、関数predictを用います。

次のようにpredictの引数typeを"response"にすることで、予測値を得ることができます。

また、予測値とテストデータの比較方法についてみていきます。

これらの予測値を実際の値と比較するために、次のようなデータフレームpredictedDataを構成します。

counts以外のoutcomeとtreatmentはvalidationDataのものと全く同じです。

次に、このpredictedDataとvalidationDataを縦方向に結合させます。

まず、どこまでのデータが予測値であるか実際の値であるかを区別するために、labelsというベクトルを作ります。

最後に、rbinddata.frameを用いて次のようなテストデータと予測値が結合したデータフレームを作成します。

ggplotを用いて、labelsの値についてcountsの値を描画します。