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

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

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

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

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

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

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

　　　　

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

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

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

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

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

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

図3　散布図と回帰直線

　図3の中の直線の関数

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

　各引数の意味を簡潔に次に示す。
　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 を用いて返還する。

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

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

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

図4　回帰分析の主な結果

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

　データ(真のモデル)を

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

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

図5　回帰直線と残差

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

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

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

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

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

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

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

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

　式の中の ST 、SR 、Se は次のように定義されている。

偏差の平方和：

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

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

(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　回帰診断図

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