関数 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) |