はじめに
こんにちは，株式会社Nospare の小林です．
今までNospareの記事で，以下のような分位点回帰を扱ったもの
- データの端に興味があるときの分位点回帰 (解説編)
- データの端に興味があるときの分位点回帰 (データ分析編)
- [R] quantregパッケージで分位点回帰をしてみる
- [R]分位点回帰をベイズ推定してみる
がいくつか公開されていますが，これらの記事では基本的には被説明変数が連続な場合をぼんやりと想定しています．
被説明変数が離散の場合には分位点回帰の取り扱いは，
- 異なる離散のタイプに対する一般的な理論が存在しない
- 推定が実用上難しい
- 同じ値が存在する場合の標本分位点の漸近的性質の扱いが難しい
という理由から，関連する研究はまだ発展途上です（Geraci and Farcomeni, 2021）．
本記事では最近のものも含め，カウントデータに対する分位点回帰のアプローチをいくつか紹介します．
離散データに対して何らかの連続性を導入
まずは離散データを連続なものに変換，あるいはパラメトリックな離散分布に対応する連続な確率分布を導入し，連続なデータに対する条件付き分位点のモデリングを行う方法を取り上げます．
ジッタリング（jittering）
Macahdo and Santos Silve (2005)によるこの方法は，離散の被説名変数$Y$に対し，連続なランダムノイズ$U$を乗せることを考えます：$Z=Y+U$．連続化された$Z$に対して，通常の分位点回帰の損失関数を最小化して分位点回帰の係数を推定します．これを繰り返し行うことでランダムな$U$について平均化します．ちなみに$Z_E=Y+E[U], E[U]=0.5$に対する分位点回帰を考えても漸近的に等しくなります．この方法のメリットはパラメトリックなモデリングに基づく効率性が挙げられますが，過適合や$Y$の値が小さいときに$U$の影響が大きくなってしまうなどの問題があります．
パラメトリック離散分布に対応する連続確率分布を導入
Padellini and Rue (2019)は，パラメトリック離散確率分布の連続バージョンを考え，連続バージョンの分位点に対するモデリングを行い，それを元の離散分布のパラメータに対してマッピングを行うということを考えました．まず連続バージョンの分布ですが，例えばポアソン分布に従う離散の$Y^*\sim POI(\lambda)$の場合，分布関数は
F_Y^*(y)=\frac{\Gamma(\lfloor y\rfloor+1,\lambda)}{\Gamma(\lfloor y \rfloor+1)}, \quad x\geq 0
と書けます．ここで不完全ガンマ関数$\Gamma(y,\lambda)=\int_\lambda^\infty e^{-s}s^{x-1}ds$．これに対する連続なポアソン分布$Y\sim ContinousPOI(\lambda)$の分布関数は
F_Y(y)=\frac{\Gamma(y+1,\lambda)}{\Gamma( y +1)}, \quad x>-1
で与えられます．
いま$Y_i^*|X_i\sim Poisson(\lambda_i)$のとき，連続ポアソンの分位点を$Q_\alpha(Y_i|X_i)=\exp(\eta_i^\alpha)$とモデル化すると，離散ポアソンのパラメータはこの連続ポアソンの分位点を使って
\lambda_i=\frac{\Gamma^{-1}(Q_\alpha(Y_i|X_i)^\alpha+1,1-\alpha)}{\Gamma(Q_\alpha(Y_i|X_i)^\alpha+1}
で表すことができます．連続バージョンの分位点に対するモデリングをすることで，推定された分位点は離散ではなくなりますが，$Q_\alpha(Y_i|X_i)=Q_\alpha(\lceil Y_i^* \rceil |X_i)=\lceil Q_\alpha(Y_i^* |X_i)\rceil$が成り立ちます．このアプローチの利点は，推定された異なるレベルの分位点が交差する（quantile crossing）ことが起きにくいということが挙げられます．
Lamarche et al. (2021)はこのアプローチをゼロ過剰カウント縦断データの場合に拡張しました．この論文でも離散確率分布の連続バージョンを扱っているのですが，連続バージョンの分位点をモデリングしてそれを離散バージョンにマッピングするということを行わず，
- 離散分布を推定
- その推定値を連続バージョンに当てはめる
- 連続バージョンの理論的な分位点$y_i^\tau$を各データ点について求める（$\tau$は分位点のレベル）
- 以下のリスクを最小化するように回帰モデルを推定する
E[L(y_i^\tau-h(\eta_i^\tau))]
ここで$L$は2乗損失などの損失関数，$h(\cdot)$は適当なリンク関数，$\eta_i^\tau$は例えば$\eta_i^\tau=x_i'\beta_\tau$といった回帰関数．
ということを行っています．
Mid-quantileアプローチ
Geraci and Farcomeni (2021)はMid-quantileによるある種の連続化のアプローチを考えました．聞き慣れないMid-quantileを導入するのに，Mid-分布関数（CDF）を以下のように定義します．
G_Y(y)=\Pr(Y\leq y)-0.5 \Pr(Y=y)
ここで$Y$は離散なので，$G_Y(y)$はステップ関数になります．もし$Y$が連続な場合は$\Pr(Y=y)=0$なのでこれは通常の分布関数になります．${y_1<\dots<y_s}$を$Y$が取りうる値とし，$S_Y$をそれらの集合とします．また$Y$はそれらの値をそれぞれ$p_1,\dots,p_s$の確率で取ることとします．またこれらの値に対して$\pi_1=p_1/2$，$\pi_j=G(y_j)=\sum_{h=1}^{j-1}p_h+p_j/2,\ j=2,\dots,s$を定義します．Mid-quantileは
H_Y(p)=\left\{
\begin{array}{ll}
y_1 & p<\pi_1\\
y_j & p=\pi_j,\quad j=1,\dots,s\\
(1-\gamma)y_j + \gamma y_{j+1} & p=(1-\gamma)\pi_j + \gamma \pi_{j+1},\quad 0<\gamma<1,\quad j=1,\dots,s-1\\
y_s & p>\pi_s\\
\end{array}
\right.
と定義されます．ところが，Mid-quantileは$y\notin S_Y$に対しては$G_Y(y)^{-1}$を得ることはできないので，$G_Y^c(y_j)=\pi_j$の間を線形補間した連続なバージョン$G_Y^c(y)$を考えることで，すべての$y\in\mathbb{R}$で$G_Y^c(y)=H^{-1}_Y(y)$とすることができます．
Mid-quantileのメリットとして
- $F_Y$が離散の場合に標本Mid-quantileは母集団Mid-quantileの一致推定量になり，標本分布が正規分布になる
- Mid-quantileはfractional order statisticsとして見ることができる
- $p=\pi_j$のとき，Mid-quantileは$Y$の分位点，そうでない場合は連続化されたものとして見られる
- 分位点との関係が複雑なexpectileなどと比べて，Mid-quantileから通常の分位点を簡単に復元できる
などが挙げられています．
分位点回帰を考える場合は，Mid-CDFは
G_{Y|X}(y|x)=F_{Y|X}(y|x)-0.5m_{Y|X}(y|x)
となります．ここで$F_{Y|X}(y|x)=\Pr(Y\leq y|X=x)$，$m_{Y|X}(y|x)=\Pr(Y= y|X=x)$．そして単調で微分可能なリンク$h(\cdot)$を通した$h(Y)$に対するMid-quantileを
H_{T(Y)|X}(\tau)=x'\beta_\tau
などといったようにモデル化します．回帰パラメータ$\beta_\tau$は
\sum_{i=1}^n\left(\tau-\hat{G}_{Y|X}^c(\eta_i|x_i)\right)^2
を最小化するように求められます．ここで$\eta_i=T^{-1}(x_i'\beta)$，$\hat{G}^c_{Y|X}$は$G_{Y|X}^c$の推定値．$\hat{G}^c_{Y|X}$に必要な条件付き分布$F_{Y|X}(y|x)$はカーネル推定によって推定します．ありがちな話ですが，$X$の次元が大きいときなどには条件付き分布の推定が困難になることが考えられます．
おわりに
通常の分位点回帰は一見カウントデータに対しても直接適用できそうな気がするかもしれないですが，意外と奥が深い問題があり，研究分野としても若干ニッチな感じがあります．分位点に対するトリートメントの効果を推定（quantile treatment effect; QTE）するためのアプローチとして分布回帰（distributional regression）が挙げられます（Chernozhukov et al., 2020）．分布回帰は被説名変数が連続，離散，混合に関わらず適用することできますが，分位点を直接モデル化したい場合には対応できません．今後の発展に期待です！
参考文献
- Chernozhukov, V., Fernandez-Val, I., Melly, B. and Wüthrich, K. (2020). Generic Inference on Quantile and Quantile Effect Functions for Discrete Outcomes. Journal of the American Statistical Association, 115, 123-137.
- Geraci, M. and Farcomeni, A. (2021). Mid-quantile regressino for discrete responses. arXiv:1907.01945.
- Lamarche, C., Shi, X. and Young, D.S. (2021). Conditional quantile functions for zero-inflated longitudinal count data. Econometrics and Statistics. https://doi.org/10.1016/j.ecosta.2021.09.003
- Machado, J.A.F. and Santos Silva, J.M.C. (2005). Quantiles for counts. Journal of the American Statistical Association, 100, 1226-1237.
- Padellini, T. and Rue, H. (2019). Model-aware quantile regression for discrete data. arXiv:1804.03714.
一緒にお仕事しませんか！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！
