3.8 Rによるモデリング
ここでは引き続きHeightsデータセットを例に単回帰分析のモデリングをRで行う方法を紹介しいく．Heightsデータセットには母親の身長と18歳以上の娘の身長が1375ペア分記録されており，以下では母親の身長から娘の身長を説明するという回帰モデルを考えていく．
3.8.1 データの準備
まず，次のコードを実行してデータセットが読み込まれることを確認して欲しい．
## mheight dheight
## 1 59.7 55.1
## 2 58.2 56.5
## 3 60.6 56.0
## 4 60.7 56.8
## 5 61.8 56.0
## 6 55.5 57.9もしエラーが出る場合はalr4パッケージがインストールできていない可能性があるので，
以下のコードを実行した上で上記のコードを再度実行してみて欲しい．
3.8.2 分布の確認
既に何度か見ているものもあるが，以下について確認しておこう．
3.8.2.1 各変数のデータ型の確認
Heightsはデータフレームになっている．
## [1] "data.frame"このオブジェクトについて，各変数の要約値などを調べるためにsummary関数がを利用して集計しておく．
## mheight dheight
## Min. :55.40 Min. :55.10
## 1st Qu.:60.80 1st Qu.:62.00
## Median :62.40 Median :63.60
## Mean :62.45 Mean :63.75
## 3rd Qu.:63.90 3rd Qu.:65.60
## Max. :70.80 Max. :73.10Heightsに含まれる二つの変数はどちらも数値として扱われていることを確認できた．
次に分布を確認していこう．ただし変数の数種類でサンプルサイズも1000程度であるため
散布図行列というものが便利である．
Rのデフォルトの関数としてpairsが用意されているので利用しよう．．
このように各変数についての散布図をマトリックスで表示してくれる．また，変数自身のヒストグラムも見るにはGGallyパッケージのggpairs関数がさらに便利である．以下のコードを実行してまずはGGallyパッケージを読み込んで欲しい．
パッケージを読み込んだらGGally関数を使ってみよう．．
このように対角成分にヒストグラムを，上三角成分には相関係数などを表示してくれている． また離散型の変数の場合は散布図ではなくグループ別の箱ひげ図にしてくれたりと何かと便利である．
こちらのページが詳しいので興味のある人は参照されたい．
上記の散布図行列をみると，相関係数は0.49程度だが，それぞれの分布も正規分布に近く異常な外れ値も見受けられないのでこのまま学習に進んでも良さそうなことがわかる．
3.8.3 学習
目的変数\(y_i\)をdheight，説明変数\(x_i\)をmheightとしたモデルを考えよう．
Heightsデータセットは`1375行あるので
\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i, \hspace{3mm} i=1,2,\ldots,1375 \]
となる．次に回帰係数パラメーター\(\beta_0, \beta_1\)を推定するが，Rにおいてはlm関数を使うことで簡単に実行できる．
上記が実行できたらfitというオブジェクト名で回帰モデルの学習結果が保存されている．
fitの中身を確認してみよう．
##
## Call:
## lm(formula = dheight ~ mheight, data = Heights)
##
## Coefficients:
## (Intercept) mheight
## 29.9174 0.5417Callという部分で実行されたモデルが表示され，Coefficients:という部分で推定された回帰係数パラメーターの値が表示されている．
ここでは(Intercept)が\(\hat \beta_0\)，mheightが\(\hat \beta_1\)に対応している．つまり\(\hat \beta_0 = 29.9174, \hat \beta_1 = 0.5417\)と推定されている．
これをモデルに置き換えると
\[ y_i = 29.9174 + 0.5417 x_i + \varepsilon_i, \hspace{3mm} i=1,2,\ldots,1375 \]
となる．さらに詳細の情報を確認するにはfitオブジェクトにsummary関数を適用する．
##
## Call:
## lm(formula = dheight ~ mheight, data = Heights)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.397 -1.529 0.036 1.492 9.053
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.91744 1.62247 18.44 <2e-16 ***
## mheight 0.54175 0.02596 20.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.266 on 1373 degrees of freedom
## Multiple R-squared: 0.2408, Adjusted R-squared: 0.2402
## F-statistic: 435.5 on 1 and 1373 DF, p-value: < 2.2e-16残差（Residuals）の情報や回帰係数パラメーターの有意水準・決定係数などの情報も得ることができる．
このモデルにおいては決定係数\(R^2 = 0.2408\)と推定されている．
Adjusted R-squaredというのは調整済み決定係数と呼ばれるもので，こちらについては今回は割愛する．
決定係数の値をみる限り，そこまで当てはまりが良くなさそうなことがわかる．
最後に回帰直線の表示の方法を紹介しよう．x,yの散布図を出力した後にabline(fit)とすることで，回帰直線を散布図に上書きできる．
以下の例ではcol="red"というパラメータを指定することで回帰直線を赤色で出力するように命令している．
Figure 3.8: データと回帰直線
3.8.4 残差の評価
それでは，学習したモデルの残差についてみていこう．既に紹介したように残差プロットとQ-Qプロットの様子も見てモデルの当てはまりを確認する．
残差プロット
残差プロットについては，残差の回帰直線はほぼ横軸と一致しているように見えるためばらつきには問題がなさそうである． しかしいくつか残差の絶対値が大きい値を見つけることができる．
Q-Qプロット
Q-Qプロットを見ると，分布の裾の方で説明しきれていない傾向があることがわかる．
3.8.5 考察
以上の結果を踏まえると，残差の傾向としては概ねモデルの仮定に沿っているようだが，決定係数\(R^2\)は\(0.24\)程度とモデルの改善の余地があると言える．実際，母親の身長とその娘の身長のデータであることを考えると，母親の身長が高ければ娘の身長も高いという傾向はありそうだが，データ上は例外も多く見つかることから，母親の身長だけでは説明しきれない変数があるのかもしれない．例えば食生活や睡眠時間なども人間の成長には大きな影響を与え，これらが身長にも影響しているような仮説を立てても良いだろう．
モデルに改善の余地が見つかれば，改良するためにどんなデータを変数として追加すべきなのか，そもそもモデルの仮定はあっているかなどを考察しモデリングと学習を繰り返していくことになる．