Rで、最小二乗法は lm(y ~ x1+x2) と指定します。
なぜコマンド名は lm なんですか？
linear models （線形モデル）の略です。
この記事では、Rを使った簡単な最小二乗法の推定法を説明します。lm関数を使うと最小二乗法が推定でき、summary関数で結果が見られます。predict関数を使うと予測ができます。
経済統計の使い方では、統計データの入手法から分析法まで解説しています。
データフレームの作成
cp <- data.frame("GDP95"=c(311988.0 ,320726.7 ,330912.0 ,338991.9 ,352924.7 ,368211.9 ,379845.5 ,398927.2 ,424294.1 ,444876.0 ,469567.2 ,480855.2 ,483022.6 ,485298.3 ,489588.5 ,504827.3 ,521364.8 ,522220.9 ,518706.4 ,520765.6 ,539161.1 ,532444.5 ,541200.0 ), "CP95"=c(174234.2 ,176912.6 ,184665.6 ,189188.9 ,194156.2 ,201565.5 ,208975.1 ,217229.6 ,229047.3 ,238716.3 ,249355.1 ,257141.5 ,261717.7 ,266382.5 ,272215.1 ,279472.3 ,286514.5 ,282965.3 ,286104.1 ,288346.2 ,290585.4 ,294798.9 ,298974.3 ))実質GDPと実質民間最終消費支出のデータを入力し、cpというデータフレームに格納しました。
matplotコマンドでグラフを描くと以下のようになります。1が実質GDP,２が実質民間最終消費です。
matplot(cp,type="l")
最小二乗法
最小二乗法のコマンドは lm です。被説明変数を最初に置き、～マークの次に説明変数を置きます。2つ以上説明変数がある場合は、＋記号を使って並べます。 y ~ x1 + x2 のような形です。
被説明変数がCP95、説明変数がGDP95の場合は以下の式になります。推計結果はresultという名前にしていますが、名前は任意に決められます。
$ result <- lm(GP95 ~ GDP95,data=cp) $
結果を見るには、summary関数を使います。summary(result)で以下の結果が表示されます。
summary(result)
Call:
lm(formula = CP95 ~ GDP95, data = cp)
Residuals:
Min 1Q Median 3Q Max
-7565 -2393 411 2672 4505
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.845e+03 4.316e+03 0.891 0.383
GDP95 5.390e-01 9.514e-03 56.646 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3531 on 21 degrees of freedom
Multiple R-squared: 0.9935, Adjusted R-squared: 0.9932
F-statistic: 3209 on 1 and 21 DF, p-value: < 2.2e-16最小二乗法の推定結果を検討するために必要な、係数、ｔ値、p値、Ｆ値、自由度修正済み決定係数などが入手できます。推定された式は以下です。
$ CP95＝3.845e+03+5.390e-01×GDP95 +e $
指数表示を普通の表示に直すと以下の通りです。
$ CP95＝3845+0.5390×GDP95 +e $
定数項である3845に対するt値は0.89でp値は0.383、GDP95に係る係数の0.539に対するt値は56.646で、ｐ値は2e-16（ほぼゼロ）です。F値は自由度（1，21）で3209で対応するp値は 2.2e-16になります。自由度修正済み決定係数は0.9932となります。
そのほかにも、様々な出力ができます。推定結果をresultとした場合の表です。
| 関数 | 内容 | 例 |
| fitted | 推定値 | result$fitted |
| predict | 予測値 | result$predict,fitted(result) |
| 簡単な結果の要約 | print(result) | |
| summary | 結果の要約 | summary(result) |
| plot | 回帰診断 | plot(result) |
推計値は、result$fittedで得ることが出来ます。fitted(result)とすると、推定値が表示されます。
fitted(result)
1 2 3 4 5 6 7 8
171992.1 176701.9 182191.3 186546.0 194055.1 202294.1 208564.1 218848.2
9 10 11 12 13 14 15 16
232519.8 243612.4 256919.8 263003.5 264171.6 265398.1 267710.3 275923.3
17 18 19 20 21 22 23
284836.3 285297.6 283403.5 284513.3 294427.6 290807.7 295526.5 printコマンドでは、式の形と係数が表示されます。
> print(result)
Call:
lm(formula = CP95 ~ GDP95, data = cp)
Coefficients:
(Intercept) GDP95
3845.475 0.539 回帰診断は、plot(result)で実行できます。4種類のグラフがリターンキーを押すと次々に出てきます。
> plot(result)
Hit <Return> to see next plot:
Hit <Return> to see next plot:
Hit <Return> to see next plot:
Hit <Return> to see next plot: 
予測
次に予測の方法を説明します。推定した方程式が以下の場合で説明します。
$ CP95＝3845+0.5390×GDP95 +e $
右辺のGDP95に新たなデータを入れ、誤差項はゼロとすると、CP95の予測値が得られます。最小二乗法による予測は、右辺にあるGDP95 はあらかじめ何らかの形で予測しておかないとできません。
現在データは23期目まであるので、24期目のGDP95が3％増えたとして、それに対応したCP95を予測します。まず新たなデータフレームとしてcpをコピーしてcpfを作成し、24期目のデータを23期目の３％増加したものとします。
> cpf <- cp
> cpf[24,1] <- cpf[23,1]*1.03この新たなデータに予測するためのpredcit関数を適用すると予測ができます。lm関数の時はdata＝としてデータフレームを指定しましたが、predict関数ではnewdata=とするところがポイントです。予測値をCP95fという名前でデータフレームcpfに格納します。
cpf$CP95f <- predict(result,newdata=cpf）予測値を含めたデータが以下のように作成されました。23期より前は推定値、24期が予測値になります。
> cpf
GDP95 CP95 CP95f
1 311988.0 174234.2 171992.1
2 320726.7 176912.6 176701.9
3 330912.0 184665.6 182191.3
4 338991.9 189188.9 186546.0
5 352924.7 194156.2 194055.1
6 368211.9 201565.5 202294.1
7 379845.5 208975.1 208564.1
8 398927.2 217229.6 218848.2
9 424294.1 229047.3 232519.8
10 444876.0 238716.3 243612.4
11 469567.2 249355.1 256919.8
12 480855.2 257141.5 263003.5
13 483022.6 261717.7 264171.6
14 485298.3 266382.5 265398.1
15 489588.5 272215.1 267710.3
16 504827.3 279472.3 275923.3
17 521364.8 286514.5 284836.3
18 522220.9 282965.3 285297.6
19 518706.4 286104.1 283403.5
20 520765.6 288346.2 284513.3
21 539161.1 290585.4 294427.6
22 532444.5 294798.9 290807.7
23 541200.0 298974.3 295526.5
24 557436.0 NA 304276.9実績値と予測値のグラフを描くと以下のようになります。実線が実績値、点線が推定値（実績がある23期まで）と予測値（24期）です。
matplot(cpf[,2:3],type="l")
最小二乗法による係数の求め方
lm関数を使えば、係数は簡単に計算されます。実際にどのような手順で作成されているかを考えると以下の式になります。重回帰の分析の場合にも使えます。
$y$は被説明変数のベクトル、$X$は、定数項を表す要素が1の列と、説明変数の列を合わせた行列です。説明変数が1つの場合は2列で、説明変数が増えると列が増えます。
行列の計算式として、Xを例にとると逆行列はsolve(X)、行列の積は%*％、行と列を入れ替えた転置行列 $X’$ はt (X)で表します。
$ \beta = (X’ X)^{-1} X’ y $
> y <- cp$CP95
> x <- cbind(1, cp$GDP95)
> solve( t(x) %*% x ) %*% t(x) %*% y
[,1]
[1,] 3845.4748610
[2,] 0.5389523 計算の結果、定数項が3845.4748610、GDP95 に係る係数が0.5389523と計算されます。
