[連載] フリーソフトによるデータ解析・マイニング 第14回

Rと回帰分析

1.統計モデルと回帰分析

 人間のみならず多くの動物は、学習機能を持っており、過去の経験から得られた知識・規則を今後の行動に生かしている。データ解析では、既知のデータから規則を導き出し、その規則を用いて未知の部分を説明したり、予測・推測したりする。例えば、体重(kg)と身長(m)のデータを用いた肥満度の指標 BMI (Body Mass Index) を考えよう。

 この式は既知のデータと医学の知識から導いた関係式である。通常 BMI の値が22前後であれば標準であり、25前後であると肥満傾向で、35を超えると極度肥満であるといわれている。このようなデータ間の関係の規則を関数式で表することができると、その数式を用いて未知なことを予測・推測することが可能である。このようなデータから導いた規則をここでは統計モデルと呼ぶ。
 図1に平成14年度に総務省が発行した「情報通信白書」の中の、わが国のインターネットに関する報告を示す。横軸が時間軸で、縦軸がその年度のインターネットの利用者数である。図の右側には平成17年の予測値が示されている。このような時間の前後の順に並べた時系列データでは、しばしば現在までのデータの統計モデルを用いて今後を予測・推測する。

図1 平成14年度の「情報通信白書」での予測の例
(平成14年度の「情報通信白書」より)

 統計モデルの説明のため、身近な例として、身長と体重との関係を何らかの式で表すことを考えることにする。その具体的な関係がわからないので、抽象的な関数記号 f で次のように表すことにする。

体重 = f (身長)

 この式の中の「身長」を説明変数、「体重」を被説明変数、あるいは目的変数と呼ぶ。もし説明変数をx 、被説明変数をy で表すと上記の式は次のような式になる。

 このような説明変数を用いて被説明変数を説明する関係を求める統計モデルを通常回帰分析と呼ぶ。説明変数が1つの場合は単回帰分析、説明変数が複数の場合は重回帰分析と呼ぶ。
  図2に単回帰分析の例として2種類の説明変数と被説明変数の関係を示す。横軸が説明変数、縦軸が被説明変数である。
  図2(a)では、説明変数と被説明変数との関係を直線で大まかな傾向を示している。このような直線関係でモデル化する回帰分析を線形回帰分析と呼ぶ。
  図2(b)の説明変数と被説明変数のような非直線的関係でモデル化する回帰分析を非線形回帰分析と呼ぶ。

    

  (a) 線形関係                    (b) 非線形関係

図2 説明変数と被説明変数との関係

2.単回帰分析

 単回帰分析は説明変数が1つである回帰分析であり、略して単回帰と呼ぶ。単回帰でも線形と非線形に分けられるが、特別な説明がない限り線形回帰を指す。
 ここで言う線形とは、説明変数x と被説明変数y の関係を式

y = a + bx

で表すことが可能であることである。式の中のa を切片、b は直線の傾きに関する値で、説明変数x の係数と呼ぶ。
 例えば、表1のような体重、身長データがあるとする。

表1 体重、身長のデータ
身長 体重
165
50
170
60
172
65
175
65
170
70
172
75
172
75
183
80
187
85
180
90
185
95

 まずデータセットを次のように作成する。

> taikei = matrix(0,10,2)
> taikei[,1]<-c(165,170,172,175,170,172,183,187,180,185)
> taikei[,2]<-c(50,60,65,65,70,75,80,85,90, 95)
> colnames(taikei)<- c("身長","体重")
> taikei

身長 体重
[1,] 165 50
<後略>

  コマンド plot(taikei) で図3のような直線がない散布図を作成することができる。図3に身長を横軸、体重を縦軸にした散布図を示す。図3の直線は、散布図の点の傾向(トレンド)を表すもので、回帰直線と呼ぶ。

図3 散布図と回帰直線

 図3の中の直線の関数

体重 = - 222.1647 + 1.6809 × 身長

を線形回帰式、あるいは回帰直線と呼び、式の中の切片a は−222.647で、説明変数x の係数b は1.6809である。回帰直線は、適当に描いたものではなく、既知のデータに最もフィット (fit) が良いと考えられる直線である。回帰分析は最適と思われる回帰式の係数を求めるのが一つの主な目的である。Rで線形回帰分析を行う関数の中で最も広く使用されているのは lm である。

