elsur.jpn.org >

« 2017年10月 | メイン | 2017年12月 »

2017年11月26日 (日)

Frank, M.R., Cebrian, M., Pickard, G., Rahwan, I. (2017) Validating Bayesian truth serum in large-scale online human experiments. PRoS ONE. 12(5).
　原稿の準備で読んだ奴。読んだ際のメモが出てきたので記録しておく。久々のベイジアン自白剤論文で、面白く読んだという記憶がある。

　第1著者はMITメディアラボの人で、Prelecとどういうつながりがあるのかわからない(謝辞にPrelecの名前はない)。第3著者はGoogle所属。

　いわく。
　調査回答者に主観的判断を求めるということが、各分野においていかに不可欠か、という前置きがあって...
　不誠実な回答を引き起こす原因のひとつは強欲である。特にAmazon Mechanical Turk(MTurk)なんかだと回答者は利益の最大化を目指すわけで、これは深刻な問題になる。
　対処策としてベイジアン自白剤(BTS)が提案されているけれど、実証実験は小規模なのしかない。そこで大規模にやりました。

　BTSの説明。
　BTSとは、正直さ、ないし得られた情報に応じて報酬を与える方法で...[←という風に、BTSを明確にrewardingの手法として紹介している。この辺は書き手によってニュアンスが違うところだ]
　その仕組みは...[中略]...まあそういうわけで、α > 0で正直さがベイジアン・ナッシュ均衡になり、α=1でゼロサムゲームになる。本研究ではα=1とする。
　
　実験。MTurkでやる。
　以下、実験群には「情報スコアが上位1/3にはいったら追加ボーナスを金で払う」と教示。しかし情報スコアの中味は教えず、かわりに「MITの研究者が開発した真実申告検出メカニズムだ」と教示する。
　実験群は2種類。(1)透過BTS群。回答から情報スコアを動的に算出して提示。(2)BTS intimidation群。情報スコアは見せないが報酬は渡す。[恫喝群とでも訳すところか]

　結果。

　考察。
　BTSによる改善は、統制群よりも報酬の期待値が大きいせいか。先行研究によれば、金銭的インセンティブの増大は作業量の増大を招くが作業の質は増大させない(Mason & Watts, 2020 ACM SigKDD Newsletter)。本研究でもそうで、統制群の報酬を増やしたけど結果はかわんなかった(補足資料をみよ)。
　調査における回答の正直さ促進の手法として、honesty pledgeとか、宗教的正直さの喚起とかを行う手法があるけど、きっとこの実験の恫喝群でも同じ事が起きたのだろう。つまり同じ効果が、怒れる神とか個人的誠実性の喪失とかへの恐怖から得られたかもしれない[おおっと... BTSが一種のbogus pipelineである可能性を認めちゃうのね...]。いっぽう透過群では回答分布がさらに正直な方向に変わった。
　云々。

　...小声で超偉そうな言い方をしちゃうと、わかりやすくよく書けている論文である(うわあ、何様だろうか)。PLoS ONEだからといってなめてはいけない。ちょっと図表が冗長な感じだがな(すいませんすいません)。
　この実験、統制群と実験群の比較じゃなくて、恫喝群と透過群の比較が一番面白いところだと思うんだけど、見た感じではそんなに明確な差じゃない気がする。

　イントロのところからメモ:

よくみると、やたらにWattsの論文を引用している。

論文：予測市場 - 読了：Frank, et al. (2017) ベイジアン自白剤 in クラウド・ソーシング

2017年11月16日 (木)

Ghosh, M., Rao, J.N.K. (1994) Small Area Estimation: An Appraisal. Statistical Science, 9(1), 55-93.
　タイトルの通り、小地域推定についてのレビュー。仕事の都合で泣きながら徹夜で読んだ。辛い。というかさ、小地域の推定なんて諦めようよ... 標本が小さいんだから...

1. イントロダクション
　小地域推定の歴史とか、活用場面の紹介とか。全部すっ飛ばして...

