パッケージのインストール

# MASSパッケージをインストール
install.packages("MASS",repos = "http://cran.us.r-project.org")
## 
## The downloaded binary packages are in
##  /var/folders/nm/m79f0x2d79768ww18g9h71340000gp/T//Rtmpbh0yTT/downloaded_packages

データの用意

# ライブラリを読み込む
library(MASS)
# ボストン住宅データセットを読み込む
data(Boston)
# データセットの内容を表示
head(Boston)
##      crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
## 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
## 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
## 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
## 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
## 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
## 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
##   medv
## 1 24.0
## 2 21.6
## 3 34.7
## 4 33.4
## 5 36.2
## 6 28.7

データセットの列の変数は次のような意味。 crim: 町ごとの一人当たりの犯罪率

zn: 25,000平方フィート以上の住宅区画の割合

indus: 非小売業の土地面積の割合

chas: チャールズ川の周囲にあるかどうか(1: 川の境界に接する、0: それ以外)

nox: 窒素酸化物の濃度

rm: 住居ごとの平均部屋数

age: 1940年より前に建てられた所有者居住ユニットの割合

dis: 5つのボストンの雇用センターまでの重み付き距離

rad: ラジアルハイウェイへのアクセス性の指標

tax: 10,000ドルあたりの固定資産税率

ptratio: 町ごとの生徒と教師の比率

black: 町ごとの黒人の比率を示す指標

lstat: 低所得者の割合

medv: 住宅価格

相関係数を求め散布図を描く

cor(Boston)
##                crim          zn       indus         chas         nox
## crim     1.00000000 -0.20046922  0.40658341 -0.055891582  0.42097171
## zn      -0.20046922  1.00000000 -0.53382819 -0.042696719 -0.51660371
## indus    0.40658341 -0.53382819  1.00000000  0.062938027  0.76365145
## chas    -0.05589158 -0.04269672  0.06293803  1.000000000  0.09120281
## nox      0.42097171 -0.51660371  0.76365145  0.091202807  1.00000000
## rm      -0.21924670  0.31199059 -0.39167585  0.091251225 -0.30218819
## age      0.35273425 -0.56953734  0.64477851  0.086517774  0.73147010
## dis     -0.37967009  0.66440822 -0.70802699 -0.099175780 -0.76923011
## rad      0.62550515 -0.31194783  0.59512927 -0.007368241  0.61144056
## tax      0.58276431 -0.31456332  0.72076018 -0.035586518  0.66802320
## ptratio  0.28994558 -0.39167855  0.38324756 -0.121515174  0.18893268
## black   -0.38506394  0.17552032 -0.35697654  0.048788485 -0.38005064
## lstat    0.45562148 -0.41299457  0.60379972 -0.053929298  0.59087892
## medv    -0.38830461  0.36044534 -0.48372516  0.175260177 -0.42732077
##                  rm         age         dis          rad         tax    ptratio
## crim    -0.21924670  0.35273425 -0.37967009  0.625505145  0.58276431  0.2899456
## zn       0.31199059 -0.56953734  0.66440822 -0.311947826 -0.31456332 -0.3916785
## indus   -0.39167585  0.64477851 -0.70802699  0.595129275  0.72076018  0.3832476
## chas     0.09125123  0.08651777 -0.09917578 -0.007368241 -0.03558652 -0.1215152
## nox     -0.30218819  0.73147010 -0.76923011  0.611440563  0.66802320  0.1889327
## rm       1.00000000 -0.24026493  0.20524621 -0.209846668 -0.29204783 -0.3555015
## age     -0.24026493  1.00000000 -0.74788054  0.456022452  0.50645559  0.2615150
## dis      0.20524621 -0.74788054  1.00000000 -0.494587930 -0.53443158 -0.2324705
## rad     -0.20984667  0.45602245 -0.49458793  1.000000000  0.91022819  0.4647412
## tax     -0.29204783  0.50645559 -0.53443158  0.910228189  1.00000000  0.4608530
## ptratio -0.35550149  0.26151501 -0.23247054  0.464741179  0.46085304  1.0000000
## black    0.12806864 -0.27353398  0.29151167 -0.444412816 -0.44180801 -0.1773833
## lstat   -0.61380827  0.60233853 -0.49699583  0.488676335  0.54399341  0.3740443
## medv     0.69535995 -0.37695457  0.24992873 -0.381626231 -0.46853593 -0.5077867
##               black      lstat       medv
## crim    -0.38506394  0.4556215 -0.3883046
## zn       0.17552032 -0.4129946  0.3604453
## indus   -0.35697654  0.6037997 -0.4837252
## chas     0.04878848 -0.0539293  0.1752602
## nox     -0.38005064  0.5908789 -0.4273208
## rm       0.12806864 -0.6138083  0.6953599
## age     -0.27353398  0.6023385 -0.3769546
## dis      0.29151167 -0.4969958  0.2499287
## rad     -0.44441282  0.4886763 -0.3816262
## tax     -0.44180801  0.5439934 -0.4685359
## ptratio -0.17738330  0.3740443 -0.5077867
## black    1.00000000 -0.3660869  0.3334608
## lstat   -0.36608690  1.0000000 -0.7376627
## medv     0.33346082 -0.7376627  1.0000000
pairs(Boston) #全部

