とりあえずデータを読み込み、OLS推定を行いましょう。 “c:/documents/”のところはd01.csvが保存されているフォルダを指定してください。
> d01<-read.csv(“c:/documents/d01.csv”)
> o01<-lm(price~.-No,d01)
> summary(o01)
Call:
lm(formula = price ~ . – No, data = d01)Residuals:
Min 1Q Median 3Q Max
-37.5588 -7.0805 -0.2964 7.1464 39.7761Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.214e+03 2.615e+02 -23.765 < 2e-16 ***
areakyoto -7.367e-02 8.357e-01 -0.088 0.92977
areaosaka 7.763e-01 8.147e-01 0.953 0.34091
dealer -8.255e-01 6.675e-01 -1.237 0.21652
modelyear 3.138e+00 1.306e-01 24.031 < 2e-16 ***
gradeGX 4.330e+00 9.242e-01 4.685 3.20e-06 ***
gradeS -5.717e+00 9.413e-01 -6.074 1.79e-09 ***
gradeT -1.236e+00 9.383e-01 -1.317 0.18798
run -2.206e-04 1.457e-05 -15.140 < 2e-16 ***
displ 2.572e-02 2.047e-03 12.565 < 2e-16 ***
AT 5.050e+00 1.546e+00 3.265 0.00113 **
audio 4.382e+00 7.837e-01 5.591 2.93e-08 ***
navi 1.122e+01 7.262e-01 15.444 < 2e-16 ***
colorbrack 1.263e+01 1.228e+00 10.284 < 2e-16 ***
colorred 1.775e+01 1.822e+00 9.743 < 2e-16 ***
colorsilver 1.015e+01 1.202e+00 8.451 < 2e-16 ***
colorwhite 6.784e+00 1.200e+00 5.653 2.06e-08 ***
damage -5.402e+00 2.730e-01 -19.788 < 2e-16 ***
precord -1.758e+01 1.525e+00 -11.531 < 2e-16 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 10.47 on 981 degrees of freedom
Multiple R-squared: 0.6787, Adjusted R-squared: 0.6728
F-statistic: 115.1 on 18 and 981 DF, p-value: < 2.2e-16> AIC(o01)
[1] 7555.906
項目変数「area」はダミー変数として扱われていおり、「areakyoto」「areaosaka」となっている。「areahyogo」がないのは、これを基準にしているから。それぞれの値は、「areahyogo」に対して「areakyoto」と「areaosaka」は推定係数だけ高いということになります。「grade」「color」も同様。
ただし変数「area」と「dealer」の係数が有意でないので、これを削除してみます。
もう一度、lm(price~.-No-area-dealer,d01)としてもよいが、updateという関数で、一度行った推定を修正することもできます。
> o02<-update(o01,.~.-area-dealer)
> AIC(o02)
[1] 7552.933
ほんのちょっとAICが改善した。gradeも抜いてみると
> o03<-update(o02,.~.-grade)
> AIC(o03)
[1] 7654.08
AICは悪化。runを自然対数をとってみると
> o04<-update(o02,.~.-run+log(run))
> AIC(o04)
[1] 7473.465
AICがかなり改善。結局、結果は、以下のようになる。
> summary(o04)
Call:
lm(formula = price ~ modelyear + grade + displ + AT + audio +
navi + color + damage + precord + log(run), data = d01)
Residuals:
Min 1Q Median 3Q Max
-35.57573 -7.01326 -0.03811 6.72185 30.11925
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.192e+03 2.509e+02 -24.676 < 2e-16 ***
modelyear 3.154e+00 1.253e-01 25.177 < 2e-16 ***
gradeGX 4.118e+00 8.860e-01 4.648 3.81e-06 ***
gradeS -5.741e+00 9.013e-01 -6.370 2.89e-10 ***
gradeT -1.256e+00 9.008e-01 -1.395 0.163392
displ 2.521e-02 1.966e-03 12.824 < 2e-16 ***
AT 4.951e+00 1.483e+00 3.338 0.000875 ***
audio 4.332e+00 7.510e-01 5.769 1.07e-08 ***
navi 1.149e+01 6.978e-01 16.462 < 2e-16 ***
colorbrack 1.246e+01 1.179e+00 10.565 < 2e-16 ***
colorred 1.786e+01 1.750e+00 10.203 < 2e-16 ***
colorsilver 9.646e+00 1.153e+00 8.365 < 2e-16 ***
colorwhite 6.334e+00 1.151e+00 5.503 4.76e-08 ***
damage -5.322e+00 2.620e-01 -20.313 < 2e-16 ***
precord -1.794e+01 1.460e+00 -12.285 < 2e-16 ***
log(run) -6.122e+00 3.370e-01 -18.165 < 2e-16 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.06 on 984 degrees of freedom
Multiple R-squared: 0.7024, Adjusted R-squared: 0.6978
F-statistic: 154.8 on 15 and 984 DF, p-value: < 2.2e-16
新たな車を評価してみたい。新たな車のデータを（この場合1件だけだが）d01.csvと同じ形式で保存し、Rに読み込みます。
その後、predict関数を使って推定値を当てはめることができます。
> o05<-predict(o04,d01)
> d02<-read.csv(“c:/home/doc/seminar/Rlesson/d02.csv”,T)
> o06<-predict(o04,d02)
> o06
1
116.6253
新たな車の期待値は116万円。100万円というのは、「どちらかといえば買いかな」程度の安さです。グラフに描いてみるとよくわかります。
> plot(o05,d01$price,pch=20)
> abline(lm(d01$price~o05),col=4)
> points(o06,100,col=2,pch=19,cex=2)
> text(o06,100,”another car”,col=2,pos=4)
R Tips
いろんな変数を抜いたり足したりして最小のAICのモデルを探す（ステップワイズ法と言います）関数がRにはあります。stepという関数です。
すべての変数で回帰した結果がo01に入っているとき、step(o01)でできます。
> o04<-step(o01)
Start: AIC=4716.03
price ~ (No + area + dealer + modelyear + grade + run + displ +
AT + audio + navi + color + damage + precord ) – No
Df Sum of Sq RSS AIC
– area 2 151 107709 4713.4
– dealer 1 168 107726 4715.6
<none> 107558 4716.0
– AT 1 1169 108727 4724.8
– audio 1 3427 110985 4745.4
– grade 3 12393 119950 4819.1
– precord 1 14579 122136 4841.1
– color 4 17178 124736 4856.2
– displ 1 17311 124869 4863.3
– run 1 25130 132688 4924.0
– navi 1 26151 133709 4931.7
– damage 1 42933 150491 5049.9
– modelyear 1 63317 170875 5176.9
Step: AIC=4713.43
price ~ dealer + modelyear + grade + run + displ + AT + audio +
navi + color + damage + precord
Df Sum of Sq RSS AIC
– dealer 1 175 107884 4713.1
<none> 107709 4713.4
– AT 1 1192 108901 4722.4
– audio 1 3377 111086 4742.3
– grade 3 12292 120001 4815.5
– precord 1 14598 122307 4838.5
– color 4 17170 124880 4853.3
– displ 1 17334 125043 4860.7
– run 1 25105 132814 4920.9
– navi 1 26154 133863 4928.8
– damage 1 42791 150500 5046.0
– modelyear 1 63179 170888 5173.0
Step: AIC=4713.06
price ~ modelyear + grade + run + displ + AT + audio + navi +
color + damage + precord
Df Sum of Sq RSS AIC
<none> 107884 4713.1
– AT 1 1194 109078 4722.1
– audio 1 3328 111212 4741.4
– grade 3 12201 120085 4814.2
– precord 1 14435 122319 4836.6
– color 4 17224 125108 4853.2
– displ 1 17284 125168 4859.7
– run 1 25171 133055 4920.8
– navi 1 26255 134139 4928.9
– damage 1 42698 150582 5044.5
– modelyear 1 63035 170919 5171.2
> summary(o04)
Call:
lm(formula = price ~ modelyear + grade + run + displ + AT + audio +
navi + color + damage + precord, data = d01)
Residuals:
Min 1Q Median 3Q Max
-37.5356 -7.0990 -0.1771 6.9833 39.7709
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.192e+03 2.611e+02 -23.712 < 2e-16 ***
modelyear 3.127e+00 1.304e-01 23.978 < 2e-16 ***
gradeGX 4.294e+00 9.218e-01 4.659 3.62e-06 ***
gradeS -5.656e+00 9.383e-01 -6.028 2.35e-09 ***
gradeT -1.222e+00 9.374e-01 -1.304 0.19259
run -2.207e-04 1.457e-05 -15.152 < 2e-16 ***
displ 2.568e-02 2.046e-03 12.556 < 2e-16 ***
AT 5.093e+00 1.543e+00 3.300 0.00100 **
audio 4.306e+00 7.815e-01 5.509 4.60e-08 ***
navi 1.123e+01 7.260e-01 15.475 < 2e-16 ***
colorbrack 1.261e+01 1.227e+00 10.273 < 2e-16 ***
colorred 1.768e+01 1.821e+00 9.704 < 2e-16 ***
colorsilver 1.011e+01 1.199e+00 8.432 < 2e-16 ***
colorwhite 6.668e+00 1.197e+00 5.568 3.32e-08 ***
damage -5.379e+00 2.726e-01 -19.734 < 2e-16 ***
precord -1.743e+01 1.519e+00 -11.474 < 2e-16 ***
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 10.47 on 984 degrees of freedom
Multiple R-squared: 0.6777, Adjusted R-squared: 0.6728
F-statistic: 138 on 15 and 984 DF, p-value: < 2.2e-16
結局、手動で求めたo02と同じモデルが得られました。しかしながら、特定の変数の対数をとるなどの工夫はやってくれませんので要注意。