2. 人口学的方法
　人口統計学で用いられてきた、一番古典的な方法。行政からもらってきた、出生とか死亡とか新築とか入学とかのデータ(これをsymptomtic変数という)を、センサスのデータと併せて使う。
　そのひとつがVital Rates(VR)法。直近のセンサス年$t=0$について、小地域の粗出生率を$r_{10}$, 粗死亡率を$r_{20}$、その上の親地域の粗出生率を$R_{10}$, 粗死亡率を$R_{20}$とする。年$t$の小地域の粗出生率$r_{1t}$と粗死亡率$r_{2t}$をそれぞれ次とする：
　$r_{1t} = r_{10}(R_{1t}/R_{10})$
　$r_{2t} = r_{20}(R_{2t}/R_{20})$
　[えーと、つまり今年の出生率なり死亡率なりのセンサス年に対する比が小地域と親地域で同じだとするわけね]
　すると、年$t$の小地域について粗出生数$b_t$、粗死亡数$d_t$がわかっているとき、人口を次のように推定できる：
　$P_t = (1/2) (b_t / r_{1t} + d_t / r_{2t})$
　[ああそうか、粗出生率とか粗死亡率ってのは出生数とか死亡数を人口で割ったものだからね]
　この手法の弱点は、いうまでもなく、$r_{1t}/r_{10}$とかを$R_{1t}/R_{10}$とかで近似しちゃうところ。
　これを拡張した手法として、composite法とか、米センサス局のCM-II法とか、Administrative Records法とかHousing Unit法とかがあって...[なんとなく雰囲気がわかったので省略]

　こういう手法はすべて重回帰の特殊ケースとみることができる。そのほかに、symptomatic変数を独立変数にとった重回帰という手もある。
　そのひとつがratio-correlation法。前々回センサス年を$\alpha=0$, 前回センサス年を$\alpha=1$、今年を$\alpha=t(>1)$としよう。小地域$i$の年$\alpha$における人口を$P_{i \alpha}$, $j$番目のsymptomatic変数を$S_{ij\alpha}$とする。
　まず全部割合にする[うわ...気色悪い...]。
　$p_{i\alpha} = P_{i \alpha} / \sum_i P_{i \alpha}$
　$s_{ij\alpha} = S_{ij \alpha} / \sum_i S_{ij \alpha}$
センサス年との比、センサス年同士の比をとる。
　$R'_{i} = p_{i1} / p_{i0}$
　$R_{i} = p_{it} / p_{i1}$
　$r'_{ij} = s_{ij1} / s_{ij0}$
　$r_{ij} = s_{ijt} / s_{ij1}$
まず、センサス年の間の人口の伸びを目的変数にして回帰する。[めんどくさいので略記するけど$\beta$にはみんなハットがついている]
　$R'_{i} = \beta'_0 + \beta'_1 r'_{i1} + \ldots + \beta'_j r'_{ij}$
この係数でセンサス年との比を当てに行く。
　$\tilde{R}_{i} = \beta'_0 + \beta'_1 r_{i1} + \ldots + \beta'_j r_{ij}$
で、今年の全地域人口(これは別途わかっているとしよう)に$p_{i1}$を掛け、さらに$\tilde{R}_{i}$もかけて、小地域人口の推定値とする。
　この手法の弱点は、センサス年の間で作った$\beta'_j$が、今年と直近センサス年の間でも変わらないという仮定である。

　そこで出てきたのが標本回帰法。
　$m$個の小地域のうち$k$個については、センサス年に対する比$R_i$の標本推定値$\hat{R}_1, \cdots, \hat{R}_k$があるとしよう。[←なるほどね、ありそうな話だ。小地域のなかにもでかい奴あるしね]
　これを説明する回帰モデルを組む。[例によって$\beta$の上のハットは略記]
　$\hat{R}_i = \beta_0 + \beta_1 r_{i1} + \ldots + \beta_p r_{ip}$
このモデルを全ての小地域にあてはめる。

3. 合成推定量(synthetic estimators)と関連する推定量
　要するに、大地域で得た推定量を使っちゃえという方法。この路線も歴史が古い。

　人口が大きな領域$g$に分割されている[えーと、性別みたいなもんですかね]。で、各領域の合計$Y_{.g}$については信頼できる推定量$\hat{Y}'_{.g}$がある。で、それとは別に人口が小地域$i$に分かれてて、セル$(i,g)$の合計$Y_{ij}$を小地域を通じて足しあげると$Y_{.g}$になる。で、なんらかの補足情報$X_{ig}$も手に入るとしよう。
　いま、$Y_{i} = \sum_g Y_{ig}$を推定したい。そこで
　$\hat{Y}^S_i = \sum_g (X_{ig}/X_{.g}) \hat{Y}'_{.g}$
とすればいいんじゃないかしら。というのが最初期の合成推定量である。
　[以下、推定量の分散とかバイアスとかの話が続くので端折って...]
　Purcell & Kish (1980)とかがSPREE (構造保存推定)というのを提案していて、これは多元クロス表における反復比例フィッティングを使っているのだけれど、このSPREEは合成推定量を一般化したものだと捉えることができる。[←うおおお。SPREEってKishの論文に出てきて意味不明だったんだけど、そういう意味だったのか...]

