はじめに

東京大学/株式会社Nospare リサーチャーの栗栖です．
この記事では前回の局所定数回帰についての記事から引き続き，ノンパラメトリックな回帰分析の統計手法の一つである局所線形回帰について紹介します．

モデル

ノンパラメトリックな回帰モデル

Y_i = m(\boldsymbol{X}_i) + U_i,\ E[U_i|\boldsymbol{X}_i] = 0

から独立同分布 (independent and identically distributed, i.i.d.) なデータ $(Y_i, \boldsymbol{X}_i)$, $i = 1,\dots,n$ が得られたとします．ただし，$Y_i \in \mathbb{R}$, $\boldsymbol{X}_i=(X_{1,i},\dots,X_{d,i})' \in \mathbb{R}^d$とします．

以下ではデータから回帰関数 $m(x) = E[Y_i|\boldsymbol{X}_i = x]$, $x = (x_1,\dots,x_d) \in \mathbb{R}^d$ を推定するという問題を考えましょう．推定にあたり，上記のモデルに対して以下を仮定します．

(A1:LL)
[1] $X_i$は密度関数 $f$ をもち，$f$ は $x$ のある近傍 $N_x$ で連続．
[2] $m$ は $N_x$ で2回連続偏微分可能．
[3] $\sigma^2(x) = E[U_i^2|\boldsymbol{X}_i = x]$ は $x$ で連続．

また $x \in \mathbb{R}^d$ で2回連続偏微分可能な関数 $g: \mathbb{R}^d \to \mathbb{R}$ に対して，

\partial_j g(x) = {\partial  \over \partial x_j}g(x),\ \partial_{jk}g(x) = {\partial^2  \over \partial x_j \partial x_k}g(x).

とします．

局所線形回帰

回帰関数 $m$ の推定にあたり，以下の最適化問題を考えましょう．

\begin{align*}
\hat{\beta} &= (\hat{\beta}_0,\dots, \hat{\beta}_d)\\ 
&= \underset{\beta \in \mathbb{R}^{d+1}}{\text{argmin}}\sum_{i=1}^n (Y_i - \beta_0 - \sum_{j=1}^{d}\beta_j(\boldsymbol{X}_i - x))^2 K\left({\boldsymbol{X}_i - x \over h}\right),\\
K\left({\boldsymbol{X}_i - x \over h}\right) &= K\left({X_{1,i} - x_1 \over h},\dots, {X_{d,i} - x_d \over h}\right)
\end{align*}

ここで，$h = h_{n} \to 0$, $n \to \infty$ はバンド幅，$K:\mathbb{R}^d \to \mathbb{R}$ はカーネル関数と呼ばれ，以下を満たすとします．

(A2:LL)
[1] $\int K(u)du=1$
[2] $K(u)$ は有界かつ $[-1,1]^d$ に台をもつ ($K(u) = 0$, $\|u\| >1$)
[3] $\check{u} = (1, u_1,\dots,u_d)'$として，$S=\int \check{u}\check{u}'K(u)du$ は正則(逆行列をもつ)

上記の最適化問題は $\boldsymbol{X}_i = x$ の近傍で「局所」的に $Y_i$ を「線形関数」に「回帰」していると解釈できるので，「局所線形回帰」と呼ばれます．また得られる解は「局所線形推定量」(local linear estimator, LL 推定量) と呼ばれることもあります．特にこの方法では，$\boldsymbol{X}_i = x$ の近傍で

Y_i \approx m(\boldsymbol{X}_i) \approx m(x) + \sum_{j=1}^{d}\partial_j m(x)(\boldsymbol{X}_i - x)

と近似できることを利用して $m(x)$ (とその偏微分 $\partial_j m(x)$) を推定しています．上記の最適化問題を解くと，$m(x)$ の LL推定量が以下で与えられることが分かります．