ひとまず、いくつかの変数は相関が高く、いくつかはそれほど高くないように見える。説明変数同士の相関が高い場合、多重共線性を考慮して、どちらかの変数を抜いた上で、分析を行う必要が出てくる。 ひとまず、厳密な解釈のための変数選択は後回しにするとして、重回帰分析の手法としては、下記の命令で実行できる。

分析

lmfit = lm(medv ~ ., data=Boston)
summary(lmfit)
## 
## Call:
## lm(formula = medv ~ ., data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.595  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
## crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
## zn           4.642e-02  1.373e-02   3.382 0.000778 ***
## indus        2.056e-02  6.150e-02   0.334 0.738288    
## chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
## nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
## rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
## age          6.922e-04  1.321e-02   0.052 0.958229    
## dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
## rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
## tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
## ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
## black        9.312e-03  2.686e-03   3.467 0.000573 ***
## lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
## F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

medv ~ .とすることで、目的変数をmedv, 説明変数はそれ以外の全てを意味する。変数が多い場合に便利。

lmfit_2 = lm(medv ~ crim+zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstat, data=Boston)
summary(lmfit_2)
## 
## Call:
## lm(formula = medv ~ crim + zn + indus + chas + nox + rm + age + 
##     dis + rad + tax + ptratio + black + lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.595  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
## crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
## zn           4.642e-02  1.373e-02   3.382 0.000778 ***
## indus        2.056e-02  6.150e-02   0.334 0.738288    
## chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
## nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
## rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
## age          6.922e-04  1.321e-02   0.052 0.958229    
## dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
## rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
## tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
## ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
## black        9.312e-03  2.686e-03   3.467 0.000573 ***
## lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
## F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

medv ~ crim+zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstatとしても同じ意味。 説明変数を+で繋げていく。

回帰の診断

このユニットを作成にあたり https://htsuda.net/stats/regression.html#%E5%9B%9E%E5%B8%B0%E3%81%AE%E8%A8%BA%E6%96%AD%E6%96%B9%E6%B3%95 を参照しています。

par(mfrow=c(2,2))
plot(lmfit)

正規性(Normality)

目的変数が正規分布に従うという条件。 もしそうなっていれば、Normal Q-Qにおいてデータは対角線上に分布 値が小さい場合と大きい場合は少し外れる結果。

residplot = function( lmfit, nbreaks = 10 ){
    z = rstudent( lmfit )
    hist( z, breaks = nbreaks, freq = F,
        xlab = "Studentized Residual", main = "Distribution of Errors" )
    rug( jitter( z ), col = "brown" )
    curve( dnorm( x, mean = mean( z ), sd = sd( z ) ), add = T, col = "blue", lwd = 2 )
    lines( density( z )$x, density( z )$y, col = "red", lwd = 2, lty = 2 )
    legend( "topleft", legend = c( "Normal Curve", "Kernel Density Curve"),
        lty = 1:2, col = c( "blue", "red" ), cex =.7 )
}

residplot( lmfit )

正規分布と実際データの確率密度関数の違いがよくわかる。

独立性(Independence)

各観測が独立でなければならないという条件。 時系列データや階層性のあるデータの場合にはこの条件が満たされない。 car パッケージの durbinWatsonTest() 関数は Durbin-Watson 検定によって系列的な相関について調べることができる。

car::durbinWatsonTest( lmfit )
##  lag Autocorrelation D-W Statistic p-value
##    1       0.4542626      1.078375       0
##  Alternative hypothesis: rho != 0
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.04687792      2.083648    0.59
##  Alternative hypothesis: rho != 0

p < 0.05 であるため自己相関がある(つまり変数選択しなけばならない)

線形性(Linearity)

目的変数が説明変数と線形な関係を持つという条件。 Component plus residual プロットを使うと線形性を確認できる。

par( mar = c( 4.5, 4.5, 3, 1 ) )
#マージンは、プロットエリア(グラフの描画領域)と図全体との間のスペースを指定。par(mar = c(4.5, 4.5, 3, 1))の引数c(4.5, 4.5, 3, 1)は、マージンの上部、下部、左側、右側のサイズを指定。
library(car)
## Loading required package: carData
crPlots( lmfit )

2つの線が重なっていれば、線形性が満たされる。 crim, chas, rm は重なっていない。

等分散性(Homoscedasticity)

目的変数の分散が位置によらず一定であるという条件。 もし等分散性が満たされていれば残差の値は予測値の値に依存しない(つまり水平な分布)になる

library(car)
ncvTest( lmfit )
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 15.24407, Df = 1, p = 9.4473e-05
par( mar = c( 4.5, 4.5, 2, 1 ) )
spreadLevelPlot( lmfit )
## Warning in spreadLevelPlot.lm(lmfit): 
## 1 negative fitted value removed

## 
## Suggested power transformation:  0.8541632