　[ここから関心ある箇所なので細かくメモ]
　合成推定量には潜在的バイアスがあり、いっぽう直接推定量は不安定なので、重みづけ平均を取ろうという路線もある。
　つまりこうだ。小地域$i$の直接推定量を$\hat{Y}_{1i}$, なんらかの間接推定量を$\hat{Y}_{2i}$として、適切な重み $w_i$ (0以上1以下)を用意して
　$\hat{Y}^C = w_i \hat{Y}_{1i} (1-w_i) \hat{Y}_{2i}$
　こういうアプローチは多いんだけど、ここでは直接推定量が不偏推定量$\hat{Y}_i$、間接推定量が合成推定量$\hat{Y}^S_i$であるときの最適な重み$w_i$(opt)の決め方について考える。

　$w_i$(opt)を$\hat{Y}^C$のMSEを最小化される重みと捉えると、ふたつの推定量の共分散が0だとして、
　$w_i$ (opt) $= MSE(\hat{Y}^C) / [MSE(\hat{Y}^C)+V(\hat{Y}_i)]$
となる。これは推定できるんだけどすごく不安定である。
　Purcell & Kish (1979)は、MSEの平均を最小化する共通の$w$を使うというのを提案している。すると次の形になる：
　$w_i$ (opt) $= 1 - \sum v(\hat{Y}_i) / \sum (\hat{Y}_i^S - \hat{Y}_i)^2$
もし$\hat{Y}_i$の分散がだいたい等しかったら、$v(\hat{Y}_i)$を$\bar{v} = \sum v(\hat{Y}_i)$で置き換えて、
　$w_i$ (opt) $= 1 - m\bar{v} / \sum (\hat{Y}_i^S - \hat{Y}_i)^2$
というわけで、James-Steinタイプのウェイトになる。[←ちょ、ちょっと待って？！ これがなぜJames-Steinタイプなの？ というかタイプってなに？ どこまでがタイプなの？]
　もっとも、個々の分散$V(\hat{Y}_i)$が大きく変動するとき、共通の$w$を決めるのが難しくなる。また、プールされる小地域のなかに、他と似ていない小地域があるとき、Jame-Stein推定量は非効率になる。

　もっと単純に$w_i$を決めようという路線もある。２つ紹介しよう。以下、地域の人口$N_i$の直接的な不偏推定を$\hat{N}_i$とする。
　ひとつめ、Drew, Singh, & Choudhry (1982)。もし$\hat{N}_i \geq \delta N_i$だったら$w_i=1$、そうでなかったら$w_i=\hat{N}_i / \delta N_i$とする。$\delta$は適当に決める。
　ふたつめ、Sarndal & Hidiroglou (1989)。もし$\hat{N}_i \geq N_i$だったら$w_i=1$、そうでなかったら$w_i=(\hat{N}_i / N_i)^{h-1}$とする。$h$は適当に決めていいんだけど、著者らは$h=2$を推奨。
　なお、$\delta=1, h=2$とするとこの2つは等しくなる。
　例として、人口$N$から$n$を単純無作為抽出する場面を考えよう。$\hat{N}_i = N(n_i/n)$である。
　もし、$n_i$がその期待値$E(n_i)=n(N_i/N)$と同じくらい大きければ、($\delta=1$とすれば)どちらのやりかたでも$w_i=1$となる。つまり、他の小地域から値を借りてくることができなくなる。$E(n_i)$がいくら小さくてもそうなるわけである。
　もし、$\hat{N}_i \lt N_i$だったら、($h=2$として)どちらのやりかたでも$w_i$は$n_i$の現象と共に減少する。合成推定量の重みがどんどん増すことになる。
　もうひとつ欠点を挙げると、こういう方法は小地域以内の変動に対する小地域間の変動の相対的サイズを無視しているわけで、小地域間の等質性がどうであれ、重みは同じになる。
　[うぐぐ... 途中から混乱しちゃって写経状態になってしまった。あとでゆっくり考えよう]

　Holt, Smith & Tomberlin(1979)はこういうのを考えた。セル$(i,g)$の単位$l$について、
　$y_{igl} = \mu_g + e_{igl}$
とする。大地域$g$の固定効果$\mu_g$のBLUP(最良線形不偏推定量)は$\hat{\mu}_g = \bar{y}_{.g}$である。ここから$Y_i$のBLUPは下式となる：
　$\hat{Y}_i^B = \sum_g \hat{Y}_{ig}^C$
ただし$\hat{Y}_{ig}^C$は、直接推定量$\hat{Y}_{ig} = N_{ig}\bar{y}_{ig}$に重み$w_{ig} = n_{ig}/N_{ig}$、合成推定量$\hat{Y}_{ig}^S=N_{ig}\bar{y}_{ig}$に重み$1-w_{ig}$を与えて足し上げた推定量である。この方法もまた小地域間変動を無視しているわけだけど、モデルに小地域のランダム効果をいれることで改善できる。これが次節以降の手法のもとになっている。[なるほど...]

