回帰分析 – TauStation

概要

単回帰が $(X, Y)$ への線形関係を扱うのに対して、重回帰は複数の説明変数 $X_k \; (k=1 \ldots m)$ と $Y$ の線形関係を扱う。

データ $(X_1, \ldots , X_m, Y) = (x_{1i}, \ldots, x_{mi}, y_i), (i = 1 \ldots n)$ に対して以下の線形式で最も説明性の高いものを求める。

(1) $\begin{equation*} \hat{y} = w_0 + w_1 x_1 + \cdots + w_m x_m \end{equation*}$

そのために、データとその推測値の残差 $\hat{y_i} - y_i$ の平方和が最小となるような係数 $w_0, w_1, \ldots, w_m$ を最小二乗法により求める。

定式化

説明変数とターゲットのデータセット $n$ 組が次のように得られているとする。

(2) $\begin{equation*} x_{1i}, \ldots , x_{mi}, y_i \quad {\rm where} \quad i = 1, \ldots , n \end{equation*}$

残差の平方和については

(3) $\begin{equation*} \min \sum_{i=1}^n (\hat{y} -y_i)^2 \quad {\rm where} \quad \hat{y}_i = w_0 + w_1 x_{1i} + \cdots w_m x_{mi} \end{equation*}$

ここで残差を最小化する $w_0, w_1, \ldots, w_m$ を求めるために、それぞれで偏微分する。

(4) $\begin{equation*} \begin{array}{c} \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial w_0} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \displaystyle \frac{\partial}{\partial w_1} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \vdots \\ \displaystyle \frac{\partial}{\partial w_m} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \end{array} \right. \end{equation*}$

変形すると、

(5) $\begin{equation*} \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) = 0 \\ \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) x_{1i} = 0 \\ \vdots \\ \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) x_{mi} = 0 \\\end{array} \right. \end{array} \end{equation*}$

更にこれを展開するのに、 $S_k = \sum_{i=1}^n x_{ki}$ , $S_{kl} = \sum_{i=1}^n x_{ki} x_{li}$ , $S_y = \sum_{i=1}^n y_i$ , $S_{ky} = \sum_{i=1}^n x_{ki} y_i$ とおいて

(6) $\begin{equation*} \left\{ \begin{array}{l} n w_0 + w_1 S_1+ \cdots + w_m S_m - S_y = 0 \\ w_0S_1 + w_1 S_{11} + \cdots + w_m S_{1m} - S_{1y} = 0 \\ \vdots \\ w_0 S_m + w_1 S_{m1} + \cdots + w_m S_{mm} - S_{my} = 0 \\ \end{array} \right. \end{array} \end{equation*}$

これを行列形式で表示すると

(7) $\begin{equation*} \left( \begin{array}{cccc} n & S_1 & \cdots & S_m \\ S_1 & S_{11} & \cdots & S_{1m} \\ \vdots & \vdots & & \vdots \\ S_m & S_{m1} & \cdots & S_{mm} \end{array} \right) \left( \begin{array}{c} w_0 \\ w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} S_y \\ S_{1y} \\ \vdots \\ S_{my} \end{array} \right) \end{equation*}$

解の導出

式(7)の連立方程式を解けば各係数を得ることができる。ここで左辺の係数行列は以下のようにも表せる。

(8) $\begin{equation*} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \end{array} \right) \left( \begin{array}{cccc} 1 & x_{11} & \cdots & x_{m1} \\ 1 & x_{12} & \cdots & x_{m2} \\ \vdots & \vdots & & \vdots \\ 1 & x_{1n} & \cdots & x_{mn} \end{array} \right) = \boldsymbol{X}^T \boldsymbol{X} \end{equation*}$

また右辺は次のように表せる。

(9) $\begin{equation*} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \end{array} \right) \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) = \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

したがって式(7)は以下のように表せる。

(10) $\begin{equation*} \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} = \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

これを解くと

(11) $\begin{equation*} \boldsymbol{w} = ( \boldsymbol{X}^T \boldsymbol{X} )^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

ベクトル・行列による表現と解

定式化の段階からベクトル・行列を用いてみる。式(1)の各変数を全てベクトル形式で表すと以下のようになる。

