axjack's blog

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.

重回帰分析をRで

やること

永田『多変量解析法入門』(以下参考書)よりp.2のデータをもとに重回帰分析を行う。

ソースコード

# データ入力 #### 
hirosa <- c(51,38,57,51,53,77,63,69,72,73)
tiku   <- c(16,4,16,11,4,22,5,5,2,1)
kakaku <- c(3,3.2,3.3,3.9,4.4,4.5,4.5,5.4,5.4,6)

# データフレーム #### 
d <- data.frame(
  hirosa
  ,tiku
  ,kakaku
)

# 重回帰式 #### 
d.lm <- lm(kakaku ~ ., data = d)

# 要約 ####
summary(d.lm)

# 係数のみ表示 ####
print(d.lm)

# 新しいデータが与えられた時の予測値 #### 
# 新しいデータをデータフレームで用意し、
givenData <- data.frame(hirosa=70,tiku=10)

# 予測区間を出す
predict(d.lm,newdata=givenData,interval = "predict")

結果の確認

(1)回帰式

以下の通り参考書と一致。

> print(d.lm)

Call:
lm(formula = kakaku ~ ., data = d)

Coefficients:
(Intercept)       hirosa         tiku  
    1.02013      0.06680     -0.08083  

(2)自由度調整済み寄与率

Adjusted R-squared: 0.9336 となった。

> summary(d.lm)

Call:
lm(formula = kakaku ~ ., data = d)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.32468 -0.20952  0.03939  0.17150  0.36196 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.020130   0.443624   2.300 0.055029 .  
hirosa       0.066805   0.007065   9.456 3.09e-05 ***
tiku        -0.080830   0.012241  -6.603 0.000303 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2636 on 7 degrees of freedom
Multiple R-squared:  0.9484,    Adjusted R-squared:  0.9336 
F-statistic:  64.3 on 2 and 7 DF,  p-value: 3.126e-05

(3)同じ地区で広さ=70平米, 築年数=10年, 価格=5.8千万円の提示は妥当か

係数を使って線形式に与えられたデータを代入し価格の予測値を出しても良い。

> coef(d.lm)
(Intercept)      hirosa        tiku 
 1.02012955  0.06680477 -0.08082993 

> sum(coef(d.lm) * c(1,70,10) )
[1] 4.888164

が、predict関数を用いてみると、

# 新しいデータをデータフレームで用意し、
givenData <- data.frame(hirosa=70,tiku=10)

> # 予測区間を出す
> predict(d.lm,newdata=givenData,interval = "predict")
       fit      lwr      upr
1 4.888164 4.214119 5.562209

と、予測値と予測区間が出るので便利。参考書と値が一致していることが確認できた。