4. 小地域モデル
　ここからは小地域のランダム効果をいれたモデル。大きく２つにわかれる。

　ひとつめ。小地域の補足データ$\mathbf{x}_i = (x_{i1}, \ldots, x_{ip})^t$をいれる。関心あるパラメータを$\theta_i$として
　$\theta_i = \mathbf{x}_i^t \beta + v_i z_i$
$z_i$は既知の正定数。$v_i$は$E(v_i)=0, V(v_1)=\sigma^2_v$、なんなら正規性を仮定してもよい。
　ここで、直接推定量$\hat{\theta}$について
　$\hat{\theta} = \theta_i + e_i$
かつ$E(e_i | \theta_i) = 0$(つまりデザイン不偏)、$V(e_i | \theta_i) = \psi_i$という風に仮定し、$\psi_i$は既知とみることが多い。実際には$\theta_i$が非線形な関数になっていたりして、不偏という仮定には無理があったりするんだけど。まあとにかく、
　$\hat{\theta}_i = \mathbf{x}_i^t \beta + v_i z_i + e_i$
ということになる。線形混合モデルの特殊ケースである。
　[あとで出てくるけど、こういうのをFay-Herriotタイプのモデルという由]

　ふたつめ。測定単位の補足データ$\mathbf{x}_{ij} = (x_{ij1}, \ldots, x_{ijp})^t$をいれる。
　$y_{ij} = \mathbf{x}_{ij}^t \beta + v_i + e_{ij}$
ここで$e_{ij} = \tilde{e}_{ij} k_{ij}$, $k_{ij}$は既知の定数で、$\tilde{e}_{ij}$は平均0, 分散$\sigma^2$、なんなら$v_i$と$\tilde{e}_{ij}$に正規性を仮定してもよい。こういうのをnested error 回帰モデルという。関心があるのは合計$Y_i$、ないし平均$\bar{Y}_i = Y_i / N_i$である。
　このモデルは選択バイアスを考慮してないことに注意。したがって標本が単純無作為抽出でないとそのままでは使えない。
　このモデルはいろんなところで使われていて[...中略...]、最後に触れるけど、拡張の提案も沢山ある。

5.1 EBLUP(分散成分)アプローチ
　BLUPというのは、固定パラメータのBLUEみたいなもので、正規性を仮定することなく、線形不偏推定量のなかでMSE最小の奴を探すアプローチ。50年代に遡る。

　さっきの$\theta_i$のほうのモデルでいうと、そのBLUPは...[中略] そのMSEは... [中略]。ところがこれらは分散成分$\sigma^2_v$を既知ととっている。実際にはふつう未知なわけで、これをまずMLとかREMLとかで推定する、というのがEBLUP。[... 中略... ]
　では、nested error 回帰モデルのほうはどうなるか ... [...中略...]

5.2 EBアプローチ
　まずデータの周辺分布からモデルのパラメータを推定し、それを使って関心あるパラメータの事後分布を得る、というアプローチ。
　[...時間が無くなってきたのでメモ省略...]

5.3 HBアプローチ
　モデルのパラメータの事前分布を決め、関心あるパラメータの事後分布を得る、というアプローチ。[... 心底面倒くさいので読み飛ばした...]

　Datta & Ghosh (1991) はHB, EB, EBLUPを比較している。[...中略。要するに、点推定についていえば大差ない]

6. 事例
[略]

7.1 モデル診断
[略。適合度の話ではなくて、q-qプロット描きましょう的な話]

7.2 制約付き推定
[小地域の推定値の合計がなにかの既知の値にぴったり合わないと困る、というような場面の話。あるある、そういうの。時間がないので飛ばしたが、結構切実な話だ。いずれ読み直そう]

7.3 拡張
Fay-Herriotタイプのモデル：

nested error回帰モデル：

8. 結論
　最後にひとこと。小地域の間接推定量を使うときには十分に気をつけるように。それはもともと、直接推定量の代替に過ぎないんだからね。

　やれやれ、やっと終わった... 疲れた...

論文：データ解析(2015-) - 読了：Ghosh & Rao (1994) 小地域推定レビュー

2017年11月14日 (火)

Molina, I., Marhuenda, Y. (2015) sae: An R Package for Small Area Estimation. The R Journal, 7(1).
　小地域推定のためのRパッケージsaeの解説。仕事の都合で読んだ。いっちゃなんだけど、面白くも何ともない話題だ... (すいません)