(12) $\begin{equation*} \hat{\boldsymbol{y}} = \boldsymbol{x}^T \boldsymbol{w} \end{equation*}$

ここで、

(13) $\begin{equation*} \boldsymbol{x} = \left( \begin{array}{c} 1 & x_1 & \vdots & x_n \end{array} \right) \quad , \quad \boldsymbol{w} = \left( \begin{array}{c} w_0 & w_1 & \vdots & w_n \end{array} \right) \end{equation*}$

複数データからなるデータセットに対して、式(8)の行列表現を用いると、以下のように表現できる。

(14) $\begin{equation*} \hat{\boldsymbol{y}} = \boldsymbol{X} \boldsymbol{w} \end{equation*}$

各データの残差については以下のようになる。

(15) $\begin{equation*} \hat{\boldsymbol{y}} - \boldsymbol{y} = \boldsymbol{X} \boldsymbol{w} - \boldsymbol{y} \end{equation*}$

ここで残差の二乗和を損失関数Lとして、これを行列で表現すると以下のようになる。

(16) $\begin{align*} L &= \left( \hat{\boldsymbol{y}} - \boldsymbol{y} \right)^T \left( \hat{\boldsymbol{y}} - \boldsymbol{y} \right) \\ &= \left( \boldsymbol{Xw} - \boldsymbol{y} \right)^T \left( \boldsymbol{Xw} - \boldsymbol{y} \right) \\ &= \left( \boldsymbol{w}^T \boldsymbol{X}^T - \boldsymbol{y}^T \right) \left( \boldsymbol{Xw} - \boldsymbol{y} \right) \\ &= \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} - \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} + \boldsymbol{y}^T \boldsymbol{y} \end{align*}$

上式の第2項、第3項については以下のように整理できる。

(17) $\begin{align*} - \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} = - \left( \boldsymbol{Xw} \right)^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} = -2 \boldsymbol{y}^T \boldsymbol{Xw} \end{align*}$

損失関数Lはスカラーであり、これを最小にするためにベクトルwで微分し、それらの値がゼロとなるような方程式とする。

(18) $\begin{align*} \frac{dL}{d\boldsymbol{w}} = \frac{d}{d\boldsymbol{w}} \left( \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} - 2 \boldsymbol{y}^T \boldsymbol{Xw} + \boldsymbol{y}^T \boldsymbol{y} \right) = 2 \boldsymbol{X}^T \boldsymbol{Xw} -2 \boldsymbol{X}^T \boldsymbol{y} = \boldsymbol{0} \end{align*}$

なお第1項については、式(7)と式(8)より $\boldsymbol{X}^T\boldsymbol{X}$ は対象行列なので転置しても同じ行列となり、行列の微分の公式から以下のようになることを利用している。

(19) $\begin{align*} \frac{d}{d\boldsymbol{x}} \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} = \left( \boldsymbol{X}^T \boldsymbol{X} + \left(\boldsymbol{X}^T \boldsymbol{X} \right)^T \right) \boldsymbol{w} = 2 \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} \end{align*}$

式(18)をwについて解いて以下を得る。

(20) $\begin{equation*} \boldsymbol{w} = \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

多重共線性

式(7)から $w_0$ を消去する。

(21) $\begin{equation*} \left( \begin{array}{ccc} S_{11} - \dfrac{{S_1}^2}{n} & \cdots & S_{1m} - \dfrac{S_1 S_m}{n}\\ \vdots & & \vdots \\ S_{m1} - \dfrac{S_m S_1}{n} & \cdots & S_{mm} - \dfrac{{S_m_}^2}{n} \end{array} \right) \left( \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} S_{1y} - \dfrac{S_1 S_y}{n} \\ \vdots \\ S_{my} - \dfrac{S_m S_y}{n} \end{array} \right) \end{equation*}$

これを分散・共分散で表すと以下の通り。ただし以下で $V_{ii}=\rm{Var}(X_i)$ 、 $V_{ij} = \rm{Cov}(X_i, X_j)$ と表している。

