R言語

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

  1. HOME >
  2. R言語 >

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

スポンサーリンク

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

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

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

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

続きを見る

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

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

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

関数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の値を描画します。

上を実行すると次のようなプロットが出力されます。

テストデータと予測値の比較

青がテストデータの値で、赤が予測値です。

だいたい予測できていることが見て取れます。

誤差構造がベルヌーイ分布(ロジスティック回帰)

次に、誤差構造がベルヌーイ分布に従う場合について考えてみます。

ここでは、二値分類で有名なロジスティック回帰の実行方法を解説していきます。

分類問題で多項ロジットなどもありますが、ここでは目的変数が二値であるものについて説明します。

ロジスティック回帰

ロジスティック回帰は次の線形モデルで表されます。

\begin{align}\log p_i = \log\left(\cfrac{p_i}{1-p_i}\right) = \beta_0 + \beta_1x_{1i} + \cdots + \beta_kx_{ki},\ \ i = 1, \ldots, n,\end{align}

ここに

\begin{align}p_i = \mathrm{Pr}\{Y_i=1\}.\end{align}

また、\(p_i ,\ i = 1, \ldots, n\)に関して次のようにも表される。

\begin{align}p_i = \mathrm{Pr}\{Y_i = 1\} = \cfrac{1}{1 + e^{\beta_0 + \beta_1x_{1i} + \cdots + \beta_kx_{ki}}}.\end{align}

0と1の二値から成る目的変数について、\(Y_i=1\)となる確率を\(p_i\)とし、この\(p_i\)に対してロジット変換を行うことで\(-\infty\)から\(\infty\)をとる変数に変換しているのが見て取れます。

難しい話はここまでにして、Rでロジスティック回帰を実行する例をみていきましょう。

データセットの用意

まず、ロジスティック回帰を行うために、パッケージISLRをインストールします。

ISRLには、次のようなDefaultというデータセットが含まれています。

このデータには"Yes"と"No"から成るdefaultという列があるのが分かります

このdefaultを、ロジスティック回帰を用いて予測していきます。

次に準備として、目的変数defaultの"Yes"と"No"をそれぞれ1と0に変換します。

注意ポイント

1と0に変換した後はfactor型ではなくnumeric型にしましょう。

ここで、balancedefaultをプロットしてみます。次のようにdefaultの値に関してbalanceをプロットしましょう。

defaultが0のときbalanceは小さく、defaultが1のときbalanceは大きい傾向があります。

defaultとbalance

defaultはbalanceによって変化すると考えられるので、defaultをbalanceで回帰することにします。

回帰分析の実行

早速、defaultをbalanceで回帰していきます。

準備として、dataをトレインデータとテストデータに分割します。後にテストデータのあてはめも行いたいので、次のようにtrainDatavalidationDataという新しい変数を与えます。

ロジスティック回帰を実行するには関数glmの引数familyを次のように与えます。

誤差構造が二項分布であり、リンク関数はロジットであるため、family = binomial(logit)と与えればよいです。

ロジスティック回帰を実行し、regressionを参照すると次のような解析結果が得られます。

関数summaryを用いることで、ポアソン回帰と同様に回帰係数やAICを得ることができます。

Coefficientsの欄を見ると、統計検定の結果が載っているのも分かります。(Intercept)とbalanceの一番右に「*」がついているのが分かります。p値が0.05未満であり、これらは統計的に有意な変数であることが言えます。

summaryRegは回帰係数やAICをリストの要素として持ちます。次のようにして、回帰係数やAICを得ることができます。

ポアソン回帰で解説した通り、次のようにして解析結果csvファイルへ出力することができます。

モデルを用いた予測

最後に、先ほど作成したvalidationDataをモデルにあてはめていきます。

テストデータの予測値を得るには、ポアソン回帰と同様に関数predictを用います。

予測値と実際の値をプロットするために、次のようテストデータと予測値をもつvalidPredDataを作ります。

validPredDataが作れたら、パッケージtidyverseggplot2を用いて、ロジスティック回帰曲線とテストデータをプロットします。

上を実行すると、次の画像が描画されます。

ロジスティック回帰曲線

青線がモデルへあてはめたときの値(予測値)で、黒点がテストデータの値です。

balanceが小さいとdefaultは0に近く、大きいと1に近いためある程度予測できていることが分かります。

さらに、判別の精度として誤判別の割合を計算してみましょう。予測値が0.5より大きいか以下であるかでTRUEとFALSEに分割したデータをclassifierとします。

次に、関数tableを用いて、defaultの値とclassifierの値に関して2×2の集計表を作成します。

集計表が作れたら、次のようにして誤判別の割合を計算することができます。

0.0248の割合で誤判別をしていることがわかります。このことから、balanceという変数でdefaultをほとんど予測することが可能であることが言えます。

誤差構造がベルヌーイ分布(多重ロジスティック回帰)

最後に、多重ロジスティック回帰について解説します。

さきほどは、defaultという二値変数をbalanceで回帰しましたが、今回はdefaultを他の全ての変数で回帰していきます。

回帰分析の実行

次のように、関数glmformulaを指定することで簡単に多重ロジスティック回帰を実行できます。

関数lmの使いかたで紹介しましたが、formulaを「default ~.」と指定するとdefault以外の変数でdefaultを回帰することができます。

また補足ですが、balanceincomeでdefaultを回帰したい場合は、次のようにします。

関数summaryを使って、次のように解析結果をリストに変換することができます。Coefficientsの欄を見ると、(Intercept)とstudentYesとbalanceの一番右に「*」がついているのが分かります。これら3つの変数はp値が0.05未満であり、有意な変数であることがいえます。

解析結果は以下のようにして取り出すことができます。

また、これらの解析結果をcsvファイルへ出力したい場合は、次のようにしてできます。

作業ディレクトリ中に次の"多重ロジスティック回帰.csv"が保存されているのが分かります。

多重ロジスティック回帰

モデルを用いた予測

最後に、モデルの予測性能についてみていきます。

ロジスティック回帰では、balanceでdefaultを単回帰していましたが、今回はdefault以外の変数を全てモデルに組み込みました。

予測性能に変化があったかどうか見ていきます。

単回帰のモデルを用いた予測と同様に、関数predictを用いてvalidationDataの値をモデルにあてはめます。

予測値predictedYの値が0.5より大きいか以下であるかでTRUEとFALSEに分類します。これをclassifierとします。

次に、関数tableを用いて、defaultの値とclassifireの値に関して集計表を作ります。

上の集計表から、誤判別の割合を算出してみます。誤判別の割合は次のようにして簡単に求めることができます。

0.0244であることから、非常に判別の精度が良いことがいえます。また、単回帰の時と比較すると誤判別の割合が若干減少したことがわかります。

ポイント

このように、モデルに新たな変数を組み込むことで予測精度を向上させることができます。また、AICについても多重ロジスティック回帰の結果の方が小さいことが分かります。

まとめ

R言語で一般化線形モデルを用いた回帰分析を行う方法を解説しました。

目的変数が正規分布に従わない場合に、回帰分析を実行する方法を学びました。

また、テストデータのあてはめ等のモデルへのデータのあてはめの方法も見ていきました。

ぜひ、今回の記事を解析に役立ててください。

スポンサーリンク

  • この記事を書いた人
  • 最新記事

usagi-san

統計学とゲームとかをメインに解説していくよ。 数式とかプログラミングコードにミスがあったり質問があったりする場合はコメントで受け付けます。すぐに対応します。

-R言語
-, ,