Rで線形単回帰分析 - matsuou1の日記

次回のTokyo.Rの開催が近づいてきたので、前回の復習を兼ねてRで回帰分析をやってみます。

今回は最も単純な線形単回帰分析を行います。

回帰分析の流れ

回帰式を求める意義があるか検討する（説明変数と目的変数のグラフを作成する等）
回帰式を求める
回帰式の精度を確認する
回帰係数の検定を行う
信頼区間と予測区間を求める

回帰式を求める意義があるか検討

無相関なデータに対しても、数学的には回帰式が求められるため、検討しておくことは重要です。

データはマンガでわかる統計学回帰分析編のデータを使用してみます。
ある喫茶店のアイスティーの売り上げとその日の最高気温についてのデータです。

> norns
     temperture icetea
8/22         29     77
8/23         28     62
8/24         34     93
8/25         31     84
8/26         25     59
8/27         29     64
8/28         32     80
8/29         31     75
8/30         24     58
8/31         33     91
9/01         25     51
9/02         31     73
9/03         26     65
9/04         30     84

単相関係数の確認

1 から 1 の間の実数値をとり、1 に近いときは正の相関があるといい、-1 に近ければ負の相関があるという。

0 に近いときは相関は弱い。
相関係数は順序尺度であり間隔尺度ではないので、例えば「相関係数が0.2と0.4であることから、後者は前者より2倍の相関がある」などと
言うことはできない。

Rではcor関数を使用して、単相関係数を算出することができる。

> cor(norns$temperture, norns$icetea)
[1] 0.906923

今回は1に近いので、強い正の相関関係があると考えられる。

散布図の確認

散布図をplotし、回帰直線を引く意味がありそうかどうかを確認します。

> plot(norns)

大体、右肩上がりとなっているため、回帰直線を引く意味がありそうな散布図と考えられる。

回帰分析

ここまでデータについて確認してきましたが、回帰分析をやる意味がありそうなので、早速回帰分析をやってみましょう。
回帰分析は、lm関数を使うことで簡単にできます。

まずは、lm関数の書式を確認します。

lm関数

書式

