はじめに

　『Pythonで学ぶ空間データサイエンス入門』の独学時のまとめノートです。「定義編・導出編」「作図編・実装編」「可視化編」の三部構成でモデルやアルゴリズムの理解を目指します。
　本の内容から寄り道・回り道しながら進めます。本を読んだ上で補助的に読んでください。

　この記事では、GMIの定義について、数式を使って解説します。

【前の内容】

www.anarchive-beta.com

【他の内容】

www.anarchive-beta.com

【今回の内容】

はじめに
2.3 グローバル・モランのIの定義式の導出：モラン散布図のOLS
参考文献
おわりに

2.3 グローバル・モランのIの定義式の導出：モラン散布図のOLS

　空間的自己相関(spatial autocorrelation)の指標であるグローバル・モランのI(GMI・global Moran's I)の定義を数式で確認します。この記事では、モラン散布図(Moran scatterplot)の関係を用いて、線形回帰モデル(linear regression model)に対する最小二乗法(OLS・ordinary least square)を求めることで、GMIの定義式を導出します。
　空間ラグについては「2.3：空間ラグの定義式【空間DS入門のノート】 - からっぽのしょこ」、GMIについては「2.3：グローバル・モランのIの定義式【空間DS入門のノート】 - からっぽのしょこ」、モラン散布図については「【Python】2.3：モラン散布図の作成【空間DS入門のノート】 - からっぽのしょこ」、線形回帰モデルのOLSについては「3.2.0：線形回帰モデルの定義式【空間DS入門のノート】 - からっぽのしょこ」と「3.2.1-3：線形回帰モデルの最小二乗法の導出【空間データサイエンス入門のノート】 - からっぽのしょこ」を参照してください。

モデルの設定

　まずは、モラン散布図における線形回帰モデルの定義(仮定)を数式で確認します。
　空間重み行列については「2.2：空間重み行列の定義式：境界の共有【空間DS入門のノート】 - からっぽのしょこ」、多変量正規分布については「多次元ガウス分布の定義式 - からっぽのしょこ」を参照してください。

　横軸が偏差 $x_n - \bar{x}$ 、縦軸が空間ラグ $\mathrm{SL}_n$ であるモラン散布図について、次のような回帰式を考えます。

$\displaystyle \mathbf{W} \mathbf{d} = \beta \mathbf{d} + \boldsymbol{\epsilon} \tag{1}$

　各種の変数やパラメータなどについて「変数・パラメータ」などで確認します。

定義式：空間統計量

　モデルの確認の前に、空間ラグとGMIの定義式を確認します。

　地域 $n$ の空間ラグ $\mathrm{SL}_n$ は、次の式で定義されます。 $N$ 個の地域をまとめて $N$ 次元ベクトル $\mathbf{SL}$ とします。

$\begin{align} \mathrm{SL}_i &= \mathbf{w}_i^{\top} (\mathbf{x} - \bar{x} \mathbf{1}) \\ &= \sum_{j=1}^N w_{ij} (x_j - \bar{x}) \tag{2}\\ \mathbf{SL} &= \mathbf{W} (\mathbf{x} - \bar{x} \mathbf{1}) \end{align}$

　地域全体のGMIは、次の式で定義されます。

$\begin{align} \mathrm{GMI} &= \frac{1}{S_0} \frac{1}{s^2} (\mathbf{x} - \bar{x} \mathbf{1})^{\top} \mathbf{W} (\mathbf{x} - \bar{x} \mathbf{1}) \tag{3}\\ &= \frac{1}{S_0} \frac{1}{s^2} \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x}) (x_j - \bar{x}) \\ &= \frac{1}{\sum_{i=1}^N \sum_{j=1}^N w_{ij}} \frac{ \sum_{i=1}^N \sum_{j=1}^N w_{ij} (x_i - \bar{x}) (x_j - \bar{x}) }{ \frac{1}{N} \sum_{i=1}^N (x_i - \bar{x})^2 } \end{align}$

