第157回： 「統計の実務」17 回帰分析《その3：重回帰分析 後編》
◀前の記事へ 次の記事へ▶
≡ はじめに
前回は、「回帰分析《その2：重回帰分析 前編》」について書きました。
内容としては、「複数の説明変数と目的変数との因果関係が直線で表せるときに、その回帰線の各説明変数の回帰係数とY軸の切片をRコマンダーで計算しよう」という話でした。最後の方で「マルチコ」（説明変数間に相関関係があると結果が正しく求まらない）に注意しよう！とも書きました。
今回は「重回帰分析 後編」です。
マーケットリサーチなどで大活躍している「どの説明変数が結果に効いているのかを調べる方法」について書きます。
様々な要因（立地条件、価格、接客態度、フロア面積、、、等々）が、結果（店舗売上等）に影響を与えているときに、どの要因の改善にリソースを集中するかの戦略立案はとても大切なことです。
そして、今回の方法を知っていると、どの要因の改善に集中すると高い効果が見込めるかの判断ができるようになります。
（もちろん良い元データがあればの話です。“Garbage In, Garbage Out ”です。☺️）
前回、
次回も本当なら「ラーメン屋の売上データ、要因付き」を使うと良いとは思うのですが、たぶん見つからないので、今回のように手順を示すところまでとなりそうです。
と書きました。
その後、何度もググって良いデータを探しました。けれど、見つかりませんでした。
今回、前回のように「フィッシャーのあやめ」で、手順を示すことだけをするか、封印していた「Boston Housing」を使うかでとても悩みました。
手順だけを示されても面白くないと思うからです。
また、重回帰分析では回帰線を求めることよりも、「どの説明変数が効いているのかを調べる」ほうが役に立つと思っています。
そこで、思い切って「Boston Housing」の封印を解くことにしました。大袈裟ですね。(笑)
ということで、今回のnoteでは、「Boston Housing」データを使って、重回帰分析によって、どの説明変数の寄与が大きいかを調べる方法について説明します。
≡ Boston Housing
“Boston Housing”は、ボストンの住宅価格という目的変数と、13個の説明変数（ほとんどが数値変数だけれど、一つだけ「チャールズ川の周辺かどうか」を0と1で示しているカテゴリ変数があります）からなるデータの集まりです。
データは506行あり、それぞれはボストンの町（タウン）のデータです。
つまり、「Boston Housing」は、ボストンの506の町ごとに、川の周辺かどうか、部屋数の平均値はいくつか、犯罪率はいくつかなどの13の説明変数と、住宅価格（その町の住宅価格の中央値です。平均値ではありません）がセットになったデータセットというわけです。
もともとは、大気汚染と住宅価格の調査で始まったので、NOX値（窒素酸化物濃度（1,000万分の1））の情報もあります。😀
ボストンというと、みなさんは何を思い浮かべますか？
私はボストンへは行ったことがないので「ボストンバッグ」と「ボストンクラブ」（＝ スタンハンセンがマーセルを倒した、逆エビ固め）を最初に思い浮かべました。続いて、「ボストン美術館」と「ボストンテリア」（犬）と「ボストン・コンサルティング」を思い出しました。
そういえば、「ボストンマラソン」も有名ですね。書いているうちに思い出しました。
あとはー。うーん🤔と頭をフル回転させて、「ボストン茶会事件」について歴史の授業で習ったことを捻り出しました。
そうそう、ボストンバッグは、ボストン大学の学生が使っていたことから名前がついた収納が大きなバッグ（数泊の旅行なら足りてしまうカバン）のことです。
で、ボストンバッグは和製英語のようです。
travelling bag、または、club bagといわないと欧米では通じないとか……。
ちなみに今回のキャッチイメージはボストンシティです。
閑話休題（無駄話はさておき）、Boston Housingに話を戻します。
Boston Housingは、13の説明変数が住宅価格にどれだけ寄与しているかを重回帰分析するのにピッタリなデータです。
元データ数は506個ですので、機械学習用のデータとしては少ないですが、重回帰分析の元データ数としては十分です。また、調査会社が行ったデータであり、、、一部精度が粗い変数や上限を丸めすぎな変数もありますが、、、そこそこ信用がおけます。
ええと、良い子は全ての変数、それぞれについてヒストグラムや箱ひげ図を描いて眺めると良いですよ。
統計解析の前に元データの分布を色々グラフ化して調べておくことは、滅茶苦茶大切です。
最初のうちは、『そんなこと面倒でやりたくないなあ』って思うんだけど、超おすすめです。
そして、パブリックドメインで自由に使える点もありがたいです。RではMASSパッケージをロードすると使えます。具体的には、メニューから、
[ツール]＞[パッケージのロード…]
を選択して出た下記のウィンドウでMASSを選択します。
MASSパッケージをロードできたら、メニューから、
[データ]＞[パッケージ内のデータ]＞[アタッチされたパッケージからデータを読み込む]
を選択して開いた以下の画面で、MASSをダブルクリック（シングルクリックでは選択できない。不思議）して、その右に現れたBostonを選べばOKです。
手順が煩雑なのと列名が英語の略語で、重回帰分析をした結果がピンとこないので、列名を日本語に直したCSVファイルを置いておきます。
≡ Boston Housingの問題点
ここまで読んで、「いいことずくめじゃないか！」と思われたと思います。私も大昔に初めて知ったときにはそう思いました。
問題は説明変数の一つ「B」にあります。
「DIS」も「雇用センターまでの距離」の意味を考えるとアレだけど……。
13の説明変数と目的変数の住宅価格について、原文とともに訳してみます。
CRIM - per capita crime rate by town
犯罪率 - 町ごとの一人当たりの犯罪発生率
ZN - proportion of residential land zoned for lots over 25,000 sq.ft.
豪邸率 - 25,000平方フィート（≒23アール）以上の区画に分類される住宅地の割合
INDUS - proportion of non-retail business acres per town.
商業 - 町ごとの小売業以外の商業が占める土地の割合。
CHAS - Charles River dummy variable (1 if tract bounds river; 0 otherwise)
川の周辺 - チャールズ川のダミー変数（1: 川の周辺、0: それ以外）。
NOX - nitric oxides concentration (parts per 10 million)
NOX - 窒素酸化物濃度（1,000万分の1）
RM - average number of rooms per dwelling
部屋数 - 住居の平均部屋数
AGE - proportion of owner-occupied units built prior to 1940
古い町並み - 1940年以前に建てられた持ち家の割合
DIS - weighted distances to five Boston employment centres
雇用センターまでの距離 - ボストンの5つの雇用センターまでの距離を加重平均したもの
RAD - index of accessibility to radial highways
高速道路 - 放射状高速道路へのアクセス指数
TAX - full-value property-tax rate per $10,000
固定資産税 - 1万ドルあたりの固定資産税全額負担率
PTRATIO - pupil-teacher ratio by town
教師比率 - 町ごとの生徒数-教師数比率
B - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
B - 1000(Bk - 0.63)^2 ただし、Bkは町ごとの黒人の割合
LSTAT - % lower status of the population
低所得者の割合 - 低い地位の人口比率
MEDV - Median value of owner-occupied homes in $1000's
住宅価格 - 1000ドル単位の持ち家の中央値
で、問題の「B」について再掲します。
B - 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
B - 1000(Bk - 0.63)^2 ただし、Bkは町ごとの黒人の割合
「1000(Bk - 0.63)^2」ってなんだよ？ と思いますよね？
なんでも、0.63は「理想的な黒人比率」だそうで、その町の黒人の割合との差を二乗して1000倍しているとのことです。
仮に黒人が1人もいなければ「1000(0 - 0.63)^2 = 396.9」です。
B変数のデータについて、ヒストグラムにするとこうなります。
なんといいますか、こういうデータを取っていること自体がありえないといいますか、腹立たしく、、、でずっと「こんなデータ使うもんか」と封印していたのです。
今回、この列を消して紹介しようか？とも思ったのですが、そんなことをするよりも、そのまま知っていることをお伝えして判断はみなさんにお任せするほうが良いかなと思い直しました。
以上で、「Boston Housing」データの話はお終いです。以降は、重回帰分析による説明変数の寄与率の求め方（テクニカルな話）となります。
≡ 重回帰分析による説明変数の寄与の調べかた
まずは、前回と全く同じ手順で重回帰分析をします。
準備として、Rコマンダーから、本noteの上の方においた「ボストン（MASSパッケージより）.csv」を読み込んでおきます。
次に、Rコマンダーのメニューから、
[統計量]＞[モデルへの適合]＞[線形回帰]
を選択します。
すると以下のウィンドウが開きますので、目的変数に「住宅価格」、説明変数に「住宅価格」と「No」以外をすべて選びます。
（住宅価格は目的変数で、Noは単なるデータ番号なので選びません。つまり、説明変数をすべて選びます。リストから複数のアイテムを選択するのは、やりにくい操作だけれど、シフトキーとコントロールキーを駆使して頑張ってください）
[OK]ボタンを押すと結果が出力されます。
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) 36.4594884 5.1034588 7.144 3.28e-12 ***
B 0.0093117 0.0026860 3.467 0.000573 ***
NOX -17.7666112 3.8197437 -4.651 4.25e-06 ***
教師比率 -0.9527472 0.1308268 -7.283 1.31e-12 ***
古い町並み 0.0006922 0.0132098 0.052 0.958229
固定資産税 -0.0123346 0.0037605 -3.280 0.001112 **
雇用センターまでの距離 -1.4755668 0.1994547 -7.398 6.01e-13 ***
高速道路 0.3060495 0.0663464 4.613 5.07e-06 ***
豪邸率 0.0464205 0.0137275 3.382 0.000778 ***
商業 0.0205586 0.0614957 0.334 0.738288
川の周辺 2.6867338 0.8615798 3.118 0.001925 **
低所得者の割合 -0.5247584 0.0507153 -10.347 < 2e-16 ***
犯罪率 -0.1080114 0.0328650 -3.287 0.001087 **
部屋数 3.8098652 0.4179253 9.116 < 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読み方は、前回と同じで、(Intercept)がY軸との切片で、あとは、それぞれの説明変数の回帰係数（今回は重回帰分析なので、偏回帰係数と呼ぶほうが正しい）です。
ここまでは前回と同じことしかしていません。
回帰式を書けば、
住宅価格 ＝ 36.461352 ＋ 0.009357×B ー 17.7666112×NOX …
です。
今回、知りたいことは、住宅価格に影響を与えている説明変数の影響度合いです。
さて、回帰係数の絶対値が一番大きな「NOX」が一番大きく関係するのか？というと、そんなことはありません。
上の表は各変数の平均値ですが、「NOX」はその平均値が、0.６と「川の周辺」の0.1（ただし、これはカテゴリ変数）の次に小さな説明変数です。小さい値だから少しの変化で結果に影響を与えるので回帰係数は大きくなります。
これで回帰係数をそのまま影響度合いの分析に使えないことがわかりました。
どうすればよいかというと、「変数の標準化」を行えば良いです。変数の標準化とは、変数の平均値と標準偏差を共通の基準に変換することで、データの単位を揃えることをイメージしてもらっても構いません。
まず、Rコマンダーのメニューから、
[データ]＞[アクティブデータセット内の変数の管理＞[変数の標準化…]
を選択します。
すると以下のウィンドウが現れますので、「No」以外の変数をすべて選んで[OK]します。
先ほどと同じく、変数の標準化がされた各変数の平均値を出してみます。
今度は桁があっていることが分かります。これで比較可能になりました。
イメージとしては試験の得点を偏差値に変えてから比較するようなものです。
（変数名の頭の「Z.」は標準化した変数につくプレフィックスです）
そもそも、偏差値は、標準化した後に平均値が50、標準偏差が10になるように調整したものです。
準備ができましたので、今度は、標準化済みのZで始まる変数に対して重回帰分析をします。以下に結果を示します。
あ
> RegModel.1 <-
+ lm(Z.住宅価格~Z.B+Z.NOX+Z.教師比率+Z.古い町並み+Z.固定資産税+Z.雇用センターまでの距離+Z.高速道路+Z.豪邸率+Z.商業+Z.川の周辺+Z.低所得者の割合+Z.犯罪率+Z.部屋数,
+ data=Dataset)
> summary(RegModel.1)
Call:
lm(formula = Z.住宅価格 ~ Z.B + Z.NOX + Z.教師比率 + Z.古い町並み +
Z.固定資産税 + Z.雇用センターまでの距離 + Z.高速道路 + Z.豪邸率 +
Z.商業 + Z.川の周辺 + Z.低所得者の割合 + Z.犯罪率 + Z.部屋数,
data = Dataset)
Residuals:
Min 1Q Median 3Q Max
-1.69559 -0.29680 -0.05633 0.19322 2.84864
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.628e-16 2.294e-02 0.000 1.000000
Z.B 9.243e-02 2.666e-02 3.467 0.000573 ***
Z.NOX -2.238e-01 4.813e-02 -4.651 4.25e-06 ***
Z.教師比率 -2.243e-01 3.080e-02 -7.283 1.31e-12 ***
Z.古い町並み 2.119e-03 4.043e-02 0.052 0.958229
Z.固定資産税 -2.260e-01 6.891e-02 -3.280 0.001112 **
Z.雇用センターまでの距離 -3.378e-01 4.567e-02 -7.398 6.01e-13 ***
Z.高速道路 2.897e-01 6.281e-02 4.613 5.07e-06 ***
Z.豪邸率 1.177e-01 3.481e-02 3.382 0.000778 ***
Z.商業 1.534e-02 4.587e-02 0.334 0.738288
Z.川の周辺 7.420e-02 2.379e-02 3.118 0.001925 **
Z.低所得者の割合 -4.074e-01 3.938e-02 -10.347 < 2e-16 ***
Z.犯罪率 -1.010e-01 3.074e-02 -3.287 0.001087 **
Z.部屋数 2.911e-01 3.193e-02 9.116 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.516 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回帰係数（標準化済みの変数を用いているので、“標準化偏回帰係数”と呼ぶのが正しい）が10のマイナス2乗（e-02）が付いている極小のものを外し、絶対値が大きい順に並べます。
Estimate Std. Error t value Pr(>|t|)
Z.低所得者の割合 -4.074e-01 3.938e-02 -10.347 < 2e-16 ***
Z.雇用センターまでの距離 -3.378e-01 4.567e-02 -7.398 6.01e-13 ***
Z.部屋数 2.911e-01 3.193e-02 9.116 < 2e-16 ***
Z.高速道路 2.897e-01 6.281e-02 4.613 5.07e-06 ***
Z.固定資産税 -2.260e-01 6.891e-02 -3.280 0.001112 **
Z.教師比率 -2.243e-01 3.080e-02 -7.283 1.31e-12 ***
Z.NOX -2.238e-01 4.813e-02 -4.651 4.25e-06 ***
Z.豪邸率 1.177e-01 3.481e-02 3.382 0.000778 ***
Z.犯罪率 -1.010e-01 3.074e-02 -3.287 0.001087 ** 「住宅価格」との因果関係が高い順に5つ抜き出すと
低所得者の割合
雇用センターまでの距離
部屋数
高速道路
固定資産税
となりました。
トップの「低所得者の割合」は、低所得者の割合が上がると高い家は売れない（だから回帰係数はマイナスとなる）という、「所得に応じた価格の家を買っている」という当たり前の話でしょう。
2番目の「雇用センターまでの距離」は、ちょっとわかりにくいですが、雇用が不安定な人が多く住む貧しさを計測しようとしたのでしょうが、意図と反して回帰係数はマイナスとなりました。雇用センターが町の中心にあるからでしょうか？（例えば役所や公園などの公共施設までの距離との疑似相関？？）
3番目の「部屋数」は高額な家は部屋が多いとも考えられますので、どちらが原因で結果かわからないものです。
4番目の「高速道路」はボストンでの生活の利便性、5番目の「固定資産税」も、「部屋数」と同じく持ち家の価値に応じた納税なのでしょうか、、、。
これだけではわからない説明変数が多いですね。😅
逆に、関係（住宅価格への寄与）が薄い方を見てみましょう。
Z.古い町並み 2.119e-03 4.043e-02 0.052 0.958229
Z.商業 1.534e-02 4.587e-02 0.334 0.738288
Z.川の周辺 7.420e-02 2.379e-02 3.118 0.001925 **
Z.B 9.243e-02 2.666e-02 3.467 0.000573 ***「古い町並み」と「商業」はp値も大きい（*がついていない）ので住宅価格にはあまり寄与しないと考えた方がよさそうです。「川の周辺」かどうかも、災害の観点で見る人と、眺望の観点で見る人がいるのかもしれません。「B」はp値が小さいので寄与率は大きくないとしても住宅価格に影響を与えているとみた方がよさそうです。
≡ おわりに
今回は、回帰分析の3回目として「重回帰分析 後編」の話をしました。統計的に説明変数が結果に影響を与えている大きさの違いが判ることは実務に役立ちます。色々な応用が出来そうです。
次回は、「回帰分析のヒント」と題して、回帰分析の広がり的な話をしたいと思っています。対数軸とか、ゼロの扱いとか、カテゴリ変数の扱いとか、ポアソン分布の場合はどうするとか、、、。
私も苦労しているところなのでこんな課題があるよねと課題の列挙で終わってしまうかもしれません。
次回で「回帰分析」のテーマは終わる予定です。残りは「クロス集計とカイ二乗検定」→「管理図とEVM」→「確率分布」となります。