\begin{align*}
\hat{\beta} &= (\hat{m}(x),\partial_1 \hat{m}(x),\dots, \partial_d \hat{m}(x))\\ 
&=　\left(\sum_{i=1}^{n}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_i\tilde{\boldsymbol{X}}'_i \right)^{-1}\sum_{i=1}^{n}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_iY_i 
\end{align*}

ここで $\tilde{\boldsymbol{X}}_i = (1, (\boldsymbol{X}_i - x)')'$ です．

LL推定量の漸近正規性

最後にLL推定量の漸近的性質について紹介しておきます．まず(A1:LL), (A2:LL) に加えて以下を仮定します．

(A3:LL)
[1] $n \to \infty$ として，$h \to 0$, $nh^d \to \infty$, $nh^{d+4} \to c \in (0,\infty]$．
[2] また $x \in \mathbb{R}^d$ に対して $f(x)>0$，
[3] ある $\delta>0$ と $x$ の近傍 $N_x$ に対して $\sup_{z \in N_x}E[|U_i|^{2+\delta}|\boldsymbol{X}_i = z] \leq U(x)<\infty$．

仮定(A1:LL)-(A3:LL)の下で，$n \to \infty$ として以下が成り立ちます．

\begin{align*}
&\sqrt{nh^d}\left(
\left(
\begin{array}{c}
\hat{\beta}_0 - m(x)\\
h(\hat{\beta}_1 - \partial_1 m(x))\\
\vdots \\
h(\hat{\beta}_d - \partial_d m(x))
\end{array}
\right) - S^{-1}\left(
\begin{array}{c}
b_{0,n}(x)\\
b_{1,n}(x)\\
\vdots \\
b_{d,n}(x)
\end{array}
\right)
\right)\\ 
&\stackrel{d}{\to}N
\left( \boldsymbol{0}, {\sigma^2(x) \over f(x)}S^{-1} V S^{-1}
\right).
\end{align*}

ここで，$V = \int \check{u}\check{u}'K^2(u)du$,

\begin{align*}
b_{0,n}(x) &= h^2\sum_{j,k=1}^{d}{\partial_{jk}m(x) \over 2}\int u_j u_k K(u)du,\\
b_{r,n}(x) &= h^2\sum_{j,k=1}^{d}{\partial_{jk}m(x) \over 2}\int u_r u_j u_k K(u)du,\ r=1,\dots,d.
\end{align*}

LC推定量とLL推定量の比較

• 線形性の保存：$e_1 = (1,0,\dots,0)' \in \mathbb{R}^{d+1}$,

W(\boldsymbol{X}_i) = e'_1\left(\sum_{i=1}^{n}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_i\tilde{\boldsymbol{X}}'_i \right)^{-1}K\left({\boldsymbol{X}_i - x \over h}\right)\tilde{\boldsymbol{X}}_i

とおくと，

\hat{\beta}_0 = \sum_{i=1}^{n}W(\boldsymbol{X}_i)Y_i

と書けることから，$W(\boldsymbol{X}_i)$ は $m(x)$ を推定する際の $Y_i$ の $\boldsymbol{X}_i$ 貢献度と解釈することができます．LL推定量は線形関係を復元することが分かります．実際，任意の $a \in \mathbb{R}$, $b \in \mathbb{R}^d$ に対して $Y_i = a + b'\boldsymbol{X}_i$ とすると，

\begin{align*}
\sum_{i=1}^{n}W(\boldsymbol{X}_i)(a+b'\boldsymbol{X}_i) &= \sum_{i=1}^{n}W(\boldsymbol{X}_i)(a+b'x + b'(\boldsymbol{X}_i-x))\\
&= \sum_{i=1}^{n}W(\boldsymbol{X}_i)\tilde{\boldsymbol{X}}_i \left(
\begin{array}{c}
a + b'x \\
b
\end{array}
\right) \\
&= a + b'x
\end{align*}

が成り立ちます．これはLC推定量にはない性質です．

• 仮定に関する比較(カーネルの対称性の有無)

LL推定量では，LC推定量で仮定したカーネルの対称性 (A2:LC)[2]は必要ないことに注意しておきます．この性質は，$x$ が空間を2つに分ける境界線上にあり，データ $\boldsymbol{X}_i$ がその境界の片側でしか観測されない場合に $m(x)$ を推定したい場合なども扱える条件になっています．このような問題は因果推論における回帰非連続デザイン(regression discontinuity design, RDD) において現れ，LL推定量を用いて因果効果の推定が行われます．RDDについては石原先生の記事「回帰非連続デザインについて」で解説されています．

• バイアスに関する比較

$m(x)$ のLC推定量のバイアスが

b_{m,n}(x) = h^2\sum_{j,k=1}^d \left({\partial_{jk} m(x) \over 2} + {\partial_j m(x)\partial_k f(x) \over f(x)}\right)\int u_ju_k K(u)du

だったのに対し，$m(x)$ のLL推定量のバイアスは

b_{0,n}(x) = h^2\sum_{j,k=1}^{d}{\partial_{jk}m(x) \over 2}\int u_j u_k K(u)du

となり，${\partial_j m(x)\partial_k f(x) \over f(x)}$ の部分だけ LL推定量ではバイアスに出てくる項が少なくなっています．これは必ずしも LL推定量の方が LC推定量よりもバイアスが小さくなることを意味しませんが，実際に推定を行うと多くの場合に LL推定量の方がバイアスが小さくなることが知られています．

まとめ

この記事ではノンパラメトリック回帰分析の１つとして局所線形回帰を紹介しました．次の記事では局所線形回帰の拡張として局所多項式回帰について紹介する予定です．
株式会社Nospareには統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．