« 読了：木山(2017) RCT至上主義 (in 開発経済学) とその問題 | メイン | 読了：尾崎(1993) エリア・マーケティングの理論と実際 in 1993 »
2018年6月 8日 (金)
Min, Y., Agresti, A. (2002) Modeling nonnegative data with clumping at zero: A survey. Journal of The Iranian Statistical Society, 1, 7-33.
仕事の都合で知りたいことがあって、タイトルがまさにジャスト・ミートだったもので読んだ論文。
掲載誌は... イラン統計学会...?! なぜにイラン?! とびびりましたが、Google Scholar上では被引用回数94, 第二著者はカテゴリカル・データ解析の神様アグレスティ先生、読んでも損はなかろうかと。
タイトルにあるClumping at Zeroという表現、なんと訳したらいいのかよくわからないんだけど、要はゼロがいっぱいあるということであろう。
1. イントロダクション
連続的分布を持っているけど値0の確率質量だけが不連続。こういうのを半連続変数と呼ぼう。これに関連するタイプのデータとしてゼロ過剰カウントデータがある。なお、これらを左打ち切り(censored)データと混同しないように。ここで値0は本物のアウトカムである。
半連続データやゼロ過剰カウントデータを扱う際の難しさは、正規分布やガンマ分布を当てはめるわけにいかないという点にある。本論文は、このような「ゼロに集まっている」データをモデリングするための方法についてレビューする。
2. 半連続データのモデル
初期の研究は主に計量経済学でなされた。Tobin(1958)は耐久財の世帯支出を記述するための打ち切りつき回帰モデルを提案した。いまでいうところのTobitモデルである。この系統のモデルは正規確率変数がなんらかのメカニズムで打ち切られると考える。
このほかに、正規分布を仮定しない系統の研究もある。
2.1 Tobitモデル
反応変数を$Y$とする。被験者$i$の観察を$y_i$とする。
Tobitモデルではこう仮定する。正規分布に従う変数$Y^*_i$を考えて、それが0より大なら$y_i = y^*_i$、でなければ$y_i = 0$。
説明変数ベクトル$x_i$があるならば
$y^*_i = x'_i \beta + u_i$
ただし$u_i$はiidに$N(0, \sigma^2)$に従う。
TobitモデルはML推定できる。$N(0,1)$のCDFを$\Phi(\cdot)$, PDFを$\phi(\cdot)$として、Tobitモデルの尤度関数は...[略]。
セミパラ推定の提案もある...[略]。そっちのほうが性能がいいという研究がある。
2.2 Two-partモデル
Tobitモデルでは、反応がゼロか正値を決める確率的過程と、反応が正値であるときにその値を決める確率的過程とが同じである。ここからはそうでないモデル。
Duan et al.(1983)はこう考えた。
$logit(P(Y_i=0)) = x'_{1i} \beta_1$
$log(y_i | y_i > 0) = x'_{2i} \beta_2 + \epsilon_i$
ただし$\epsilon_i \sim N(0, \sigma^2)$。
これはわりかし簡単にML推定できる、というか、別々に推定すればいいのである。尤度関数は...[略]。
このモデルは、アウトカムが本当に正である時のアウトカムのレベルを記述する、条件付きの式を推定している。
2.3 標本選択モデル
Heckman(1974, 1979)の提案も、Tobitモデルをtwo-partモデルに拡張することであった。その後いろんなバージョンが出てきたんだけど、ここではvan de Ven & van Praag (1981)に基づいて説明しよう。
観察$i$について、$\{u_{1i}, u_{2i}\}$がiidに二変量正規分布$N(0, \Sigma)$に従うとする。$\Sigma$の対角成分を$\sigma^2_1, \sigma^2_2$とし、非対角を$\sigma_{12}$とする。
モデルは以下の通り。
$I_i = x'_{1i} \beta_1 + u_{1i}$
$y^*_i = x'_{2i} \beta_2 + u_{2i}$
で、$y_i$は$I_i > 0$のときに$\exp(y^*_i)$、そうでないときに$0$になる。
このモデルはML推定できる。尤度関数は...[略]。
Heckmanの2段階推定でも推定できる...[説明略]。MLより良いとはいえないけど、すっごくかんたんなので、標準的な推定方法となっている。
もっとも、Duan et al.(1983 J.Business&Econ.Stat., 1984同誌)いわく、このモデルは数値的にも統計的にもあまりよくない。尤度関数は局所解に落ちるかもしれないし、計算も大変。打ち切られたデータは観察不能だという想定を置いているわけだけど、それを検証する方法がない。2段階推定だと、結局、下位標本について
$\log(Y_i) = x'_{2i} \beta_2 + (\sigma_{12}/\sigma_{1})\lambda_i + \epsilon_i$
ただし$\lambda_i = \phi(z_i)/\Phi(z_i), z_i = x'_{1i} \beta_i/\sigma_1$という式になるんだけど[ああ...逆ミルズ比っていうんだっけか...]、$\lambda_i$が$x_2$と強く相関していたら推定量が不安定になる。$x_1$と$x_2$は別の変数にしなさいとよく言われるが、そんなの現実的じゃない。
このモデルは、「もしすべての人がアウトカムを持っていたらそのレベルはなんだったか」を記述する、条件付きでない式を推定している。two-partモデルとの違いにご注目。どっちがいいかという点については議論がある。上述のDuan et al.や, Manning et al.(1987 J.Econometrics), Leung & Yu (1996 J.Econometrics)をみるがよろしい。
2.4 Compound Poisson exponential dispersionモデル
一般化線型モデルでは指数型分布族
$f(y_i; \theta_i, \phi) = c(y_i, \phi) \exp(\frac{\theta_i y_i - b(\theta_i)}{\phi})$
を使う。平均$\mu_i$, その分散を$v(\mu_i)$として、$v(\mu)=\mu$とするとポアソン分布である。[ううう...この辺、良く理解できていない...]
Jorgensen(1987,1997)はこう考えた。$v(\mu)=\mu^p$かつ$1 < p <2$とするとこれは複合ポワソン分布になって...なんたらかんたら... なんたらかんたら...
[本節、私の能力を超えて大空を飛んでいったので、ぼんやり空を見上げたまま飛び去るのを待ちますが、要するに、半連続変数をひとつの確率分布でうまく表現できるような謎の確率分布を考える、というアプローチらしい。なんだか知らんが推定できるらしい。Tweedie分布という言い回しを聞いたことがあるけど(SASのGENMODとかで)、これがそれかしらん。勉強する気もあまりないんだけど、Rのcplmパッケージというのがこれなんじゃないかと思う]
2.5 順序閾値モデル
Saei, Ward, & McGilchrist(1996)は、反応変数を順序カテゴリ変数に落とせばいいじゃんと考えた。順序カテゴリ変数$Y_g$の値が、連続潜在変数$Z$と固定閾値で決まると考えるわけね。で、最初のカテゴリを値ゼロを表すカテゴリだってことにすりゃあいい。あとは$Z$のCDFを好きに決めるがよい。ロジスティックとか正規とか。
two-partモデルとは異なり係数がひとつなので、群間でモデル比較しやすい、というのがメリット。カテゴリに落とすときの恣意性と情報損失がデメリット。
2.6 既存アプローチの得失
ある種の打ち切り反応変数であれば、Tobitモデルや標本選択モデルが向いているだろう。でも、ゼロがほんとのアウトカムなら、たいていの場合正規性仮定に無理があるから、two-partモデルが向いている。合成ポアソン指数型分布モデルも悪くないけど(単一モデルで分析できるという意味ではシンプルだともいえる)、あんまり使われてない。順序反応モデルは簡単でいいけど元データを分析してないというのが欠点。
というわけで、いちばん適用範囲が広いのはtwo-partモデルだといえましょう。
3. ゼロ過剰カウントデータのモデル
ゼロが過剰なカウントデータの場合、ただのポアソン分布だと非ゼロの値の変動を説明できない。負の二項回帰モデル(これはポアソン平均のガンマ混合)なり、ポアソン平均の対数のモデルに正規誤差をいれるなりで対処することもあるけど(Agresti本の13章をみよ)、うまくいかないことが多い。たとえば、分布がゼロとどこかの二峰になっている場合、やはり混合分布を考えるほうが自然だろう。
3.1 ゼロ過剰離散分布
Lambert(1992)はゼロ過剰ポアソン回帰モデル(ZIP回帰モデル)を考えた。このモデルでは、データをゼロとポアソン変量のアウトカムとの混合とみなす。つまり、$Y_i$は確率$p_i$で$0$となり、そうでない場合に$Poisson(\lambda_i)$に従う。二値潜在クラス$Z_i$があるのだと考えても良い。でもって、
$logit(p_i) = x'_{1i} \beta_1$
$\log(\lambda_i) = x'_{2i} \beta_2$
とする。
尤度関数は...[略]。EMアルゴリズムでML推定できる。
Hall(2000)は、$Poisson(\lambda_i)$じゃなくて$binomial(n_i, \pi_i)$と考えた。こうすると$Y_i$の上界が$n_i$になる。でもって、
$logit(p_i) = x'_{1i} \beta_1$
$\log(\pi_i) = x'_{2i} \beta_2$
とする。これもEMアルゴリズムでML推定できる。
実際のところ、カウントデータにはoverdispersionがつきものである。非ゼロかどうかとか潜在クラスとかで条件づけてもやっぱりそうである。ZIPモデルだと、$Z_i=0$だったら全員同じと考えているわけで、これはやっぱり無理があることが多い。その意味ではZIPモデルよりゼロ過剰な負の二項分布モデルのほうがよい面がある。どっちがよいかを決めるスコア検定というのもある。
すごく確率が高い値が複数個あるときに注意。たとえば女性の出産数は0か2が多いので、多項ロジットモデルで解いた研究がある。
3.2 ハードルモデル
これはMullahy(1986)が提案したtwo-partモデル。ひとつは非ゼロかどうかを説明する二値モデルで、これがハードル。もうひとつは、ハードルを越えたという条件の下でアウトカムのレベルを説明するモデルなんだけど、ここで切断(truncate)を考える。切断つきポアソンとか、切断つき負の二項とか。
たとえばこう考える。[説明がないけど$\lambda$というのはポアソンパラメータであろう]
$logit[P(Y_i=0)] = x'_{1i} \beta_i$
$\log(\lambda_i) = x'_{2i} \beta_i$
すると尤度関数は...[略]。
カウントの右裾がすごく長いこともある。そういう場合のためのセミパラモデルの提案もある。
3.3 有限混合モデル
3.1に述べたように、ゼロ過剰離散分布は一種の有限混合分布である。これをもっと緩くして、ゼロの母集団と非ゼロの母集団をはっきり区別しないようにすることもできる... [以下メモは省略するけど、なるほどね、正面からゼロ過剰という問題に取り組むのではなくて、ふつうの潜在クラスポワソン回帰モデルを組むわけね。Wedel, DeSarbo, Bult, & Ramaswamy(1993, J.App.Econ.)というのを読むのがよさそう]
3.4 Neyman type A分布
ゼロ過剰カウントデータを、Neyman type A分布というのをつかってモデル化するという提案がある。これはポアソン-ポアソン混合分布で...[この節、一読したんだだけど、実現値$Y_i$はあるポアソン分布に独立に従う$N_i$個の変数の和、$N_i$は別のポアソン分布に従う、という悪い冗談のような話であった。人生においてそんなややこしい話と関わり合いを持ちたくないので、パス]
3.5 既存アプローチの得失
ゼロ過剰モデルとハードルモデルは似ている。母集団は反応がゼロの人々と非ゼロの人々の混合だ、という考え方が自然ならばゼロ過剰モデルがよい。いっぽうハードルモデルはシンプル。さらに、ゼロが過剰なんじゃなくて、標準的な仮定に基づく期待よりも少ないという場合にも、ハードルモデルは使える。
有限混合モデルは、観察が異なる母集団から抽出されているとみなすのが現実的ならば魅力的。ただし、モデルが適合してないときクラス数が過大評価されるかもしれない。Newman type Aモデルは推定が大変。
4. zero-clumpedな反復測定データのモデリング
ここまでは横断データの話であった。ここからは、縦断データみたいにクラスタのあるデータ(観察間に相関のあるデータ)について。研究は多くない。
4.1 半連続データの反復測定
Tobitモデルや標本選択モデルを縦断データに拡張しようという研究はあるんだけど、ゼロが本物のアウトカムだった場合には正規性仮定に無理があるので、本論文では扱わない。
Olsen & Schafer(2001)は、Duanらのtwo-partモデルを縦断データに拡張している。被験者(ないしクラスタ)$i$, 時点$j$の半連続反応を$y_{ij}$、それが0である確率を$p_{ij}$とする。まずゼロか正値かについてのロジスティック・ランダム効果モデルを組む。
$logit(p_{ij}) = x'_{1ij} \beta_1 + z'_{1ij} c_i$
$x_{1ij}$と$z_{1ij}$は共変量、$\beta_1$が固定係数、$c_i$がランダム係数。$c_i$はベクトルにしておくけど、通常はランダム切片モデルで十分であろう($c_i$はスカラーで$z_{1ij}$はなし)。
で、$y_{ij} > 0$のとき$V_{ij} = y_{ij}$として、
$\log(V_{ij}) = x'_{2ij} \beta_2 + z'_{2ij} d_i + \epsilon_{ij}$
ただし$\epsilon_{ij} \sim N(0, \sigma^2)$。ここでも、通常はランダム切片モデルで十分であろう($d_i$はスカラーで$z'_{2ij}$はなし)。
$c_i$と$d_i$を縦に繋いで$b_i$とし、$N(0, \Sigma)$に従うということにする。で、ここが工夫なんだけど、$\Sigma$は対角行列じゃなくて、全要素を推定するわけ。
Olsenらは推定手法をいろいろ比べている。MCMCとかEMとか罰則付き擬似尤度とかGH求積とか高次ラプラス近似とか。[ああ、そうか...よくわかんないけど、これって一般化線形混合モデルに近い話なのではなかろうか...]
順序閾値モデルを拡張するというアイデアもある。Saei, et al.(1996)は...[略]
4.2 ゼロ過剰データの反復測定
Hall(2000)はゼロ過剰ポアソンモデルを時系列に拡張している。
例によって、$Y_{ij}$は確率$p_{ij}$で0、そうでなかったら$Poisson(\lambda_{ij})$と考える。で、
$logit(p_{ij}) = x'_{1ij}\beta_1$
$log(\lambda_{ij}) = x'_{2ij} \beta_2 + b_i$
ただし$b_i \sim N(0, \sigma^2)$と考える。尤度は...[略]。推定は...[略]。なお、ポワソンじゃなくて$binomial(n_i, \pi_{ij})$とし、二本目で$logit(pi_{ij})$のモデルを組むという手もある。
一本目のほうにはランダム効果が入っていない点にご注目。
You & Lee(2001)はハードルモデルにランダム効果をいれるというのを提案していて...[略]。
5. clumpが2つあるデータへの適用
たとえばゼロと最大値がやたらに多いというような場合のように、clumpがふたつある場合にはどうするか。これは患者の順守(compliance)の研究で問題になる。全然いうことを聞かない患者と、完璧に言うとおりにする患者が多いから。ロジット変換して0-1データを実数データにしたところで、clumpに対処したことにはならない。というわけで、我々はいまこの問題を研究してます。
6. 今後の課題
まず、時系列データの研究が足りない。[いくつか「あのモデルをこう変えるのはどうか」的なことが書いてある。略]
半連続データの手法はたいがい、正値であれば対数正規だと仮定している。でもたとえば医療費はもっと右に裾を引く。セミパラな手法のほうがいいかもしれない。
云々。
。。。やれやれ、読み終えたぞ。
そんなに重量級のレビューって感じではなくて、修論かなにかの1章を切り出したような感じだけど、勉強になりました。セレクション・バイアスに対するヘキットみたいなアプローチと、ゼロ過剰カウントデータに対するアプローチはたいてい別の文脈で説明されていると思う。私のような初学者は学習時の文脈に過度に制約されちゃうきらいがあるので、こういう風に、多様な文脈にまたがる話を統一的な記法でレビューしてもらえると、すごくすっきりする。ありがとう著者のみなさん、ありがとうイランのみなさん。
半連続データのところ、意外にもヘキットでなくてTwo-partモデル推しであった。考えてみると、データが半連続であるということの背後にある発生メカニズムを特定せず、一般的にレビューすれば、それはそうなるんだろうなと思う。ヘキットというのは誤差に対して結構きつい仮定を置いているわけで、どこでも使えるわけじゃない。やっぱし、いろんなツールについて知っていることと、目の前の現象について都度都度真剣に考えること、この両方がないといかんのだなあ...と、これは反省の弁である。
論文：データ解析(2018-) - 読了：Min & Agresti (2002) ゼロがやたらに多い非負データの分析方法レビュー