１．統計モデルと回帰分析
人間のみならず多くの動物は、学習機能を持っており、過去の経験から得られた知識・規則を今後の行動に生かしている。データ解析では、既知のデータから規則を導き出し、その規則を用いて未知の部分を説明したり、予測・推測したりする。例えば、体重(kg)と身長(m)のデータを用いた肥満度の指標BMI(Body Mass Index)を考えよう。
この式は既知のデータと医学の知識から導いた関係式である。通常BMIの値が22前後であれば標準であり、25前後であると肥満傾向で、35を超えると極度肥満であるといわれている。このようなデータ間の関係の規則を関数式で表することができると、その数式を用いて未知なことを予測・推測することが可能である。このようなデータから導いた規則をここでは統計モデルと呼ぶ。
図１に平成14年度に総務省が発行した「情報通信白書」の中の、わが国のインターネットに関する報告を示す。横軸が時間軸で、縦軸がその年度のインターネットの利用者数である。図の右側には平成17年の予測値が示されている。このような時間の前後の順に並べた時系列データでは、しばしば現在までのデータの統計モデルを用いて今後を予測・推測する。
図1 平成14年度の「情報通信白書」での予測の例
(平成14年度の「情報通信白書」より)
統計モデルの説明のため、身近な例として、身長と体重との関係を何らかの式で表すことを考えることにする。その具体的な関係がわからないので、抽象的な関数記号
で次のように表すことにする。
この式の中の「身長」を説明変数、「体重」を被説明変数、あるいは目的変数と呼ぶ。もし説明変数を
、被説明変数を
で表すと上記の式は次のような式になる。
このような説明変数を用いて被説明変数を説明する関係を求める統計モデルを通常回帰分析と呼ぶ。説明変数が１つの場合は単回帰分析、説明変数が複数の場合は重回帰分析と呼ぶ。
図2に単回帰分析の例として2種類の説明変数と被説明変数の関係を示す。横軸が説明変数、縦軸が被説明変数である。
図2(a)では、説明変数と被説明変数との関係を直線で大まかな傾向を示している。このような直線関係でモデル化する回帰分析を線形回帰分析と呼ぶ。
図2(b)の説明変数と被説明変数のような非直線的関係でモデル化する回帰分析を非線形回帰分析と呼ぶ。
(a)線形関係
(b)非線形関係
図2 説明変数と被説明変数との関係
２．単回帰分析
単回帰分析は説明変数が1つである回帰分析であり、略して単回帰と呼ぶ。単回帰でも線形と非線形に分けられるが、特別な説明がない限り線形回帰を指す。
ここで言う線形とは、説明変数
と被説明変数
の関係を式
で表すことが可能であることである。式の中の
を切片、
は直線の傾きに関する値で、説明変数
の係数と呼ぶ。
例えば、表1のような体重、身長データがあるとする。
表1 体重、身長のデータ
|
身長 |
体重 |
|
165 |
50 |
|
170 |
60 |
|
172 |
65 |
|
175 |
65 |
|
170 |
70 |
|
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.647で、説明変数
の係数
は1.6809である。回帰直線は、適当に描いたものではなく、既知のデータに最もフィット(fit)が良いと考えられる直線である。回帰分析は最適と思われる回帰式の係数を求めるのが一つの主な目的である。Ｒで線形回帰分析を行う関数の中で最も広く使用されているのはlmである。
lm(formula,
data, weights, subset, na.action)
各引数の意味を簡潔に次に示す。
formulaは回帰式に用いる被説明変数と説明変数、定数項を用いるか用いないかなどのモデルの形式である。例えば、回帰式を
にしたいときにはformulaの部分はy~xと指定し、定数項を用いない回帰式
にしたいときには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計算された結果の表示及び回帰分析に関する関連操作行うコマンドを表２に示す。
表２ 回帰分析の関連操作を行うコマンド
|
コマンド |
内 容 |
使 用 例 |
|
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)
図４ 回帰分析の主な結果
まず関数summaryが返した結果について順を追って説明を行い、次に予測値の求め方やグラフによる視覚的な考察などについて説明する。
(1) 残差
データ(真のモデル)を
回帰直線(推測モデル) を
としたとき、説明変数
が
を取った場合の両値の差
を残差(residuals)と呼ぶ。図4にデータ、回帰直線、残差を示す。
図5 回帰直線と残差
summary(takei.lm)のResidualsでは残差の最小値、第1四分位数、中央値、第2四分位数、最大値が返される。次のコマンドで全ての残差を返すことができる。
(2) 回帰式と回帰係数
回帰直線
の係数(Coefficients)である
、
は、残差の2乗の和
上記のsummary(takei.lm)で返されている結果のCoefficientsの部分が回帰係数及びそれに関連する情報である。Interceptの行が切片
の推測値、その標準誤差、
値、
値である。その次の行が説明変数「身長」の係数と関連の情報である。
回帰式は返されている回帰係数を用いて次のように構築する。
回帰係数の標準誤差、
値は次の式で定義されている。
残差の平方和： 
残差の不偏分散：
式の中
はデータの標本数、
は説明変数の数である。
係数
の標準誤差：
係数
の標準誤差：
係数
の
値： 
係数
の
値： 
値に対応する
値は関数ptで求めることが可能である。
この
値は「回帰係数がゼロである」という仮説検定の統計量である。
値が通常よく使われている有意水準0.05(5％)、0.1(10％)、0.05(0.5％)より小さいときには、出力結果の
値の右にそれぞれ1つの星”*”,2つの星”**”、3つの星”***”で印をつける。
回帰係数は次のコマンドで返される。
>
coefficients(taikei.lm)#あるいは
>
taikei.lm$coefficients
(3) 決定係数
回帰直線がどの程度データにフィットしているかを評価する指標として決定係数(coefficient of determination) がある。決定係数を通常
で示す。関数lmでは決定係数(Multiple
R-Squared)と調整済み決定係数(Adjusted
R-squared)が返される。決定係数と調整済み決定係数が1に近づくほど回帰直線がデータによくフィットしていると判断する。
決定係数と調整済み決定係数は次の式で定義されている。式の中の
は標本の数で、ここの
は説明変数の数である。単回帰では説明変数が１つであるので
である。
決定係数： 
調整済み決定係数：
式の中の
、
、
は次のように定義されている。
偏差の平方和： 
偏差の平方和： 
残差の平方和： 
平方和の関係： 
関数lmが返す結果の中の決定係数、調整済み決定係数の下部に
値とその
値が返される。
値と
値は「全ての回帰係数がゼロである」という帰無仮説の検定等計量である。
値は決定係数から求めることができ、この
値は自由度
、
の
分布に従う。
(5) グラフによる分析
単回帰の場合は、説明変数と非説明変数の散布図に求めた回帰直線を加えることでデータの傾向を概観することができる。次にそのコマンドと結果を示す。
>plot(taikei.F)
>abline(taikei.lm)
図６ データの散布図と回帰直線
また回帰分析ではしばしば残差を視覚的に分析を行う。Rでは回帰分析の残差などを視覚的に考察するための環境が用意されている。回帰分析の結果を次のコマンドのように関数plotに代入するだけで図７に示す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の残差が相対的に大きいことがわかる。
② 残差の正規Q-Qプロット
正規Q-Qプロットはデータの正規性を考察するためにデータを視覚化する方法である。Rでは関数qqplotを用いてQ-Qプロットを作成することができる。データが正規分布に従うと、点が直線上に並べられる。通常の回帰分析では、残差が標準正規分布に従うという仮定の下で行っている。
図7(b)は標準化した残差のQ-Qプロットである。ここで用いた例では標本データの数が少ないのでその正規性に関する議論には大きな意味がない。
③ 残差の平方根プロット
図7(c)は標準化した残差の絶対値の平和根を縦軸にし、予測値を横軸とした散布図である。この図の目的も残差の変動状況を考察することである。図から個体6、8、9の変動が相対的に大きいことが読みとられる。
④ Cookの距離のプロット
Cookの距離は一種の距離の測度であり、Rには関数cooks.distanceが用意されている。Cookの距離は回帰分析における影響度が大きいデータの検出などに多く用いられている。Cookの距離は全てデータ用いた場合と１つのデータを除いた後求めた回帰式による予測値を用いた場合との食い違いに関する距離の測度である。Cookの距離が大きいとそのデータが回帰式による予測値に大きく影響していることを意味する。よって、Cookの距離が大きいデータは異常値である可能性がある。Cookの距離が0.5以上であれば大きいと言われている。図7(d)では個体8のCookの距離が比較的に大きいが0.5以下であるので個体8が異常値であるとはいえない。
診断図を作成するplotではwhich=nの引数を用いて、4種類の図の中の1つのみを作成することができる。ここのnは1、2、3、4でそれぞれ上記の①、②、③、④に対応する。例えばplot(taikei.lm,which=1)であれば残差対予測値の散布図が作成される。
回帰診断図をさらに加工したいときには、次のように個別のデータを求めて散布図を作成したほうが便利である。
①
残差対予測値の散布図
plot(resid(taikei.lm)) #plot(taikei.lm$resid)
abline(h=0)
②
Q-Qプロット
qqnorm(resid(taikei.lm))
#qqnorm(taikei.lm$resid)
qqline(resid(taikei.lm)) #qqline(taikei.lm$resid)
③
Cookの距離
plot(cooks.distance(taikei.lm))