　有限母集団$U$が$D$個の地域に相互排他的かつ網羅的に分割されている。サイズを$N_1, \ldots, N_D$とする。[ここで身構えたけど、有限母集団であることはこの後の話には顔を出さない。やれやれだぜ]
　地域$d$における個人$j$の、ある関心ある変数の測定値の$Y_{dj}$とする。その地域の全員について$\mathbf{y}_d = (Y_{d1}, \ldots, Y_{dN_d})^t$とベクトルで書いちゃうこともある。目標パラメータを$\delta_d = h(\mathbf{y}_d)$とする(たとえば平均とか)。
　各地域についてサイズ$n_d$の下位標本$s_d$が手に入っている。残りの部分を$r_d$とする。

　saeパッケージが提供する小地域推定手法は以下の通り。

1. Fay-Herriotモデルに基づくEBLUP。
　EPLUPってのはあれね、経験最良線形不偏予測量のことね。
　Fay-Herriotモデルってのは79年にアメリカの所得推定のために提案されたモデル。$\delta_d$の直接推定量を$\hat{\delta}_d^{DIR}$とすると、それは不偏であって、
　$\hat{\delta}_d^{DIR} = \delta_d + e_d, \ \ e_d \sim N(0, \psi_d)$
と書ける(これを抽出モデルという)。$\psi_d$は既知とする。いっぽう、地域レベルの補足変数ベクトルを$\mathbf{x}_d$として
　$\delta_d = \mathbf{x}_d^t \mathbf{\beta} + u_d, \ \ u_d \sim N(0, A)$
とする(これをリンキング・モデルという)。$\beta$は全地域共通。$A$は当面既知とする。
　これを併せると
　$\hat{\delta}_d^{DIR} = \mathbf{x}_d^t \mathbf{\beta} + u_d + e_d$
という混合モデルになる。これをFHモデルという。[←意外にあたりまえな話でびっくり...これにわざわざ人名をつけるかね...]

　さて、このモデル$\delta_d$のBLUP
　$\tilde{\delta}_d^{BLUP} = \mathbf{x}_d^t \tilde{\beta}(A) + \tilde{u}_d(A)$
はすでに75年に得られていて、それはこれこれこうである[メモ省略]。問題はそれが$A$の関数だということ。じゃあ$A$はどうすんだと。しょうがない、$\hat{A}$をデータから推定しましょう、というのが、経験BLUP, 略してEBLUPである。
　EBLUPは結局こういう形になる。
　$\hat{\delta}_d^{EBLUP} = \hat{\gamma} \hat{\delta}_d^{DIR} + (1-\hat{\gamma}_d) \mathbf{x}_d^t \hat{\beta}$
ただし$\hat{\beta} = \tilde{\beta}(\hat{A})$。さて、謎の$\hat{\gamma}$というのが出て、よく見ると2つの項に重みをつけているんだけど、これは[...メモ省略...] $\hat{\gamma} = \hat{A}/(\hat{A} + \psi_d)$である。つまりこういうことだ。直接推定量$\hat{\delta}_d^{DIR}$が十分あてになる場合(=$\psi_d$が小さい場合)、EBLUPは直接推定量に近づく。そうでないとき、EBLUPは回帰による推定量に近づく。なるほど、うまいことできてんね。

　$\hat{A}$の一致推定量はいろいろ提案されていて...[メモ略]
　モデル比較にはAIC, BICが使えて...[メモ略]

　というわけで、saeパッケージはEBLUPを求める関数 eblupFH() とそのMSEを解析的に求める関数 mseFH()をご提供しております。ただし、$\psi_d$は教えてやらないといけないので、事前になんらか推定しておくこと。surveyパッケージでもなんでも好きなパッケージでやるがよろしい。saeにもdirect()という関数がある。
　計算例...[略]
　[モデルとしてはごくふつうの階層回帰モデルなので、地域数が少なければMplusでどうにかなっちゃいそう。どのくらいスケールするのかが知りたいところだ...]

2. 空間Fay-Herriotモデルに基づくEBLUP。
　せっかく空間の話してんのに、なぜただの階層モデルなんだよ、と御不満のみなさん。お待たせしました。地域効果$\mathbf{u} = (u_1, \ldots, u_D)^t$にSAR(1)をいれます。イエー。そう来なくっちゃー。
　$\mathbf{u} = \rho_1 \mathbf{Wu} + \mathbf{\epsilon}$
　$\mathbf{\epsilon} \sim N(\mathbf{0}_d, \sigma_1^2\mathbf{I}_D)$
　$\mathbf{0}$は0の列ベクトル、$\mathbf{I}$は単位行列、添え字はサイズ。隣接行列$\mathbf{W}$は対角が0で隣接に1がはいっている行列を行ごとに標準化したものだと思いねえ。
　このモデルのEBLUPも、その解析的なMSEも得られている。saeパッケージではeblupSFH(), mseSFH()をご提供している。ブートストラップMSEとして、パラメトリック版 pbmseSFH(), ノンパラ版 npbmseSFH()もご用意している。[どっち使えばいいんだろう？ 何も書いてない]
　計算例... [パス]

