関数 lm() と predict() を用いて単回帰分析を行うことができる．例えば以下のようなアメリカ桜のデータ trees があったとする． このデータに関して散布図を描いた後，単回帰分析を行う．
Girth ： 木の直径（数値，単位はインチ）
Height ： 木の高さ(数値，単位はフィート）
Volume ： 長さ 1 フィートの立方体に切ったときの重さ
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() により，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) |
単回帰分析・直線回帰