　ここで、 $\bar{x}$ は変数 $\mathbf{x}$ の標本平均、 $s^2$ は変数 $\mathbf{x}$ の標本分散です。また、空間重み行列 $\mathbf{W}$ の総和を $S_0$ とおきます。

　GMIの定義式を求めていきます。

変数・パラメータ

　線形回帰モデルの変数やパラメータなどの表記を確認します。

　地域数(データ数)を $N$ 、地域番号(データ番号・インデックス)を $n \in \{1, 2, \dots, N\}$ とします。 $n$ 番目の地域(データ)や地域 $n$ と呼びます。場合によっては、 $i, j$ などでも地域を表します。

　地域 $n$ のデータ値を $x_n$ 、 $N$ 個の地域をまとめて $N$ 次元ベクトル $\mathbf{x}$ とします。
　地域 $n$ のデータ値の偏差を $d_n = x_n - \bar{x}$ 、 $N$ 個の地域をまとめて $N$ 次元ベクトル $\mathbf{d}$ とします。

$\begin{align} \mathbf{d} &= \mathbf{x} - \bar{x} \mathbf{1} \tag{4}\\ &= \begin{pmatrix} x_1 - \bar{x} \\ x_2 - \bar{x} \\ \vdots \\ x_N - \bar{x} \end{pmatrix} \\ &= \begin{pmatrix} d_1 \\ d_2 \\ \vdots \\ d_N \end{pmatrix} \end{align}$

　( $\mathbf{x}$ ではなく) $\mathbf{d}$ を説明変数とします。

　地域 $n$ に関する重みを $\mathbf{w}_n$ 、 $N$ 個の地域に関する重みをまとめて $N \times N$ の空間重み行列 $\mathbf{W}$ とします。

$\displaystyle \mathbf{W} = \begin{pmatrix} \mathbf{w}_1^{\top} \\ \mathbf{w}_2^{\top} \\ \vdots \\ \mathbf{w}_N^{\top} \end{pmatrix} = \begin{pmatrix} w_{11} & w_{12} & \cdots & w_{1N} \\ w_{21} & w_{22} & \cdots & w_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ w_{N1} & w_{N2} & \cdots & w_{NN} \end{pmatrix}$

　 $\mathbf{W}$ の各行( $\mathbf{w}_n$ の各要素)は、0から1の値で総和が1になる値の条件を満たす必要があります。

$\displaystyle 0 \leq w_{ij} \leq 1 ,\ \sum_{j=1}^N w_{ij} = 1$

　つまりこのモデルでは、 $\mathbf{W}$ が行ごとに正規化されている必要があります。

　回帰パラメータを $\beta$ とします。
　地域 $n$ の被説明変数を $y_n$ 、 $N$ 個の地域をまとめて $N$ 次元ベクトル $\mathbf{y}$ とします。
　地域 $n$ の空間ラグ $\mathrm{SL}_n$ を被説明変数 $y_n$ とします。

$\begin{align} \mathbf{y} &= \mathbf{W} \mathbf{d} \tag{5}\\ y_n &= \mathbf{w}_n^{\top} \mathbf{d} \\ &= \sum_{j=1}^N w_{nj} d_j \\ &= \sum_{j=1}^N w_{nj} (x_j - \bar{x}) \\ &= \mathrm{SL}_n \tag{2'} \end{align}$

　地域 $n$ の誤差項を $\epsilon_n$ 、 $N$ 個の地域をまとめて $N$ 次元ベクトル $\boldsymbol{\epsilon}$ とします。
　誤差項 $\boldsymbol{\epsilon}$ は、平均ベクトル $\boldsymbol{\mu} = \mathbf{0}$ ・分散共分散行列 $\boldsymbol{\Sigma} = \sigma^2 \mathbf{I}$ の多変量正規分布(多次元ガウス分布)に従って生成されると仮定します。

$\displaystyle \begin{aligned} \boldsymbol{\epsilon} &\sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}) \\ \epsilon_n &\sim \mathcal{N}(0, \sigma^2) \end{aligned}$