3. 時空間Fay-Herriotモデルに基づくEBLUP。
[いま関心ないのでパス]

4. BHFモデルに基づくEBLUP。
　今度は補足変数が単位レベルで手に入っている場合について考える。地域$d$の単位$j$について$(Y_{dj}, \mathbf{x}_{dj}^t)$が得られているわけだ。
　$Y_{dj} = \mathbf{x}_{dj}^t \mathbf{\beta} + u_d + e_{dj}$
　$u_d \sim N(0, \sigma_u^2)$ (iid)
　$e_{dj} \sim N(0, \sigma_e^2)$ (iid)
とする。Bettesesさんという人たちの88年の論文の著者頭文字をとってBHFモデルと呼ぶ。
　このモデルのEBLUPも得られていて...[メモ省略]。eblupBHF(), pbmseBHF()をご用意している。
　計算例... [パスしちゃったけど、使うときには読んだほうがよさそう]

5. BHFモデルに基づく非線形パラメータの経験最良推定量。
　[いよいよ大ボスの登場だ...]
　$\delta_d = h(\mathbf{y}_d)$が非線形である場合について考える。以下、$\mathbf{y}_d$を並び替え、$(\mathbf{y}_{ds}^t, \mathbf{y}_{dr}^t)^t$としておく[標本を前に持ってくるってことね]。
　もはやBLUPではなく、線形性も不偏性も捨てて、MSEを最小化したい。すなわち、手元の$y_{ds}$に条件づけられた、$y_d$のパラメータ$h(y_d)$の、$y_{dr}$の分布を通じた期待値
　$\tilde{\delta}_d^{B} = E_{y_{dr}}[h(y_d)|y_{ds}]$
が欲しいわけだ。
　こういう場合は、まず$\beta, \theta=(\sigma_u^2, \sigma_e^2)^t$の一致推定量$\hat{\beta}, \hat{\theta}$をMLなりREMLなりで推定する。次にモンデカルロ法で期待値を近似し、$\delta_d$を推定する。これを経験最良(EB)推定量と呼ぶ。
　具体的な手順は次の通り。まず$\hat{\beta}, \hat{\theta}$を推定する。次に、$\hat{\beta}, \hat{\theta}$の下での$h(y_{dr}|y_{ds})$を求める(BHFモデルの場合で言えば正規分布である)。で、そこからベクトル$y_{dr}$をたくさん生成する。それぞれの頭に手元の$y_{ds}$をくっつけて、たくさんの架空のセンサス・ベクトル$y_d$を得る。それぞれの$y_d$から$h(y_d)$を得る。これを平均し、$\hat{\delta}_d^{EB}$とする。

　ここで大変なのは、長ーいベクトル$y_{dr}$をたくさん生成するというところ。実際にはもうちょっと楽な方法がある。BHFモデルより、$j \in r_d$について
　$Y_{dj} = \mathbf{x}_{dj}^t \hat{\beta} + \hat{u} + v_d + \epsilon_{dj}$
　$v_d \sim N(0, \hat{\sigma}_u^2(1-\hat{\gamma}))$
　$\epsilon_{dj} \sim N(0, \hat{\sigma}_e^2)$
なので、この単変量を生成すれば良い。

　さらにこういう手もある。たとえば、いま各地域の貧困生起率$\delta_d$に関心があるとしよう。$O_{dj}$を年収、貧困線を$z$として
　$\delta_d = \frac{1}{N_d} \sum^{N_d} I(O_{dj} \lt z)$
と定義しよう。$O_{dj}$はどうみても歪んでいるのでBFHモデルは仮定できないが、たとえば$Y_{dj} = \log(O_{dj}+c)$というようにして正規近似できるかもしれない。$O_{dj}=T^{-1}(Y_{dj})$とすれば
　$\delta_d = \frac{1}{N_d} \sum^{N_d} I(T^{-1}(Y_{dj}) \lt z)$
というわけで、一般化していえば$O_{dj}$のBox-Cox変換で正規性を得ることでどうにかなるかもしれない。

　saeパッケージではebBHF(), pbmseebBHF()を御用意しており、Box-Cox変換とかも指定できます。...[などなど。途中で疲れて読み飛ばした]
　計算例[パス]。

　最後に、小地域推定のためのそのほかのRパッケージをご紹介。

[ちなみにこのsaeパッケージは2015年で止まっている。CRAN Task ViewのOfficial Statistics & Survey Methodologyには小地域推定というセクションがあって、載っているのはrsae, hbsae, josae, そしてnlmeとlme4。このパッケージ載ってないじゃん....]

