70. 単回帰分析


単回帰分析・直線回帰

関数 lm() と predict() を用いて単回帰分析を行うことができる.例えば以下のようなアメリカ桜のデータ trees があったとする. このデータに関して散布図を描いた後,単回帰分析を行う.

 
 plot(Volume ~ Girth, data = trees)                           # 散布図を描く
 result <- lm(Volume ~ Girth, data = trees)                   # 回帰分析を行う
 abline(result)                                               # 推定回帰直線を描く
 summary(result)                                              # 分析結果の要約

Call:
lm(formula = Volume ~ Girth, data = trees)
Residuals:
    Min      1Q  Median      3Q     Max 
-8.0654 -3.1067  0.1520  3.4948  9.5868 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -36.9435     3.3651  -10.98 7.62e-12 ***
Girth         5.0659     0.2474   20.48  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
Residual standard error: 4.252 on 29 degrees of freedom
Multiple R-Squared: 0.9353,     Adjusted R-squared: 0.9331 
F-statistic: 419.4 on 1 and 29 DF,  p-value: < 2.2e-16

結果は -36.9435 が切片項,5.0659 が傾きの推定値となっている.ここで,予測区間と信頼区間をプロットしてみる.

 
 new <- data.frame(Girth = seq(8, 24, 0.5))
 result.pre <- predict(result, new, interval="prediction")    # 予測区間
 result.con <- predict(result, new, interval="confidence")    # 信頼区間
 matplot(new$Girth, cbind(result.pre, result.con), lty=c(1,2,2,3,3), type="l")

プロット結果は以下のようになる.

 

 

関数 lsfit() による最小二乗法

関数 lsfit() により,y = X b + e の b を推定する.y には従属変数(被説明変数)のベクトルを与え,X には独立変数(計画行列,説明変数;それぞれがベクトル) を cbind したものを与える.つまり,この行列の各列が説明変量で各列が観測値となる.lsfit() の書式は以下の通り.

 
 lsfit(x, y, wt=NULL, intercept=TRUE, tolerance=1e-07, yname=NULL)  

引数の説明は以下の通り.

引数

機能

x

説明変数の行列.

y

被説明変数のベクトル.

wt

重みをつける場合はここに指定する.

intercept

デフォルトは TRUE で,この場合は自動的に計画行列(説明変数の行列)に定数項が加わる.定数項を加えたくない場合は FALSE を指定する.

例として,データ trees に関して最小二乗法を行う.

 
 x <- trees$Height;  y <- trees$Volume
 z <- lsfit(x, y)      # 推定結果を z に格納
 print(z)              # coefficients, residuals, intercept, qr を出力

出力の説明は以下の通り.

引数

機能

coef

回帰係数ベクトル.intercept=T  (デフォルト)  ならば第 1 要素が回帰式の定数項とる.

residuals

残差ベクトル.

intercept

引数interceptと同じ値で,説明行列が定数項に対応する列を含んでいるか否かを示す.

qr

x に定数項に対応する列を追加した行列を QR分解した結果.この要素の副要素 qt は t(Q) %*% y を表し ( t(Q) は Q の転置行列) ,この要素の qt 以外の副要素は関数 qr が返す要素をそのまま保持している.

回帰による予測値を求める場合は,出力 z について以下の様に入力すればよい.

 
 y - z$residuals

 [1] 20.91087 13.19412 10.10742 23.99757 37.88772 40.97442 14.73747 28.62762
 [9] 36.34437 28.62762 34.80102 30.17097 30.17097 19.36752 28.62762 27.08427
[17] 44.06112 45.60447 22.45422 11.65077 33.25767 36.34437 27.08427 23.99757
[25] 31.71432 37.88772 39.43107 36.34437 36.34437 36.34437 47.14782

関数 ls.diag(z) で残差の標準偏差やハット行列 H の対角成分などが得られる.

また,重相関係数 R =  (予測値の分散) / (被説明変数の分散)  を求める場合には以下の様に入力すればよい.

 
 sqrt(var(y - z$residuals)/var(y)) 

[1] 0.5982497

以下のように lsfit() の結果を abline() に渡すことによって,回帰直線を描くことも出来る.

 
 plot(x, y)
 abline(z, col=2)