# 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)
目的変数が正規分布に従うという条件。 もしそうなっていれば、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 )
正規分布と実際データの確率密度関数の違いがよくわかる。
各観測が独立でなければならないという条件。 時系列データや階層性のあるデータの場合にはこの条件が満たされない。 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 であるため自己相関がある(つまり変数選択しなけばならない)
目的変数が説明変数と線形な関係を持つという条件。 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 は重なっていない。
目的変数の分散が位置によらず一定であるという条件。 もし等分散性が満たされていれば残差の値は予測値の値に依存しない(つまり水平な分布)になる
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 未満であり、等分散ではないという結果
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