論文：データ解析(2015-) - 読了：Molina & Marhuenda (2015) 小地域推定のためのRパッケージsae

2017年11月13日 (月)

　仕事でデータの分析していて、この場面ではこれこれこういうわけでこのような方法で推定しております、これは当該分野の常識でございます、なあんていかにも専門家づらでにこやかに語っているんだけど、心のなかでは、(ああ、こういうのってちゃんとした名前があるんだろうな... 俺がきちんと勉強してないだけで...) と不安を抱えていることが結構ある。
　これは辛い。かなり辛い。人生間違えたな、生まれ変わってやり直したいな、と思う瞬間が毎日の生活の中に沢山あるけれど、これもそのひとつである。

　たとえば、なにかについて推定していて、ああ、これって要するにスタイン推定量って奴じゃないのかしらん、って思うことがたまにある。でもそういうの、誰にもきちんと教わったことがない。辛い。
　あんまり辛いので、月曜朝に早起きしてノートをとった。本来、仕事が溜まっていてそれどころじゃないのであって、こういう現実逃避をしているから人生ぱっとしないのだともいえる。まあとにかく、以下は自分用の覚え書きであります。
　
　なんでもいいんだけど、いま$N$個の母集団特性$\mu_1, \ldots, \mu_N$があって、それぞれの$\mu_i$について$z_i \sim N(\mu_i, \sigma^2_0)$が独立に観測されるとする。$z_i$は1個の観測値でもいいし標本平均でもいい。また分散は$N$個の特性を通じて等しければなんでもよい。
　個別の$\mu_i$について推定する場面を考える。その推定量$\hat{\mu}_i$として何を使うのがいいか。
　もちろん$z_i$そのものであろう。$z_i$は最小分散不偏推定量であり、最尤推定量であり、MSEを損失としたときの許容的推定量である(=「$\mu_i$の値がなんであれ$z_i$よりも誤差の二乗の期待値が小さい推定量」は存在しない)。

　ところが。
　今度は$\mathbf{\mu}=(\mu_1, \ldots, \mu_N)^t$について同時推定する場面を考える。いま、$N$個の母集団特性を通じたMSEの和
　$\sum^N E[(\hat{\mu}_i-\mu_i)^2]$
が小さいと嬉しい、としよう。これを最小化する推定量はなにか。
　$N$が3以上の時、驚いたことに、それは$\mathbf{z} = (z_1, \ldots, z_N)^t$ではない。実は、$\mathbf{\mu}$がなんであれ、MSEの和が$\mathbf{z}$より小さくなるような推定量が存在する。それが有名なJames-Stein推定量
　$\displaystyle \mathbf{\delta} = \left( 1 - \frac{N-2}{||z||^2}\right) \mathbf{z}$
である。
　母集団特性$\mu_1, \ldots, \mu_N$の間にはなんの関係もない。従って、たとえば$\mu_1$を推定する際に役立つのは$z_1$だけであるはずであって、$z_2, \ldots, z_N$を使うのはおかしい。なぜこんなことが起きるのか？
　1950年代の統計学を揺るがせた、スタイン・パラドクスの登場である。

　この奇妙な現象を説明する方法はいくつかある。そのひとつが、James-Stein推定量を経験ベイズ推定量として捉える説明である。

　母集団特性を確率変数と見なし、$\mu_i \sim N(0, A)$としよう。
　話を簡単にするために、$z_i$の分散は当面 $1$ としておく。すなわち $z_i|\mu_i \sim N(\mu_i, 1)$。
　ベイズの定理より、$z_i$の下での$\mu_i$の事後分布は
　$\mu_i | z_i \sim N(Bz_i, B), \ \ B = A/(A+1)$
であることが示せる。
　複数の母集団特性について一気に書こう。$I$を単位行列として
　$\mathbf{\mu} \sim N_N(0, AI)$
　$\mathbf{z} | \mu \sim N_N(\mathbf{\mu}, I)$
　$\mathbf{\mu} | \mathbf{z} \sim N(B\mathbf{z}, BI), \ \ B = A/(A+1)$

　$\mathbf{\mu}$の推定誤差を最小二乗誤差で表すことにしよう。
　$L(\mathbf{\mu}, \hat{\mathbf{\mu}}) = || \hat{\mathbf{\mu}} - \mathbf{\mu}||^2 = \sum^N(\hat{\mu}_i - \mu_i)^2$
　リスク関数を、所与の$\mathbf{\mu}$の下での推定誤差の期待値としよう。
　$R(\mathbf{\mu}) = E_\mu[L(\mathbf{\mu}, \hat{\mathbf{\mu}}) ]$
$\mathbf{\mu}$は固定で$\mathbf{z}$が動くということをはっきりさせるために$E_\mu$と書いている。