lm(formula, data, weights, subset, na.action)

 各引数の意味を簡潔に次に示す。
 formula は回帰式に用いる被説明変数と説明変数、定数項を用いるか用いないかなどのモデルの形式である。例えば、回帰式を y = a 0+a 1x にしたいときには formula の部分はy~x と指定し、定数項を用いない回帰式 y = a 1x にしたいときには formula の部分は y~ 1+x (あるいは y~ x -1 ) と指定する。
 data は回帰分析に用いたデータセットの名前である。例えば、data = taikei、あるいは data = を省略して taikei のみでもよい。
 wights は用いる説明変数に重みをつける引数である。初心者は無視してもよい。特に設定しない場合は重みを付けない。
  subset はデータセットの中の一部分を用いる際に、用いる部分を明示するための添字ベクトルの引数である。指定しない場合は、全てのデータを用いる。
 na.action は欠損値扱いを指定する引数である。指定がない場合は、欠損値のデータを除いたデータを用いて計算を行う。よって、na.action を指定しない場合は、予測値や残差の数はデータセットから欠損値を除いた数と等しい。
  関数 lm に用いるデータ形式はデータフレームである。マトリックス形式からデータフレーム形式には次のように関数 data.frame を用いて返還する。

> taikei.F<- data.frame(taikei)

 takei.F を用いた線形回帰関数 lm の使用例を示す。

> taikei.lm<- lm(体重~身長,data=taikei.F)

 関数 lm 計算された結果の表示及び回帰分析に関する関連操作行うコマンドを表2に示す。

表2 回帰分析の関連操作を行うコマンド
 コマンド  内容 使用例
 print  要約より簡単な結果  print(taikei.lm)
 summary  回帰分析結果の要約  summary(taikei.lm)
 coef  回帰係数  cofe(taikei.lm)、taikei.lm$coef
 fitted  用いたデータの予測値  fitted(taikei.lm)、taikei.lm$fitted 
 deviance  残差の平方和  deviance(taikei.lm)
 anova  回帰係数の分散分析  anova(taikei.lm)
 predict  新たのデータに対する予測値   predict(taikei.lm)
 plot  回帰診断プロット  plot(taikei.lm)

 次に関数 summary による結果について説明する。

> summary(taikei.lm)

図4 回帰分析の主な結果

 まず関数 summary が返した結果について順を追って説明を行い、次に予測値の求め方やグラフによる視覚的な考察などについて説明する。

(1) 残差

 データ(真のモデル)を

y = a + bx + e

回帰直線(推測モデル)を

としたとき、説明変数xxi を取った場合の両値の差 を残差 (residuals) と呼ぶ。図4にデータ、回帰直線、残差を示す。


図5 回帰直線と残差

 summary(takei.lm) の Residuals では残差の最小値、第1四分位数、中央値、第2四分位数、最大値が返される。次のコマンドで全ての残差を返すことができる。

> residuals(taikei.lm)
#あるいは
> taikei.lm$residuals

(2) 回帰式と回帰係数

 回帰直線 の係数 (Coefficients) であるab は、残差の2乗の和

を最小化する連立方程式を解くことで次の解が得られる。

 上記の summary(takei.lm) で返されている結果の Coefficients の部分が回帰係数及びそれに関連する情報である。Intercept の行が切片a の推測値、その標準誤差、t 値、p 値である。その次の行が説明変数「身長」の係数と関連の情報である。
 回帰式は返されている回帰係数を用いて次のように構築する。

体重 = - 222.1647 + 1.6809 × 身長

 回帰係数の標準誤差、t 値は次の式で定義されている。

  残差の平方和:    

  残差の不偏分散:  

 式の中n はデータの標本数k は説明変数の数である。

  係数a の標準誤差:  

  係数b の標準誤差:  

  係数at 値:       

  係数bt 値:       

 t 値に対応するp 値は関数 pt で求めることが可能である。
 このt 値は「回帰係数がゼロである」という仮説検定の統計量である。p 値が通常よく使われている有意水準0.05(5%)、0.1(10%)、0.05(0.5%)より小さいときには、出力結果のp 値の右にそれぞれ1つの星“*”,2つの星“**”、3つの星“***”で印をつける。
 回帰係数は次のコマンドで返される。

> coefficients(taikei.lm)  #あるいは
> taikei.lm$coefficients

(3) 決定係数

 回帰直線がどの程度データにフィットしているかを評価する指標として決定係数 (coefficient of determination) がある。決定係数を通常R 2 で示す。関数 lm では決定係数 (Multiple R-Squared) と調整済み決定係数 (Adjusted R-squared) が返される。決定係数と調整済み決定係数が1に近づくほど回帰直線がデータによくフィットしていると判断する。
 決定係数と調整済み決定係数は次の式で定義されている。式の中のn は標本の数で、ここのk は説明変数の数である。単回帰では説明変数が1つであるのでk =1 である。

  決定係数:         

  調整済み決定係数:  

 式の中の STSRSe は次のように定義されている。

  y 偏差の平方和:    

  偏差の平方和:    

  残差の平方和:      

  平方和の関係:      

 関数 lm が返す結果の中の決定係数、調整済み決定係数の下部にF 値とそのp 値が返される。F 値とp 値は「全ての回帰係数がゼロである」という帰無仮説の検定等計量である。F 値は決定係数から求めることができ、このF 値は自由度kn -k -1 のF 分布に従う。