p 値が 0.05 未満であり、等分散ではないという結果

総当たり法でAICを求める。

http://stat.sys.i.kyoto-u.ac.jp/titech/class/dk2005/san08.pdfを参照しています

library(leaps)
best<-regsubsets(medv~., data = Boston, nvmax = 13, method = "seqrep")
sbest<-summary(best)
plot(best,scale="adjr2")

sbest$rss #これらのモデルのRSS(残差変動)
##  [1] 19472.38 15439.31 13727.99 13228.91 12469.34 12141.07 11868.24 11678.30
##  [9] 11526.12 11308.58 14082.96 13489.62 11078.78
n<-ncol(Boston) #Bostonの変数の数をnとする
p<-nrow(sbest$which) #ベストモデルsbestの数
aicbest<- n*(1+log(2*pi*sbest$rss/n))+2*(2+(1:p)) #各モデルのAICの計算
aicbest #AICの表示
##  [1] 147.0580 145.8089 146.1642 147.6457 148.8179 150.4444 152.1262 153.9003
##  [9] 155.7167 157.4499 162.5216 163.9189 163.1625
order(aicbest) #AICの小さい順に並べると2番目のモデルが最小
##  [1]  2  3  1  4  5  6  7  8  9 10 11 13 12
sbest$which[2,] #2番目のモデルの変数の有無を表示。TRUEが採用する変数。
## (Intercept)        crim          zn       indus        chas         nox 
##        TRUE       FALSE       FALSE       FALSE       FALSE       FALSE 
##          rm         age         dis         rad         tax     ptratio 
##        TRUE       FALSE       FALSE       FALSE       FALSE       FALSE 
##       black       lstat 
##       FALSE        TRUE
aicbest[2] #2番目のモデルのAIC
## [1] 145.8089

medv~rm+lstatのモデルが最良という結果になった。 そのAICは145.8089

多重共線性の評価

bestfit<-lm(medv~rm+lstat,data=Boston)
summary(bestfit)
## 
## Call:
## lm(formula = medv ~ rm + lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.076  -3.516  -1.010   1.909  28.131 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.35827    3.17283  -0.428    0.669    
## rm           5.09479    0.44447  11.463   <2e-16 ***
## lstat       -0.64236    0.04373 -14.689   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.54 on 503 degrees of freedom
## Multiple R-squared:  0.6386, Adjusted R-squared:  0.6371 
## F-statistic: 444.3 on 2 and 503 DF,  p-value: < 2.2e-16
library(car)
vif(bestfit)
##      rm   lstat 
## 1.60452 1.60452
cor(Boston$rm, Boston$lstat)
## [1] -0.6138083
plot(Boston$rm, Boston$lstat)

rm 1.60452 lstat 1.60452 となり、10を超えていないので多重共線性は見られない。

ただし、散布図と相関係数を見ると、-0.6程度で中程度の負の相関が見られる。

rm: 住居ごとの平均部屋数

lstat: 低所得者の割合

なので、住宅価格= 5.09479×住居ごとの平均部屋数-0.64236×低所得者の割合-1.35827

部屋数が増えれば住宅価格が上がり、低所得者の割合が増えれば、住宅価格が下がるという関係が見られる。

標準化回帰係数

z = as.data.frame( scale(Boston) ) #scale関数で各列の値は平均0、標準偏差1の分布となる
zfit = lm( medv ~ rm + lstat, data = z )
options( scipen = 20 ) #表示される数値の指数表記の桁数(表示精度)を設定
summary( zfit )
## 
## Call:
## lm(formula = medv ~ rm + lstat, data = z)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.9654 -0.3823 -0.1098  0.2075  3.0587 
## 
## Coefficients:
##                         Estimate           Std. Error t value
## (Intercept)  0.00000000000000129  0.02677956880524802    0.00
## rm           0.38921875253899446  0.03395515341075680   11.46
## lstat       -0.49875703630006041  0.03395515341075676  -14.69
##                        Pr(>|t|)    
## (Intercept)                   1    
## rm          <0.0000000000000002 ***
## lstat       <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6024 on 503 degrees of freedom
## Multiple R-squared:  0.6386, Adjusted R-squared:  0.6371 
## F-statistic: 444.3 on 2 and 503 DF,  p-value: < 0.00000000000000022

住宅価格= 0.389218752538994461343 ×住居ごとの平均部屋数-0.498757036300060407896 ×低所得者の割合+ 0.000000000000001289686

標準化することで、低所得者層の割合の方が平均部屋数より住宅価格にマイナスの影響を与えていることがわかる。

その他、各検定の結果は有意水準5%の棄却域にあるので、検定の結果は有意である。

その他注意点

このデータでは線形性と等分散性が見られない部分があった。 このような場合、説明変数の選択を行うほか、データの外れ値を抜いて回帰、目的変数の指数変換や対数変換を行う必要がある。 授業で取り扱うレベルを超えるので、下記に参考資料を示します。 https://htsuda.net/stats/regression.html#%E5%AF%BE%E5%87%A6%E6%96%B9%E6%B3%95