　このとき、各誤差項 $\epsilon_n$ は、平均 $\mu = 0$ ・分散 $\sigma^2$ (標準偏差 $\sigma$ )の一変量正規分布(1次元ガウス分布)に従って、独立に生成されると言えます。

　以上の変数やパラメータなどを用いてモデルを定義します。

定義式：回帰モデル

　続いて、線形回帰モデルの定義式を確認します。

　モラン散布図における線形回帰モデルは、次の式で定義されます。この式は、各データの被説明変数 $y_n = \mathrm{SL}_n$ の計算式を表します。

$\begin{align} \mathbf{y} &= \beta \mathbf{d} + \boldsymbol{\epsilon} \tag{1'}\\ y_n &= \beta d_n + \epsilon_n \\ &= \beta (x_n - \bar{x}) + \epsilon_n \end{align}$

　以上で、モデルを定義しました。

推定値

　線形回帰モデルのパラメータ推定における推定値などの表記を確認します。

　回帰パラメータの推定値を $\hat{\beta}$ とします。
　地域 $n$ の被説明変数の推定値を $\hat{y}_n$ 、 $N$ 個の地域をまとめて $N$ 次元ベクトル $\hat{\mathbf{y}}$ とします。
　被説明変数の推定値 $\hat{\mathbf{y}}$ は、回帰パラメータの推定値(推定したパラメータ) $\hat{\beta}$ を用いて求まる予測値です。

$\begin{align} \hat{\mathbf{y}} &= \hat{\beta} \mathbf{d} \tag{6}\\ \hat{y}_n &= \hat{\beta} d_n \\ &= \hat{\beta} (x_n - \bar{x}) \end{align}$

　観測値(実現値) $\mathbf{y}$ と予測値 $\hat{\mathbf{y}}$ の差を残差 $\mathbf{e}$ とします。

$\begin{align} \mathbf{e} &= \mathbf{y} - \hat{\mathbf{y}} \tag{7}\\ e_n &= y_n - \hat{y}_n \end{align}$

　以上の推定値なども用いてパラメータを推定します。

計算式

　さらに、線形回帰モデルに関する計算式を求めます。
　式変形については「GMIの定義式」も参照してください。

　標本分散 $s^2$ は、次のように計算できます。

$\begin{align} s^2 &= \frac{1}{N} \mathbf{d}^{\top} \mathbf{d} \\ &= \frac{1}{N} \sum_{n=1}^N d_n^2 \\ &= \frac{1}{N} \sum_{n=1}^N (x_n - \bar{x})^2 \tag{8} \end{align}$

　空間重み行列 $\mathbf{W}$ の総和を $S_0$ とおき、次のように計算できます。

$\begin{align} S_0 &= \mathbf{1}^{\top} \mathbf{W} \mathbf{1} \\ &= \sum_{i=1}^N \sum_{j=1}^N w_{ij} \tag{9}\\ &= \sum_{i=1}^N 1 \\ &= N \end{align}$

　ただし、このモデルで用いられる空間重み行列 $\mathbf{W}$ は、行(地域)ごとに正規化されている(各要素が0から1の値 $0 \leq w_{ij} \leq 1$ で行和が1になる値 $\sum_{j=1}^N w_{ij} = 1$ の条件を満たす)とします。

　残差 $\mathbf{e}$ の式(7)の予測値 $\hat{\mathbf{y}}$ に式(6)を代入します。

$\begin{align} \mathbf{e} &= \mathbf{y} - \hat{\mathbf{y}} \tag{7}\\ &= \mathbf{W} \mathbf{d} - \hat{\beta} \mathbf{d} \tag{7'} \end{align}$

　観測値 $\mathbf{d}$ を用いた式に変形しました。

　残差 $\mathbf{e}$ の内積(残差平方和)を計算します。

$\begin{align} \mathbf{e}^{\top} \mathbf{e} &= \Bigl( \mathbf{W} \mathbf{d} - \hat{\beta} \mathbf{d} \Bigr)^{\top} \Bigl( \mathbf{W} \mathbf{d} - \hat{\beta} \mathbf{d} \Bigr) \\ &= (\mathbf{W} \mathbf{d})^{\top} \mathbf{W} \mathbf{d} - (\mathbf{W} \mathbf{d})^{\top} \hat{\beta} \mathbf{d} - (\hat{\beta} \mathbf{d})^{\top} \mathbf{W} \mathbf{d} + (\hat{\beta} \mathbf{d})^{\top} \hat{\beta} \mathbf{d} \\ &= \mathbf{d}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{d} - \hat{\beta} \mathbf{d}^{\top} \mathbf{W}^{\top} \mathbf{d} - \hat{\beta} \mathbf{d}^{\top} \mathbf{W} \mathbf{d} + \hat{\beta}^2 \mathbf{d}^{\top} \mathbf{d} \\ &= \mathbf{d}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{d} - 2 \hat{\beta} \mathbf{d}^{\top} \mathbf{W} \mathbf{d} + \hat{\beta}^2 \mathbf{d}^{\top} \mathbf{d} \tag{10} \end{align}$

途中式の途中式(クリックで展開)

1: 式(7')より、 $\mathbf{e}$ の内積の式を立てます。
2: 転置の性質 $(\mathbf{A} + \mathbf{B})^{\top} = \mathbf{A}^{\top} + \mathbf{B}^{\top}$ より、括弧を展開します。
3: 転置の性質 $(\mathbf{A} \mathbf{B})^{\top} = \mathbf{B}^{\top} \mathbf{A}^{\top}$ より、括弧を展開します。
4: 中2つの項は二次形式(スカラ)なので、転置できます。

　2つの項が一致するのでまとめます。

$\displaystyle \begin{aligned} \mathbf{d}^{\top} \mathbf{W}^{\top} \mathbf{d} &= \Bigl( \mathbf{d}^{\top} \mathbf{W}^{\top} \mathbf{d} \Bigr)^{\top} \\ &= \mathbf{d}^{\top} (\mathbf{W}^{\top})^{\top} \mathbf{d} \\ &= \mathbf{d}^{\top} \mathbf{W} \mathbf{d} \end{aligned}$

転置の性質 $(\mathbf{A}^{\top})^{\top} = \mathbf{A}$ より、2重の転置は元に戻ります。

　以上で、線形回帰モデルの定義を確認しました。

スポンサードリンク

パラメータの推定

　次は、線形回帰モデルに対する最小二乗法によるパラメータの計算式(最小二乗推定量・最小二乗解)を導出します。

回帰パラメータの最小二乗推定量

　回帰パラメータの最小二乗推定量(残差平方和が最小となる推定値)を求めます。

　残差平方和 $\mathbf{e}^{\top} \mathbf{e}$ を回帰パラメータの推定値 $\hat{\beta}$ に関して微分します。

$\begin{align} \frac{\partial \mathbf{e}^{\top} \mathbf{e}}{\partial \hat{\beta}} &= \frac{\partial}{\partial \hat{\beta}} \Bigl\{ \mathbf{d}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{d} - 2 \hat{\beta} \mathbf{d}^{\top} \mathbf{W} \mathbf{d} + \hat{\beta}^2 \mathbf{d}^{\top} \mathbf{d} \Bigr\} \\ &= \frac{\partial}{\partial \hat{\beta}} \Bigl\{ \mathbf{d}^{\top} \mathbf{W}^{\top} \mathbf{W} \mathbf{d} \Bigr\} - 2 \frac{\partial \hat{\beta}}{\partial \hat{\beta}} \mathbf{d}^{\top} \mathbf{W} \mathbf{d} + \frac{\partial \hat{\beta}^2}{\partial \hat{\beta}} \mathbf{d}^{\top} \mathbf{d} \\ &= 0 - 2 \hat{\beta}^{1-1} \mathbf{d}^{\top} \mathbf{W} \mathbf{d} + 2 \hat{\beta}^{2-1} \mathbf{d}^{\top} \mathbf{d} \\ &= - 2 \mathbf{d}^{\top} \mathbf{W} \mathbf{d} + 2 \hat{\beta} \mathbf{d}^{\top} \mathbf{d} \tag{11}\\ &= - 2 \mathbf{d}^{\top} \Bigl( \mathbf{W} \mathbf{d} - \hat{\beta} \mathbf{d} \Bigr) \\ &= - 2 \mathbf{d}^{\top} \mathbf{e} \end{align}$

途中式の途中式(クリックで展開)

1: 式(10)より、 $\mathbf{e}$ の内積の偏微分の式を立てます。 $\hat{\beta}$ に関する微分なので、 $\hat{\beta}$ 以外の項は定数として扱います。
2: 和の微分 $\frac{d \{f(x) + g(x)\}}{d x} = \frac{d f(x)}{d x} + \frac{d g(x)}{d x}$ より、微分の和に分割します。
2: 定数倍の微分 $\frac{d \{a f(x)\}}{d x} = a \frac{f(x)}{d x}$ より、係数を微分の外に出します。
3: べき乗の微分 $\frac{d x^n}{d x} = n x^{n-1}$ より、微分を計算します。
4: 指数の性質 $x^0 = 1$ より、0乗の項が1になり消えます。
5: 共通の項と $-1$ を括り出します。
6: 式(7')より、丸括弧を残差に置き換えます。

　残差平方和の微分 $\frac{\partial \mathbf{e}^{\top} \mathbf{e}}{\partial \hat{\beta}}$ が0となる回帰パラメータの推定値 $\hat{\beta}$ を求めます。

$\begin{align} && - 2 \mathbf{d}^{\top} \mathbf{W} \mathbf{d} + 2 \hat{\beta} \mathbf{d}^{\top} \mathbf{d} &= 0 \\ \Rightarrow && 2 \hat{\beta} \mathbf{d}^{\top} \mathbf{d} &= 2 \mathbf{d}^{\top} \mathbf{W} \mathbf{d} \\ \Rightarrow && \hat{\beta} &= \frac{ \mathbf{d}^{\top} \mathbf{W} \mathbf{d} }{ \mathbf{d}^{\top} \mathbf{d} } \equiv \hat{\beta}_{\mathrm{OLS}} \tag{12} \end{align}$

途中式の途中式(クリックで展開)

1: 式(11)より、 $\frac{\partial \mathbf{e}^{\top} \mathbf{e}}{\partial \hat{\beta}} = 0$ とおいて、 $\hat{\beta}$ について解きます。
2: $\hat{\beta}$ と無関係な項を右辺に移します。
3: 両辺を $2 \mathbf{d}^{\top} \mathbf{d}$ で割ります。

　回帰パラメータ $\beta$ の最小二乗推定量 $\hat{\beta}_{\mathrm{OLS}}$ の計算式が得られました。

　以上で、モラン散布図における線形回帰モデルのパラメータの推定値の計算式が求まりました。

スポンサードリンク

グローバル・モランのIとモラン散布図の関係

　最後は、モラン散布図における回帰パラメータの最小二乗推定量とグローバル・モランのIの関係を導出します。

　回帰パラメータ $\beta$ の最小二乗推定量の式(12)を変形します。

$\begin{align} \hat{\beta}_{\mathrm{OLS}} &= \frac{ \mathbf{d}^{\top} \mathbf{W} \mathbf{d} }{ \mathbf{d}^{\top} \mathbf{d} } \tag{12}\\ &= \frac{N}{N} \frac{1}{\mathbf{d}^{\top} \mathbf{d}} \mathbf{d}^{\top} \mathbf{W} \mathbf{d} \\ &= \frac{1}{N} \frac{1}{\frac{1}{N} \mathbf{d}^{\top} \mathbf{d}} \mathbf{d}^{\top} \mathbf{W} \mathbf{d} \\ &= \frac{1}{S_0} \frac{1}{s^2} (\mathbf{x} - \bar{x} \mathbf{1})^{\top} \mathbf{W} (\mathbf{x} - \bar{x} \mathbf{1}) \\ &= \mathrm{GMI} \tag{3'} \end{align}$

途中式の途中式(クリックで展開)

1: 式(12)を再掲しています。
2: $\frac{N}{N} = 1$ を掛けます。
3: 逆数の逆数 $N = \frac{1}{\frac{1}{N}}$ より、分子の $N$ を分母に移します。
4: 式(9)より、地域数 $N$ を $\mathbf{W}$ の総和とみなして置き換えます。
4: 式(8)より、2つ目の項を標本分散に置き換えます。
5: 式(3)より、式全体をGMIに置き換えます。

　回帰パラメータ $\beta$ の最小二乗推定量とGMIが一致しました。
　また、予測値の式(6)より、モラン散布図に対する回帰直線は $\hat{y}_n = \hat{\beta} (x_n - \bar{x})$ なので、GMIが回帰直線の傾きであることが分かります。

　以上で、モラン散布図における線形回帰モデルの最小二乗推定量とグローバル・モランのIの関係式が得られました。

　この記事では、グローバル・モランのIの定義式を数式で確認しました。

参考文献

Pythonで学ぶ空間データサイエンス入門: 地域の特徴を発見する方法

作者:植杉　大
Independently published

Amazon

おわりに

　2冊目の本のGMIの章に「モラン散布図の回帰直線の傾きはGMIに対応している」という一文があり、気になったのでやってみました。初稿を書いている(1周目を読んでいる)ときには、相関係数の定義式に重み行列を組み込む形でGMIの定義式を導出したいと思っていたのですが、結局できませんでした。そもそもこれはできることなんですかね。
　今回というかここ暫く実装や作図の記事が続いていたこともあり、少しは歯応えのある数式弄りをと思ってやってみたのですが、偏差をd、空間ラグをyと置いたら(既に修正を終えている3章の)素の線形回帰モデルよりも簡単で拍子抜けでした。
　逆に言うと、本には載っていない内容ですが、3章の前に読んでおくといい準備運動になると思います。今回のは1次元版の線形回帰なので内積の計算が主ですが、3章のは多次元版の線形回帰なので行列の計算ががっつり出てきます。
　ただ若干筆が乗らなくて微妙に解説をサボりました。空間ラグやGMIの定義式と線形回帰モデルの定義や導出の記事も合わせて読んでください。最後に書いてもアレですが。

　最後に、先日公開された鈴木愛理の新曲のライブ映像をどうぞ♪

　楽曲についての感想ではないですが、ショートカット可愛えぇ。圧巻のパフォーマンスについて触れるべきかもしれませんが、ショートカット可愛ぇぇ。

【次の内容】

　GMIの計算をプログラムで確認します。

www.anarchive-beta.com

　モラン散布図とGMIの関係を図で確認します。

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

2.3：グローバル・モランのIの定義式の導出：モラン散布図のOLS【空間DS入門のノート】

はじめに

2.3 グローバル・モランのIの定義式の導出：モラン散布図のOLS

モデルの設定

定義式：空間統計量

変数・パラメータ

定義式：回帰モデル

推定値

計算式

パラメータの推定

回帰パラメータの最小二乗推定量

グローバル・モランのIとモラン散布図の関係

参考文献

おわりに