　さて、$\mathbf{\mu}$の最尤推定量は$\mathbf{z}$そのものである。
　$\hat{\mathbf{\mu}}^{(MLE)} = \mathbf{z}$
そのリスクは、さきほど分散を$1$にしておいたので
　$R^{(MLE)}(\mathbf{\mu}) = N$
　いっぽう、$\mu_i \sim N(0, A)$というベイズ的信念の下では、最小二乗誤差の期待値を最小化する推定量は事後分布の平均である。
　$\hat{\mathbf{\mu}}^{(Bayes)} = B\mathbf{z} = (1-\frac{1}{A+1}) \mathbf{z}$
そのリスクは、$\mathbf{\mu}$を固定した状態では
　$R^{(Bayes)}(\mathbf{\mu}) = (1-B)^2||\mathbf{\mu}||^2+NB^2$
$A$を固定して$\mathbf{\mu}$を動かした全域的なリスクは
　$R_A^{(Bayes)} = E_A[R^{(Bayes)}(\mathbf{\mu})] = N \frac{A}{A+1}$
$R^{(MLE)}(\mathbf{\mu}) = N$と比べると、$\frac{A}{A+1}$倍に減っているわけである。

　ところが問題は、$A$が未知だという点である。そこで、$A$を$\mathbf{z}$から推測することを考える。いやぁ、大人はずるいなあ。
　$\mathbf{z} | \mu \sim N_N(\mathbf{\mu}, I)$の周辺分布をとると
　$\mathbf{z} \sim N_n(0, (A+I)/I)$
となる。ということは、$\mathbf{z}$の二乗和$S=||\mathbf{z}||^2$は、自由度$N$のカイ二乗分布を$A+1$倍した分布に従い
　$S \sim (A+1) \chi^2_N$
ここから下式が示せる：
　$E[\frac{N-2}{S}] = \frac{1}{A+1}$
やれやれ、というわけで、これを$\hat{\mathbf{\mu}}^{(Bayes)}$に代入して得られるのが、James-Stein推定量
　$\hat{\mathbf{\mu}}^{(JS)} = \left(1-\frac{N-2}{S}\right) \mathbf{z}$
である。
　その全域的なリスクは下式となる。
　$R_A^{(JS)} = N \frac{A}{A+1} + \frac{2}{A+1}$
$R_A^{(Bayes)}$よりもちょっと大きくなるけど、たとえば$N=10, A=1$のときにはたった2割増しである。

　$N \geq 4$の場合について、もっと一般的に書き直しておこう。
　$\mu_i \sim N(M, A)$ (iid)
　$z_i|\mu_i \sim N(\mu_i, \sigma^2_0)$ (iid)
として、
　$z_i \sim N(M, A+\sigma^2_0$ (iid)
　$\mu_i | z_i \sim N(M+B(z_i - M), B\sigma^2_0), \ \ B = \frac{A}{A+\sigma^2_0}$
となり、
　$\hat{\mu}_i^{(Bayes)} = M + B(z_i - M)$
だが$M, B$がわからない。そこでJames-Stein推定量の登場である。$\bar{z} = \sum z_i/N, S=\sum(z_i - \bar{z})^2$として
　$\hat{\mu}_i^{(JS)} = \bar{z} + \left(1-\frac{(N-3)\sigma^2_0}{S}\right) (z_i - \bar{z})$

　このタイプの推定量が役に立つのはどんな場面か。
　まずいえるのは、たくさんの同種類の量を同時に推定する場面であること。さらに、$X_i$の分散が大きくて困っている場面であること(でなければ、不偏性をなくしてまで$\hat{\mu}_i$を改善しようとは思わない)。
　さらに付け加えると、なんらかの先験情報が存在すること。上記の説明だと、事前に$\mu_i \sim N(M, A)$という信念があるわけだけど、こういう風に、$\mu_i$がなにかに近い、という風に考えることができる場合に役に立つ。だから、形式的にいえば$X_1, X_2, \ldots, X_N$がお互いに全く無関係な事柄についての値であってもJames-Stein推定量は使えるんだけど、役に立つといえるのは、やはり、なんらかの意味で同じ種類の量についての同時推定の場面である。

　以上、主に次の2つの資料の、それぞれ最初の数ページだけを読んで取ったメモである。続きはまたいずれ。
　篠崎信雄(1991) Steinタイプの縮小推定量とその応用. 応用統計学, 20(2), 59-76.
　Efron, B. (2012) Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction. Cambridge University Press.

雑記：データ解析 - 覚え書き：スタイン推定量とはなんぞや

« 2017年10月 | メイン | 2017年12月 »

rebuilt: 2020年11月16日 22:39
validate this page