(22) $\begin{equation*} \left( \begin{array}{ccc} V_{11} & \cdots & V_{1m}\\ \vdots & & \vdots \\ V_{m1} & \cdots & V_{mm} \end{array} \right) \left( \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} V_{1y} \\ \vdots \\ V_{my} \end{array} \right) \end{equation*}$

ここで $x_j = a x_i + b$ と2つの間に完全な線形関係がある場合、分散・共分散の性質から以下の関係が成り立つ。

(23) $\begin{equation*} V_{jj} = a^2V_{ii}, \; V_{ji} = V_{ij} = aV_{ii}, \; V_{jk} = V_{kj} = aV_{ji} = aV_{ij} \end{equation*}$

これらを式(22)の係数行列に適用すると、

(24) $\begin{align*} &\left[ \begin{array}{ccccccc} V_{11} & \cdots & V_{1i} & \cdots & V_{1j} & \cdots & V_{1m}\\ \vdots && \vdots && \vdots && \vdots\\ V_{i1} & \cdots & V_{ii} & \cdots & V_{ij} & \cdots & V_{im}\\ \vdots && \vdots && \vdots && \vdots\\ V_{j1} & \cdots & V_{ji} & \cdots & V_{jj} & \cdots & V_{jm}\\ \vdots && \vdots && \vdots && \vdots\\ V_{m1} & \cdots & V_{mi} & \cdots & V_{mj} & \cdots & V_{mm} \end{array} \right]\\ &= \left[ \begin{array}{ccccccc} V_{11} & \cdots & V_{1i} & \cdots & aV_{1i} & \cdots & V_{1m}\\ \vdots && \vdots && \vdots && \vdots\\ V_{i1} & \cdots & V_{ii} & \cdots & aV_{ii} & \cdots & V_{im}\\ \vdots && \vdots && \vdots && \vdots\\ aV_{i1} & \cdots & aV_{ii} & \cdots & a^2V_{ii} & \cdots & aV_{im}\\ \vdots && \vdots && \vdots && \vdots\\ V_{m1} & \cdots & V_{mi} & \cdots & aV_{mi} & \cdots & V_{mm} \end{array} \right] \end{align*}$

係数行列の行(列)が一次従属となっていて、その行列式がゼロとなることから、連立方程式は解を持たない。変数間の関係が線形でなくても相関がかなり高い場合は行列式がゼロに近づき、解が不安定になる。

概要

Pythonを使った統計計算と図示の練習のため、コンビニエンスストアで単回帰分析をやってみた。

コンビニの店舗数は「商業動態統計年報」の2016年データを使い、説明変数として2015年の国勢調査人口、2018年の国土地理院による都道府県面積、2017年道路統計年報の道路実延長データを使った。

計算コードは以下の通り。

import numpy as np
import matplotlib.pyplot as plt

master_data = np.loadtxt('conv-store-data.csv', delimiter=',',
        skiprows=2, usecols=(1, 2), encoding='utf8')
print(master_data)

y = num_stores = master_data[:,0]
x = population = master_data[:,1]

sxy = np.sum((x - x.mean())*(y - y.mean()))
sxx = np.sum((x - x.mean())**2)
syy = np.sum((y - y.mean())**2)
a = sxy / sxx
b = y.mean() - a*x.mean()
r = sxy / np.sqrt(sxx * syy)
n = len(num_stores)
ESS = n * sxy**2 / sxx
TSS = n * syy
R2 = ESS / TSS

max_population = 15000000
max_stores = 8000
plt.xlim(0, max_population)
plt.ylim(0, max_stores)
plt.xlabel('population')
plt.ylabel('number of stores')

plt.scatter(population, num_stores, label='#stores')

x = np.linspace(0, max_population)
plt.plot(x, a*x + b, color='m')

plt.text(1000000, 7600, 'a = ' + str(a))
plt.text(1000000, 7300, 'b = ' + str(b))
plt.text(1000000, 7000, 'r = '+str(r))
plt.text(1000000, 6700, 'R2 = '+str(R2))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

master_data = np.loadtxt('conv-store-data.csv', delimiter=',',

skiprows=2, usecols=(1, 2), encoding='utf8')

print(master_data)

y = num_stores = master_data[:,0]

x = population = master_data[:,1]