lm(formula, data, subset, weights, na.action,method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,singular.ok = TRUE, contrasts = NULL, offset, ...

引数

引数名	必須	説明
formula	○	当てはめられるモデルのシンボリックな記述式
data	○	モデル中の変数を含むオプションのデータフレーム
subset		当てはめ過程で使われる観測値の部分集合を指定
weights		当てはめ過程で使われるオプションの重みベクトル
na.action		データが NA を含むときそれを処理する関数
method		使われるメソッド
model		x, y, qr 論理値．もし TRUE なら当てはめの成分 (モデルフレーム，モデル行列，目的変数，分解) が返される
singular.ok		論理値. もし FALSE なら特異な当てはめはエラーとされる
contrasts		オプションのリスト
offset		これは当てはめの途中で，線形予測子に含められるべきことが事前に知られている成分を指定するのに使える

formulaの書式について

y~x
y~x-1

1は、通常の線形単回帰分析で使用するが、2は定数項がない場合に使用する。
自動車のガソリンの消費量と走行距離など、説明変数が0の場合、目的変数が必ず0になる強い仮定がある場合にのみ使用する。

戻り値

戻り値名	説明
coefficients	係数の名前付きベクトル
residuals	残差，つまり目的変数から当てはめ値を引いたもの
fitted.values	当てはめられた平均値
rank	当てはめ線形モデルの数値計算によるランク
weights	(重みつき当てはめだけに該当) 指定された重み
df.residual	残差自由度
call	呼出し式
terms	使われた terms オブジェクト
contrasts	(関係するときだけ) 使われた対比
xlevels	(関係するときだけ) 当てはめで使われた因子の水準の記録
y	もし要求されたなら使われた目的変
x	もし要求されたなら使われたモデル行列
model	もし要求されたなら (既定)，使われたモデルフレーム

実行結果

では、実際にlm関数を使って、回帰分析をやってみます。

> norns.lm <- lm(icetea~temperture , data=norns)
> summary(norns.lm)
Call:
lm(formula = icetea ~ temperture, data = norns)

Residuals:
   Min     1Q Median     3Q    Max 
-8.037 -5.693  2.094  4.409  8.225 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -36.3612    14.6873  -2.476   0.0292 *  
temperture    3.7379     0.5012   7.457 7.66e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 5.709 on 12 degrees of freedom
Multiple R-squared: 0.8225,	Adjusted R-squared: 0.8077 
F-statistic: 55.61 on 1 and 12 DF,  p-value: 7.661e-06

コード自体は2行しかないので、実行は超簡単です。
分析結果のサマリは、summary関数で見ることができます。

以下、summaryの各項目について見て行きましょう。

Residuals

残差の四分位数
（四分位数とは、分布の左側からの相対頻度の合計値（累積相対頻度）が 25%、50%、75% になる値のこと）

   Min     1Q Median     3Q    Max 
-8.037 -5.693  2.094  4.409  8.225

Coefficients

            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -36.3612    14.6873  -2.476   0.0292 *  
temperture    3.7379     0.5012   7.457 7.66e-06 ***

Estimate：係数
Std.Error：標準誤差
t value：t値
pr(>|t|)：p値

回帰分析においては y=a+bx において　b=0　となることはあってならず、b=0となると　y=a　となり、x は y に影響を与えないことになってしまいます。
それでは回帰分析が成立しないので（ここでいう b=0 のことを帰無仮説とよぶ）、b=0となることを棄却できるかどうか考えなくてはなりません。

そこで　ｐ値を見ます。ｐ値は帰無仮説が発生する確率なので、0.01以下ならばその推定値は99%の有意水準、0.1以下ならば90%の有意水準ということになります。

また、有意であるならば横に　*（星）がつきます。何パーセントかは
Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
をみればわかります。

つまり、星は'***'で、p値が小さければ小さいほど良いということです。

Multiple R-squared , Adjusted R-squared

Multiple R-squared: 0.8225,	Adjusted R-squared: 0.8077

Multiple R-squared：寄与率、決定係数
Adjusted R-squared：調整済み寄与率、調整済み決定係数
寄与率には説明変数が多くなる程、値が大きくなる欠点があるので、自由度調整済み寄与率を使用する。
（線形単回帰分析では、そんなに気にする必要はないが、重回帰分析では意識して扱う必要がある。）

これは１に近いほど予測値と観測値が近くなることを表し、モデルの正確性を表す。

値に対する基準はないが、「0.5以上」を一つの目安として扱う。

回帰式

上記結果から以下の回帰式が得られる。

y = 3.7379x -36.3612

plotしてみると、こんな感じになる。

> plot(norns)
> abline(norns.lm , lwd=1 , col="blue")

予測値

回帰式で計算される値を予測値（推測値）は、summary関数では返却されない。
予測値を返すためには、predict関数を用い算出する。

> norns.predict <- predict(norns.lm) # 予測値

> norns.residuals <- residuals(norns.lm) # 残差

> data.frame(norns, norns.predict , norns.residuals)
     temperture icetea norns.predict norns.residuals
8/22         29     77      72.03744        4.962555
8/23         28     62      68.29956       -6.299559
8/24         34     93      90.72687        2.273128
8/25         31     84      79.51322        4.486784
8/26         25     59      57.08590        1.914097
8/27         29     64      72.03744       -8.037445
8/28         32     80      83.25110       -3.251101
8/29         31     75      79.51322       -4.513216
8/30         24     58      53.34802        4.651982
8/31         33     91      86.98899        4.011013
9/01         25     51      57.08590       -6.085903
9/02         31     73      79.51322       -6.513216
9/03         26     65      60.82379        4.176211
9/04         30     84      75.77533        8.224670

得られた回帰式の妥当性の検証

残差を視覚的に分析するために、回帰診断図を表示してみます。
回帰診断図は、回帰分析結果をplotすることで表示できます。

> par(mfrow=c(2,2))
> plot(norns.lm)

以下に、それぞれのグラフが何を示しているかを確認しましょう。

残差とフィット値のプロット（左上）

横軸が予測値、縦軸が残差。
残差の全体像を概観するために使用する。

残差の正規Q-Qプロット（右上）

データの正規性考察するために使用する。
データが正規分布に従っている場合は、直線上に並ぶ。

残差の平方根プロット（左下）

残差の変動状況を考察するために使用する。
標準化した残差の絶対値の平方根を縦軸にし、予測値を横軸にした散布図。

残差と影響力プロット（梃子値とクック距離）（右下）

１つのデータがモデルの当てはまりへの影響力を測るために使用する。
クックの距離が0.5を超えると影響力あり、1を超えると特異に大きい。
横軸は梃子値で、縦軸は標準化した残差。点線でクックの距離0.5を示している。

信頼区間と予測区間を求める

信頼区間とは、「測定・解析を全く同じやり方で100回繰り返したら、95回は回帰直線はこの範囲を通ると考えられる」区間のことです。

信頼区間は、predict関数にinterval="confidence"を指定することで求められる。

> norns.con <- predict(norns.lm , interval="confidence")

> plot(x.plot , norns.con[,1] , type="l",xlim=c(23,35),ylim=c(40,110),ylab="" , xlab="")
> par(new=T)
> plot(x.plot , norns.con[,2] , type="l",xlim=c(23,35),ylim=c(40,110),ylab="" , xlab="" , col="red")
> par(new=T)
> plot(x.plot , norns.con[,3] , type="l",xlim=c(23,35),ylim=c(40,110),ylab="" , xlab="" , col="blue")
> par(new=T)
> plot(norns$temperture , norns$icetea , type="p",xlim=c(23,35),ylim=c(40,110),ylab="icetea" ,xlab="temperture")

グラフで見るとこんな感じです。

この結果からいえることは、例えば、同じやり方で気温が32度の日のアイスティー売り上げの平均値を見てみると、95%の確率で78.69 〜 87.81の間に収まりますということになります。

予測区間とは、「新たな測定を１００回繰り返したら、95回はこの範囲に収まると考えられる」区間のことです。

予測区間は、predict関数にinterval="prediction"を指定することで求められる。

> new <- data.frame(temperture= seq(23 , 35 , 1 ))
> norns.pre <- predict(norns.lm , new , interval="prediction")
> plot(x.plot , norns.pre[,1] , type="l",xlim=c(23,35),ylim=c(40,110),ylab="" , xlab="")
> par(new=T)
> plot(x.plot , norns.pre[,2] , type="l",xlim=c(23,35),ylim=c(40,110),ylab="" , xlab="" , col="red")
> par(new=T)
> plot(x.plot , norns.pre[,3] , type="l",xlim=c(23,35),ylim=c(40,110),ylab="" , xlab="" , col="blue")
> par(new=T)
> plot(norns$temperture , norns$icetea , type="p",xlim=c(23,35),ylim=c(40,110),ylab="icetea" ,xlab="temperture")

この結果から言えることは、例えば、気温が32度の日のアイスティー売り上げは、95%の確率で70.00 〜 96.49の間に収まりますということになります。

単位根検定

ランダムウォーク（単位根過程）同士を回帰分析すると、見せかけの回帰が発生するため、単位根検定を行い、変数が単位根でないかどうか確認しておく。
使用するのはPhillips-Perron検定。

> PP.test(norns2[,1])

	Phillips-Perron Unit Root Test

data:  norns2[, 1] 
Dickey-Fuller = -8.0442, Truncation lag parameter = 2, p-value = 0.01

> PP.test(norns2[,2])

	Phillips-Perron Unit Root Test

data:  norns2[, 2] 
Dickey-Fuller = -6.0224, Truncation lag parameter = 2, p-value = 0.01

それぞれのp-valueが0に近いため単位根過程ではなく、問題なさそう。

回帰式から定数項をはずした場合

次に、回帰式から定数項をはずした場合についても、やってみます。
formulaで指定する式を代えるだけなので簡単です。

> norns2.lm <- lm(icetea~temperture - 1 , data=norns)
> summary(norns2.lm)

Call:
lm(formula = icetea ~ temperture - 1, data = norns)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.591  -4.358  -1.103   5.888   8.890 

Coefficients:
           Estimate Std. Error t value Pr(>|t|)    
temperture  2.50366    0.06149   40.72 4.26e-15 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 6.742 on 13 degrees of freedom
Multiple R-squared: 0.9922,	Adjusted R-squared: 0.9916 
F-statistic:  1658 on 1 and 13 DF,  p-value: 4.263e-15

サマリを確認してみても、定数項は表示されいないです。
得られる回帰式は以下の通り。

y = 2.50366x

注意点

定数項をはずした場合（formula=y~x-1）の寄与率が、定数項を含んだ場合（formula=y~x）の寄与率より
高くなる場合もあるが、ある場合とない場合で計算方法が違っているので比較はできない。
実際に今回の場合は、前者の寄与率が0.9916、後者の寄与率0.8077がとなっているが、Residualsを見ても
散布図＋回帰式を見ても後者の方が精度が高い。

    Min      1Q  Median      3Q     Max 
-8.037   -5.693   2.094   4.409   8.225 (formula=y~x-1)
-11.591  -4.358  -1.103   5.888   8.890 (formula=y~x)

青：定数項なし、赤：定数項あり