Marra, G. & Radice, R. (2017). A joint regression modeling framework for analyzing bivariate binary data in R. Dependence Modeling, 5(1), 268-294.
先日から非確率標本で母集団を推定するという話についてあれこれ調べているんだけど、たいていの話は共変量の下での標本選択と目的変数の独立性、つまり欠損データ分析でいうところのMARを前提としていて、MNARモデルの話はなかなか出てこない。でも、そういうのもありそうじゃんね？ 考えてみたらHeckmanモデルとかそうじゃんね？
などと思いながらRパッケージを眺めていたら、GJRMパッケージというのが怪しげである。というわけで、実戦投入の予定はないが、vignetteに相当するであろう論文を読んでみた。掲載誌はみたこともない謎の雑誌だが、版元はDe Gruyterなので、そんなに変な雑誌ではないだろう。
1. イントロダクション
GJRMパッケージはたくさんの多変量反応回帰モデルを多様な標本抽出スキーマの下であてめる同時モデリング・フレームワークを実装しております。
主な関数は2つ。gjrm()で2変量回帰(ほか)、gamlss()で単変量回帰。
本論文では、以下にあてはまる2変量回帰に焦点を当てます: (1)内生性がある、(2)非ランダム標本抽出、(3)部分観察。二値変数のケースに絞ります。
本論文で検討するモデルの適用範囲は広いです。Chen & Astebro(2012 J.OR.Soc), Goldman, et al.(2001 JASA), Jeliazkov & Yang (2014 Book), Latif (2009, HealthEcon.), Li & Jensen (Inquiry), Pianzola(2014 Elect.Styd.), Shideler et al.(2015 Fish.Res.)あたりをみてね。
[なじみのない雑誌が多くてビビる。JASAの奴のabstractをみたところ、HIV感染者の死亡確率と保険の関係の研究で、単純に回帰すると保険に入っているほうが死亡率が高いんだけど、2変数の同時モデリングをやると保険のおかげで死亡率が下がっているのだそうだ。要するに道具変数をつかった研究じゃないかと思う]
1.1 内生性
未観察の交絡変数があるとき、処理は内生的だといわれる。そこでHeckman(1978)は2変量プロビットモデルを用いた。他にもいろんな方法があるけど似た結果になる(Clarke & Windmeijer, 2912 JASAをみよ)。
拡張はいろいろあるけど、Radice, Marra, Wojtys(2016)は非ガウス依存性と非線形的な共変量効果に拡張している[題名にcopulaとか書いてある… 私の手には負えなさそう…]。GJRMはこれを積んでます。
RパッケージならほかにVGAM, mvProbitパッケージがある。
1.2 非ランダム標本抽出
[…] 非ランダム標本選択バイアスへの対処としては選択モデルとパターン混合モデルがある。Fitzmaurice, et al.(2008 Book)をみよ。ここでは前者に焦点を当てる。アウトカムが二値の時は二変量プロビットモデルが使われる。Dubin & Rivers(1989 Sociol.MethodRes), Marra, et al.(2017 JASA), Marra & Radice(2013 ElectronJ.Stat.), McGovern et al.(2015 Epidemilogy)をみよ[最後のやつ、またcolupaとか書いてある…ううう…]。RだとsampleSelectionパッケージもある。
1.3 部分観察
ある二値アウトカムが2つの意思決定の選択の同時実現である場合を指す。二変量プロビットモデルでは4つの可能なアウトカムのうちひとつだけが観察されると仮定する。Poirier(1980 J.Econometrics)が提案した。
2. 方法論
\(i=1, \ldots, n\)について2つの二値変数\((Y_{i1}, Y_{i2})\)がある。変数\(j=1,2\)の標準化単変量分布(ここではガウシアンかロジスティックかガンベル)のcdfを\(F_j(\cdots)\), 加法予測子を\(\eta_{ji} \in \mathcal{R}\)として \(P(Y_{ij}=1) = 1-F_j(-\eta_{ji})\)とする。なんでこんな変な定義をするかというと、\(Y_{ij}\)を連続的潜在変数 \(Y^*{ij} = \eta_{ji} + \epsilon_i\)が0より大きいかどうかのインジケータとして捉えたいからで、つまり\(P(Y_{ij} = 1) = P(Y^*_{ij} \geq 0) = 1-F_j(-\eta_{ji})\) だからである。ほんとは周辺CDFは\(\eta_{1i}, \eta{2i}\)に条件づけられているんだけど略記する。
\(C\)をtwo-place copula関数、\(\theta_i\)を2つの確率変数の間の依存性を表すパラメータとし、以下のように定義する。$$ p_{11i} = P(Y_{i1} = 1 , Y_{i2} = 1) = C(P(Y_{i1} = 1), P(Y_{i2} = 1); \theta) $$ ここで\(\theta_i = m(\eta_{ci})\)とする。\(m\)は\(\theta_i\)がその範囲に入るように調整する1対1変換である。[\(\eta_{ci}\)ってなによ!? うぐぐぐ、さっぱりわからんが、とにかく線形予測子で決まるなにかなのであろう。コピュラについて全然知識がないので、そういうもんですか、という感想しかない…]
GJRMに搭載されているコピュラと\(m(\cdot)\)のリストはMarra & Radice(2017, Comput.Stat.DataAn.)に載せておいた。
標本の対数尤度関数はこうなる。\(a, b\)を二値として、\(y_{i1}=a, y_{i2} = b\)であることを示すインジケータを\(\mathbb{I}_{abi}\)として、$$ \ell = \sum_{i=1}^n \left\{ \sum_{a, b=0} \mathbb{I}_{abi} \log(p_{abi}) \right\} $$ [原文ではインジケータは \(1\)のmathbbなんだけど、Mathjaxでは出せないと思うので \(I\)に変更した。
えーと、要するに個体対数尤度は \(\log(P(Y_{i1} = y_{i1}, Y_{i2} = y_{i2}))\)だってことですよね？ ところで、上の\(P(Y_{i1} = 1 , Y_{i2} = 1) = C(P(Y_{i1} = 1), P(Y_{i2} = 1); \theta)\)という定義から、たとえば\(P(Y_{i1}=1, Y_{i2}=0)\)も決まるってこと？ よくわからんがこのコピュラというのは、単に2つの確率分布からその積の確率分布を出す関数というより、2つの周辺確率分布を同時確率分布にマップするような関数らしい…]
2.1 加法予測子の特徴づけ
\(1, 2, c\)をとる添字を\(v\)として、以下のように定義する。$$ \eta_{vi} = \beta_{v0} + \sum_{k=1}^{K_v} s_{v} (\mathbf{z}_{vki}) $$ [総和記号のインデクス\(k\)は原文では\(k_v\)なんだけど、読みにくいので変えた。だってそのまえに\(k\)が付いているから誤解しようがないじゃんね]
\(\mathbf{z}_{vi}\)は共変量ベクトルで、\(s_{vk}(\mathbf{z}_{vki})\)がそのジェネリックな効果を表す。その効果は\(Q_{vk}\)個ある基底関数 \(b_{vkq}(\mathbf{z}_{vki})\)の線形結合で近似できるとしよう[インデクス\(q\)は原文では\(q_{vk_v}\)なんだけど読みにくいので変えた]。すなわち、$$ s_{vk}(\mathbf{z}_{vki}) = \sum_{q = 1}^{Q_{vk}} \beta_{vkq} b_{vkq}(\mathbf{z}_{vki}) $$ これを全対象者について \(\mathbf{Z}_{vk} \beta_{vk} \)と書こう。ただし\(\mathbf{Z}\)は計画行列で\(\mathbf{Z}_{vk}[i, q_{vk}] = b_{vkq}(\mathbf{z}_{vki}) \)である[基底関数で変換済みということね]。とすると、全対象者について$$ \eta_v = \beta_{v0} \mathbf{1}_n + \mathbf{Z}_{v1} \beta_{v1}+ \ldots + \mathbf{Z}_{v K_v} \beta_{v K_v} $$ と書ける。さらに行列にまとめて$$ \eta_v = \mathbf{Z}_v \beta_v $$ とも書ける。
それぞれの\(\beta_{vk}\)は、フィッティングにおいては二次罰則項\(\lambda_{vk} \beta_{vk}^\top D_{vk} \beta_{vk}\)がかかることにしよう。[ひいいいいい。超めんどくせええええええ。おそらくこの話は本筋には影響しないと思うのでまるごとスキップする。正直つきあいきれねええええ]
以上の定式化によって多様な共分散の効果を表現できる。なお、GJRMパッケージではmgcvパッケージで利用できるすべての項を利用できる。
2.2 パラメータ推定
\(\delta = (\beta_1^\top, \beta_2^\top, \beta_c^\top)^\top, \mathbf{S} = diag(\mathbf{D}_1, \mathbf{D}_2, \mathbf{D}_c)\)として$$ \ell_p(\delta) = \ell(\delta) -\frac{1}{2} \delta^\top \mathbf{S} \delta$$を最大化する。[…]
2.3 その他の考慮事項
係数の信頼区間をどうやって出すかというと[…パス…]。
モデル構築には、コピュラ関数の選択、リンク関数のペアの選択、共変量選択が必要になる。AIC, BIC, 仮説検定の使用をお勧めします。
3. モデル
[ここが本題、というか最も関心ある箇所である。意外に短いぞ、頑張ろう]
ここからは、(1)内生性、(2)非ランダム標本選択、(3)部分観察をどう扱うかについて論じる。(1)(2)についてはRaddice, Marra, Wojtys(2006 Stat.Comput.), Marra, et al.(JASA)をみよ。
3.1 内生的処理を持つ二変量二値モデル
内生的処理を持つ二変量二値モデルは、主に、二値の処理が二値のアウトカムに及ぼす効果に関心があるんだけど未観察の交絡がある、というときに用いられる。
計量経済学では、この問題は「重要な共変量が無視されているせいでそれがモデルの誤差項の一部になっている回帰モデル」として定式化されることが多い。この文脈では、処理は観察された交絡因子で条件づけたときに誤差項と関連しないのなら外生的、そうでなければ内生的と呼ばれる。
二変量モデルは、未観察な交絡を構造的潜在変数フレームワークでコントロールしようとする。すなわち、１本目の方程式で二値アウトカム\(Y_{i2}\)を二値処理\(Y_{i1}\)の関数として表現し、2本目の方程式で処理を受けるかどうかを決める。で、2つの方程式の潜在誤差が連関パラメータ\(\theta_i\)を持つ二変量分布に従うと仮定する。加法予測子は $$ \eta_1 = \beta_{10} \mathbf{1}_n + \mathbf{Z}_{11} \beta_{11} + \ldots = \mathbf{Z}_1 \beta_1 $$ $$ \eta_2 = \beta_{20} \mathbf{1}_n + \beta_{21} \mathbf{y}_1 + \ldots = \mathbf{Z}_2 \beta_2 $$ $$ \eta_c = \beta_{c0} \mathbf{1}_n + \mathbf{Z}_{c1} \beta_{c1} + \ldots = \mathbf{Z}_c \beta_c $$ となる。対数尤度は$$ \ell = \sum_{i=1}^n \{ \mathbb{I}_{11i} \log(p_{11i}) + \mathbb{I}_{10i} \log(p_{10i}) + \mathbb{I}_{01i} \log(p_{01i}) + \mathbb{I}_{00i} \log(p_{00i}) $$ となる。
[なるほど。\(Y_{i2} = f(Y_{i1}, \mathbf{X}_i) + e_{i2}, Y_{i1} = g(\mathbf{X}_i) + e_{i1}, \mathbf{e}_i \sim hoge(\theta_i)\)という感じね。\(f, g\)が線形で\(hoge\)がMVNならMplusでも書けるモデルだ。これって共変量のなかに道具変数になる奴(どっちかにしか効かないとわかっている奴)が含まれていないといけないのだろうか、それとも完全に重複していてもどうにかなるのだろうか]
典型的には、主たる関心の対象は\(Y_{i2} = 1\)の確率に対する\(Y_{i1}\)の効果である。その指標はいろいろあるが、ここでは母集団における平均処理効果ではなく、手元の標本における平均処理効果$$ SATE(\beta_2, \mathbf{Z}_{2i}) = \frac{1}{n} \sum_{i=1}^n \{ (P(Y_{i2} = 1 | Y_{i1} = 1) – P(Y_{i2} = 1 | Y_{i1} = 0)\} $$ を用いる。なお、イベント生起確率への効果を出す関数もちゃんと用意してますよ。
3.2 非ランダム標本選択を伴う二変量二値モデル
非ランダム標本選択を扱うひとつの方法は構造的潜在変数フレームワークである。すなわち、１本目の方程式で標本選択されたときの二値アウトカム\(Y_{i2}\)を、2本目の方程式で標本選択\(Y_{i1}\)をモデル化し、2つの方程式の潜在誤差が連関パラメータ\(\theta_i\)を持つ二変量分布に従うと仮定する。
加法予測子は、$$ \eta_1 = \beta_{10} \mathbf{1}_n + \mathbf{Z}_{11} \beta_{11} + \ldots = \mathbf{Z}_1 \beta_1 $$ $$ \eta_2 = \beta_{20} \mathbf{1}_{n_s} + \mathbf{Z}_{21} \beta_{21} + \ldots = \mathbf{Z}_2 \beta_2 $$ $$ \eta_c = \beta_{c0} \mathbf{1}_{n_s} + \mathbf{Z}_{c1} \beta_{c1} + \ldots = \mathbf{Z}_c \beta_c $$ となる。\(\eta_1\)は\(n\)人について求めるが、\(\eta_2, \eta_c\)は選ばれた\(n_c\)人に対してしか求めない。
対数尤度関数は $$ \ell = \sum_{i=1}^n \{ \mathcal{I}_{11i} \log(p_{11i}) + \mathcal{I}_{10i} \log(p_{10i}) + (1-y_{i1}) \log(p_{0i}) \} $$ となる。
[おーっと…!!! つまり、\(n\)は確率標本のサイズってことだよね。非確率選択って、「サイズ\(n\)の確率標本からさらに非確率選択が起きて\(n_s\)人になった」という場合(確率調査における無回答)、「非確率標本から確率抽出して\(n_s\)人になった」という場合(確率調査における台帳不備)、「母集団からなんだかよくわからん標本\(n_s\)が得られた」という場合(任意参加型web調査)、などがあるけれど、ここでは最初のケースを想定しているわけだ。
でも、適用範囲が広いのは3番目のケースだと思う。2番目や3番目のケースをこのモデルで扱う場合、\(n\)はなにに相当するの？ 母集団サイズ？ ううう、きっとそうだな。その場合は、\(\eta_2, \eta_c\)のモデルが組めるくらいの共変量が、母集団の全個体について測定されていないといけないことになる。辛いなあ…
むしろ、非確率標本の傾向スコア調整のときのように、このモデルを非確率標本と参照確率標本の結合標本から推定できたりはしないのだろうか？ その場合は\(n\)は結合標本のサイズだよね。でも、その場合\(w_i\)はどうすればいいの？ 仮に適切な\(w_i\)があるとして、そのとき結合標本でのPREVは真のPREVの一致推定量になるの？ 結局、非確率標本の傾向スコア調整の場合と同じ問題になりそうだ。ううう、辛い…]
推定の対象は母集団における有病率とかであろう。これは次のように推定できる: 標本ウェイトを\(w_i\)として、$$ PREV (\hat{\beta}_2, \mathbf{Z}_2) = \frac{\sum_{i =1}^n w_i \{ 1-F_2(\hat{\eta}_{2i}) \}}{\sum_{i=1}^n w_i} $$ [びびるな、単に\(P(Y_{ij}=1)\)のHijek推定だ]
なお、標本選択モデルではふつうモデル識別のための除外制約が必要になる。すなわち、選択を予測するがアウトカムを予測する変数が必要になる。
[標本選択の道具変数が必要なわけだ。ということは、3.1節のケース、すなわち単なる内生性を扱いたい場合には、道具変数はいらないということか…]
3.3 部分観察を扱う二変量プロビットモデル
2つの二値変数\(Y_{i1}, Y_{i2}\)があって\(Y_{i1}Y_{i2}\)しか観察できない場合、加法予測子は $$ \eta_1 = \beta_{10} \mathbf{1}_n + \mathbf{Z}_{11} \beta_{11} + \ldots = \mathbf{Z}_1 \beta_1 $$ $$ \eta_2 = \beta_{20} \mathbf{1}_n + \mathbf{Z}_2 \beta_{21} + \ldots = \mathbf{Z}_2 \beta_2 $$ $$ \eta_c = \beta_{c0} \mathbf{1}_n + \mathbf{Z}_{c1} \beta_{c1} + \ldots = \mathbf{Z}_c \beta_c $$ となり、対数尤度は$$ \ell = \sum_{i=1}^n \{ \mathbb{I}_{11i} \log(p_{11i}) + (1-\mathbb{I}_{11i}) \log(1-p_{11i}) \} $$ となる。[なるほど。モデルは変わらん、対数尤度関数が変わる、ってことね]
関心の対象は\(p_{11i}\)の推定値、ならびにそれらに対する共変量の効果であろう。
なお、このモデルはガウシアン周辺分布とガウシアン・コピュラによって定義される。[どういうこと? このモデルに限り、\(\eta_{1}, \eta_{2}\)の誤差項がMVNだという仮定が必要になるということ？]
上の対数尤度関数は非線形なので、ふつう局所識別の問題が生じる。また、\(\eta_1, \eta_2\)を入れ替えても等価なモデルが生じるので(ラベリング問題)、ふつうは共変量にすくなくともひとつの除外制約をつける。2つのアウトカムに影響する未観察変数の相関がないならば、\(\theta = 0\)と仮定してモデルを単純化できる。このとき\(p_{11i} = (1-F_1(-\eta_{1i}))(1-F_2(-\eta_{2i}))\)となる。
4. GRJMパッケージの gjrm()関数
[使用方法と使用例。これは実際に使うときに読めばいいや…]
5. 事例
5.1 健康サービスへの利用に対する民間健康保険の効果
[内生性の問題であろう。スキップ]
5.2 HIVの有病率
ザンビアの2007年のHIVについてのデータセットを使う。HIV有病率を推定したいんだけど、HIVの検査は世帯調査の一環として行われ、有病者は罹患がばれるのを恐れて調査に参加しない傾向がある。代入ベースの方法のようなよくあるアプローチでは補正できない。さらに、どの下位集団に介入すべきかを知りたいので、HIVの予測子とか空間的依存性とかを取り入れてモデル化したい。
調査では、最後の項目でHIV検査への同意をとり、同意者の血液を採取した。
[ああ、なるほど、やっぱり「確率調査からの非確率選択」というケースだ。腑に落ちたので、残りはスキップ]
5.3 市民戦争開始の決定要因
市民戦争は反政府派と政府との相互作用の結果として理論化されることが多い。つまり、両方の意思決定者の同時決定だけが観察される。そこで[… なるほど、こういう風に使うモデルなのか。腑に落ちたのであとはスキップ]
6. 考察
[…]
GJRMパッケージではリンク関数はprobitで、他の関数についても調べたけど結局載せていない。数値的に困難だし結果はたいして変わんないからだ。交換不能コピュラも同様の理由で載せていない。
最初に述べたけど、多変量反応もモデル化できるからね。
云々。
——
正直、コピュラについての知識が全くなく、コピュラが出てくる論文はこれまで避けていたもので、2節でコピュラが出てきたときには「やられた、騙された…」という感じであった。誰も騙してないってば。
ま、細かいところはよくわからんが、どういう場面で使うもんなのかなんとなく見当がついたので、これでよしとしよう。
web調査の標本選択バイアスを修正したい場合は、傾向スコア調整やふつうのアウトカム・モデリングとは異なり、標本選択がMNARでもオッケー。でも共変量は母集団の全個体について測定されていないといけないし(いいかえると全共変量の同時分布がわかんないといけない)、そのなかに良い道具変数がないといけない。ということだろうと思う。