sxy = np.sum((x - x.mean())*(y - y.mean()))

sxx = np.sum((x - x.mean())**2)

syy = np.sum((y - y.mean())**2)

a = sxy / sxx

b = y.mean() - a*x.mean()

r = sxy / np.sqrt(sxx * syy)

n = len(num_stores)

ESS = n * sxy**2 / sxx

TSS = n * syy

R2 = ESS / TSS

max_population = 15000000

max_stores = 8000

plt.xlim(0, max_population)

plt.ylim(0, max_stores)

plt.xlabel('population')

plt.ylabel('number of stores')

plt.scatter(population, num_stores, label='#stores')

x = np.linspace(0, max_population)

plt.plot(x, a*x + b, color='m')

plt.text(1000000, 7600, 'a = ' + str(a))

plt.text(1000000, 7300, 'b = ' + str(b))

plt.text(1000000, 7000, 'r = '+str(r))

plt.text(1000000, 6700, 'R2 = '+str(R2))

plt.show()

人口との関係

コンビニ店舗数と人口の散布図と回帰式を以下に示す。

相関係数が極めて高いのは当然で、やはりコンビニの出店計画には人口ファクターが強く作用していることがわかる。

定数項bが店舗数のオーダーに比べてほぼゼロというのも興味深い。

係数aの値からは5万人弱で1店舗ということになるが、人口10万くらいの都市で2店舗しかないことになり、ちょっと少ないような気がする。もしかすると、一定規模以上の市町村単位や都市単位くらいで層化して出店計画を立てているのかもしれない。

面積との関係

店舗数と面積の関係

次に店舗数と都道府県面積の関係を見てみた。

この結果はかなりはずれで、データを見ても人口が少なく面積が群を抜いて大きい北海道の影響を大きく受けている。

ここで、面積が極端に大きい北海道（83423.83㎢～1位、2906店舗～5位）と面積が小さいが集積度が極端に高い東京都（2193.96㎢～45位、7003店舗～1位）の2つを除いて計算してみる。

これはさらにおかしな結果で、面積が小さいほど店舗数が多いことになる。

考えてみれば、集客数を期待するなら人口が集積している地域が有利だから、人口密度に比例する可能性を考えた方がいいのかもしれない。もし面積が小さい県の方が集積度が高いと想定すると、面積だけを取り出したときに逆の関係になるとも考えられるが、相関係数や決定係数が小さすぎるので考察は難しい。

人口密度

以下は人口密度との関係。

今度はかなりきれいに相関の高さが出ている。

直接的な計算式に入れているかどうかわからないが、GISなどで出店計画を立てるとしたら、人口密度の高いエリアを選んでいくだろうことが想定される。

ただ、店舗数は人口などといった売り上げに直結するデータから導かれるのが普通で、人口密度が高くても人口が少なければ出店インセンティブにはならない。

人口と人口密度の関係

試しに人口を説明変数、人口密度を被説明変数として両者の関係を見てみると、驚くことに「人口が多いほど人口密度が高くなる（あるいはその逆）」という関係になる。

ここから先は人口論や地域論になりそうなので置いておくが、少なくとも日本においては、「狭いところほど人が集まっている傾向がある」ということになりそうである。

もちろんこれは他の国でも一般に当てはまることかもしれないが、朝のラッシュ時に特定の車両に無理やり乗り込んでいる割に離れた車両がすいているとか、1本電車を遅らせたらガラガラだったとか、そのあたりの行動パターンを見ていると、何となく日本に特有のような気がする。

道路延長との関係

最後に道路延長との関係を見てみる。

東京のように稠密な都市は例外とすると、概ね関係はありそうである。ただし相関係数、決定係数は高くない。

コンビニ店舗が道路の利便性に依っていることは推測できるが、やはり人口という売り上げ直結のデータに比べると関係は弱い。

高速道路の延長についても見てみたが、こちらはほとんど関係は見られなかった。

ただ、高速道路の延伸に伴ってコンビニエンスストアの店舗数が伸びているようであり、マクロな延長というよりも物流上のインパクトが大きいことは予想される。

TauStation

タグ: 回帰分析

回帰分析～重回帰

概要

定式化