(4) 予測値

 回帰式で計算される値を予測値、あるいは推測値と言う。関数 summary では予測値は返さない。予測値を返すためにはRでは関数 predict を用いる。考察の便利のため回帰分析に用いたデータ、予測値、残差を次のように一覧表にして示す。

> 予測値<-predict(taikei.lm)
> 残差<-residuals(taikei.lm)
> data.frame(taikei.F,予測値,残差)

身長 体重 予測値 残差
1
165 50 55.17854 -5.178535
2
170 60 63.58288 -3.582877
<中略>
9
180 90 80.39156 9.608440
10
185 95 88.79590 6.204098

(5) グラフによる分析

 単回帰の場合は、説明変数と非説明変数の散布図に求めた回帰直線を加えることでデータの傾向を概観することができる。次にそのコマンドと結果を示す。

> plot(taikei.F)
> abline(taikei.lm)

図6 データの散布図と回帰直線

 また回帰分析ではしばしば残差を視覚的に分析を行う。Rでは回帰分析の残差などを視覚的に考察するための環境が用意されている。回帰分析の結果を次のコマンドのように関数 plot に代入するだけで図7に示す4つの異なるグラフが返される。
 関数 par は作図の環境設定関数である。用いた引数 mfrow=c(2,2) は1画面に2行2列の図を作成する指定である。回帰分析の結果を関数 plot で作成した図を回帰診断図と呼ぶ。

> par(mfrow=c(2,2)) #par(mfrow=c(2,2),oma=c(2,2,2,1),mar=c(5,4,3,2))
> plot(taikei.lm) # plot(taikei.lm,cex=1.5,pch=21,bg="blue",col="blue")

図7 回帰診断図

@ 残差とフィット値のプロット

  図7(a)は残差とフィット値の散布図である。ここで言うフィット値は予測値である。図の横軸が予測値、縦軸が残差となっている。図から残差の全体像を概観することができる。この図では個体 6、8、9 の残差が相対的に大きいことがわかる。

A  残差の正規Q-Qプロット

  正規Q-Qプロットはデータの正規性を考察するためにデータを視覚化する方法である。Rでは関数 qqplot を用いてQ-Qプロットを作成することができる。データが正規分布に従うと、点が直線上に並べられる。通常の回帰分析では、残差が標準正規分布に従うという仮定の下で行っている。
  図7(b)は標準化した残差のQ-Qプロットである。ここで用いた例では標本データの数が少ないのでその正規性に関する議論には大きな意味がない。

B 残差の平方根プロット

 図7(c)は標準化した残差の絶対値の平和根を縦軸にし、予測値を横軸とした散布図である。この図の目的も残差の変動状況を考察することである。図から個体 6、8、9 の変動が相対的に大きいことが読みとられる。

C Cookの距離のプロット

 横軸は梃子値で、縦軸は標準化残差である。赤い点線でクックの距離0.5を示している。Cook の距離は一種の距離の測度であり、Rには関数 cooks.distance が用意されている。Cook の距離は回帰分析における影響度が大きいデータの検出などに多く用いられている。Cook の距離は全てデータ用いた場合と1つのデータを除いた後求めた回帰式による予測値を用いた場合との食い違いに関する距離の測度である。Cook の距離が大きいとそのデータが回帰式による予測値に大きく影響していることを意味する。よって、Cook の距離が大きいデータは異常値である可能性がある。Cook の距離が0.5以上であれば大きいと言われている。図7(d)では個体8の Cook の距離が比較的に大きいが0.5以下であるので個体8が異常値であるとはいえない。
  診断図を作成する plot では which = n の引数を用いて、4種類の図の中の1つのみを作成することができる。ここのnは1、2、3、4でそれぞれ上記の @、A、B、C に対応する。例えば plot (taikei.lm,which=1) であれば残差対予測値の散布図が作成される。
 回帰診断図をさらに加工したいときには、次のように個別のデータを求めて散布図を作成したほうが便利である。

@ 残差対予測値の散布図

> plot(resid(taikei.lm)) #plot(taikei.lm$resid)
> abline(h=0)

A Q-Qプロット

> qqnorm(resid(taikei.lm)) #qqnorm(taikei.lm$resid)
> qqline(resid(taikei.lm)) #qqline(taikei.lm$resid)

B Cookの距離

> plot(cooks.distance(taikei.lm))