elsur.jpn.org >

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$を決めようという路線もある。2つ紹介しよう。以下、地域の人口$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. 小地域モデル
 ここからは小地域のランダム効果をいれたモデル。大きく2つにわかれる。

 ひとつめ。小地域の補足データ$\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回帰モデル:

  • 多変量化するという話。[あああ... そうか... これは結構切実だ... Fuller & Harter (1987)]
  • 小地域のなかで2段階抽出しているときのモデル。
  • $\mathbf{x}_{ij} = 1$のときのモデル。[? それってどう特別なのだろうか。Kleffe &: Rao(1992)]
  • 反応が二値の場合のモデル。ロジスティック回帰モデルを使う。
  • ポワソン回帰を使うモデル。死亡力とか疾病割合とかに使う。
  • 二変量のモデル。2つのガンの同時死亡力について使う。

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

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

論文:データ解析 - 読了: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パッケージをご紹介。

  • rsaeパッケージ。ロバスト・フィッティングの関数を提供。[2014年で止まっている]
  • JoSaeパッケージ。分散行列がブロック対角な単位レベルモデルのみ。[2015年で止まっている]
  • hbsaeパッケージ。REMLと階層ベイズ推定を提供。[2012年で止まっている]
  • mmeパッケージ。多項線形混合モデル。[2014年で止まっている]
  • saeryパッケージ。時間効果のはいった地域レベルのEBLUP。[2014年で止まっている]
  • sae2パッケージ。時系列地域レベルモデル。[2015年で止まっている]

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

論文:データ解析 - 読了: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月27日 (金)

Bivand, R.S., Gomez-Rubio, V., Rue, H. (2015) Spatial data analysis with R-INLA with some extensions. Journal of Statistical Software, 63(20).
 MCMCではない謎の方法(integrated nested Laplace approximation; 積分段階的ラプラス近似)で空間統計モデルを推定しちゃう謎のRパッケージ R-INLA の解説書。できればそういう難しい話とは関わり合いを持ちたくなかったんだけど。人生とはままならないものだ。

 モデル。
 観察$\mathbf{y} = \{y_i\}_{i=1}^{n}$が指数型分布族に従っていて、$y_i$の平均$\mu$が線形予測子$\eta_i$となんらかのリンク関数でつながっているとする。
 $\eta_i = \alpha + \sum_j^{n_f} f^{(j)}(u_{ij}) + \sum_k^{n_\beta} \beta_k z_{ki} + \epsilon_i$
 $f^{(j)}$は共変量$\mathbf{u}$のランダム効果を表す$n_f$個の関数。$\beta_k$は共変量$\mathbf{z}$の関数。

 潜在効果のベクトルを$\mathbf{x} = \{\{\eta_i\}, \alpha, \{\beta_k\}, \ldots\}$と書く。
 [←この論文、よくわからなくて困っているんだけど、このメモをとった翌日になってじっくり読み直してみて、混乱の根底はここだと気がついた。ここに$\eta_i$と$\alpha, \beta_k$の両方がはいる理由がよく分からない。$\mathbf{x}$は地点を要素にもつ、共変量の効果抜きの効果なのではないかと思うのだけれど...$\mathbf{x} = \{\epsilon_i\}$というのなら納得するんだけど...]

 $\mathbf{x}$について、平均0, 共分散行列$Q(\theta_1)$の正規分布を仮定する。さらに、「$\mathbf{x}$はガウシアン・マルコフ確率場(GRMF)であると仮定する。これは$Q(\theta_1)$が数多くのマルコフ特性を満たすことを意味する。」[←素人を殺しに来たぞ... なにをいっているのかさっぱりだが、諦めずに読み進めることに]
 $y_i$の分布は潜在効果$\mathbf{x}$に依存するが、他のハイパーパラメータ$\theta_2$にも依存しうる。[←どういうこと? $y_i$の方程式の右側に$\eta_i$以外の項があるってこと?] 以下では$\theta=(\theta_1, \theta_2)$としよう。
 観察$y_i$は所与の$x_i, \theta$の下で互いに独立である。なぜなら$\mathbf{x}$はGRMFだから。[←わ、わからん...いったい何を云っているのだ...]

 $\mathbf{x}$と$\theta$の事後分布はこうなる。
 $\pi(\mathbf{x}, \theta | \mathbf{y})$
 $\propto \pi(\theta) \pi(\mathbf{x} | \theta) \prod_{i \in I} \pi(y_i | x_i, \theta)$
 $\propto \pi(\theta) |\mathbf{Q}(\theta)|^{\frac{1}{2}} \exp (-\frac{1}{2} \mathbf{x}^T \mathbf{Q}(\theta) \mathbf{x} + \sum_{i \in I} \log \pi (y_i | x_i, \theta) )$
 ただし、$I$は観察データのインデクス($1, \ldots, n$)。$\mathbf{Q}(\theta)$は$\theta$の精度行列。最後に出てくる$ \log \pi (y_i | x_i, \theta)$は$y_i$の対数尤度。
 [えーと、数式2行目から3行目に行くところで完全にギブアップなんだけど、信じましょう]

 なお、INLAは指数型分布族だけではなく、たとえば混合分布であっても対応できる。$\mathbf{x}$について上記と異なる定式化も可能。事前分布$\pi(\theta)$についてもとても柔軟である。

 さて、$\mathbf{x}$と$\theta$の周辺分布は
 $\pi(x_i|\mathbf{y}) = \int \pi(x_i|\theta, \mathbf{y}) \pi(\theta|\mathbf{y}) d\theta$
 $\pi(\theta_j|\mathbf{y}) = \int \pi(\theta|\mathbf{y}) d\theta_{-j}$
となるのだが、この右辺の項をだね、うまいこと近似する方法があるのだ。おまえのようなアホに説明することは到底不可能だが、ラプラス近似という方法をさらにごにょごにょしちゃうことによってなにかがどうにかなっちゃうのだ。[←こんな風には書いてないが、私としてはこんな感じ。正直、ここの数段落は全く理解できない]

 というわけで、この方法を実装したソフトINLA、そしてそれをRで使うパッケージR-INLAを作ったので、使っても苦しゅうないぞ。

 パッケージの具体的な使い方。最初はラティスデータの分析の話から。
 まず次の点に注意。無相関なランダム効果を入れる場合、ランダム効果は
 $u_i \sim N(0, \tau_u)$
 と定義され、$\tau_u$じゃなくて$\log(\tau_u)$に事前分布が与えられる(デフォルトでは対数ガンマ分布)。

 空間的相関をモデル化するためには隣接エリアを定義しないといけない。空間的自己相関は、平均0, 共分散が
 $\Sigma = (1/\tau) Q^{-1}$
の正規分布でモデル化される。$Q$は、エリア$i$と$j$が隣接じゃない時には$Q_{ij}=0$となる。$Q$はふつうすごくスパースになる。[←あー!それでマルコフ場っていうのか!おとなりさんとしか相関しないというモデルなのね。やっと意味が分かった...そりゃあ計算も普通のクリギングより速くなるかもな]

 こういうときの$Q$の決め方は、いろいろご用意してある。いくつかご紹介しよう。
 ひとつめ、intrinsic CARモデル。$Q_{ii}$を隣接領域数$n_i$とし、 $Q_{ij}(i \neq j)$は隣接していれば$-1$, でなければ$0$。これは結局次のモデルになる。空間的ランダム効果を$v_i$、精度を$\tau_v$として, $i \neq j$のとき
 $\displaystyle v_i | v_j, \tau_v \sim N\left( \frac{1}{n_i} \sum_{i \sim j} v_j, \frac{1}{\tau_v n_i} \right) $
 $\log(\tau_v)$には事前分布として対数正規分布を与える。
 その2、proper CARモデル。上のモデルをちょっと変えて分布を正則にする。$d$は正の値。
 $\displaystyle v_i | v_j, \tau_v \sim N\left( \frac{1}{n_i+d} \sum_{i \sim j} v_j, \frac{1}{\tau_v (n_i+d)}\right)$
 その3、もっと一般的なアプローチ。$I$を単位行列、$\rho$を空間的自己相関パラメータ、$C$を隣接行列, $\lambda_{max}$を$C$の最大固有値として
 $\displaystyle Q = (I - \frac{\rho}{\lambda_{max}} C)$
で、$\log(\rho/(1-\rho))$に事前分布として正規分布を与える。
 [こうやって列挙されてもねえ... それぞれのアプローチのモチベーションとか長短とかを教えてくださいよ...]

 ここからは、具体的な分析例がたくさん。
 大事そうなことがいろいろ書いてあるんだけど、力尽きたのでいったん読了としておく。ま、雰囲気がわかったのでよしとしよう。あれだね、バリオグラム関数はモデル化せず、隣との相関だけをモデル化するわけね。

論文:データ解析 - 読了:Bivand, Gomez-Rubio, Rue (2015) R-INLAパッケージで楽しい楽しい空間統計モデリング

2017年10月26日 (木)

Jing, L., De Oliveira, V. (2015) geoCount: An R package for the analysis of geostatistical count data. Journal of Statistical Software. 63(11).
 カウントデータに特化した空間モデリングパッケージ geoCountについての紹介。実戦投入前の儀式として読んだ。

 カウントデータの空間モデリングにはふつう一般化線形空間モデル(GLSM)が使われるが、実際にはいろいろ大変である。(1)なにより計算が大変。メモリが足らん。(2)MCMCがめっちゃ遅い。(3)MCMCが収束しない。
 そこで新しいパッケージをお届けしよう。このパッケージは、C++で書いてあり速い。並列計算に対応しておる。そして効率的な新アルゴリズム。それではご紹介しましょう! geoCountパッケージです!

 モデル。
 $\{S(\mathbf{x}): x \in A\}$を正規確率場とする。$\mathbf{S} = (S(\mathbf{x}_i), \ldots, S(\mathbf{x}_n))^T$とする。
  $Y_i | S(\mathbf{x}_i) \sim p(\cdot | \mu_i), \ \ i=1, \ldots, n$
  $\mu_i = t_i g^{-1} (S(\mathbf{x}_i))$
  $\mathbf{S} \sim N_n (D\beta, \Sigma)$
  $(\beta, \theta) \sim \pi(\beta, \theta)$

 説明しよう。
 一本目。$\mu_i = E(Y_i | S(\mathbf{x}_i))$である。$p(\cdot | \mu_i)$については後述。
 二本目。$g(\cdot)$はリンク関数である。
 三本目。$D$はフルランク計画行列で、なかに共変量がはいっておる。$\Sigma = (\sigma_{ij})$は正定な共分散行列で、要素$\sigma_{ij}=\sigma^2 \rho(u_{ij})$はユークリッド距離$u_{ij}=||\mathbf{x}_i - \mathbf{x}_j||$で決まる(つまり等方性がある)。
 四本目。$\pi(\beta, \theta)$はパラメータの事前分布。

 $\rho(u)$としては、球面ファミリー、Maternファミリー、power exponentialファミリーをご用意しておる。
 $p(\cdot | \mu_i)$とリンク関数$g(\cdot)$については、良く用いられる次の2つをご用意。
 その1、ポワソン対数正規モデル:
  $Y_i | S(\mathbf{x}_i) \sim Poisson(\mu_i)$
  $\mu_i = t_i \exp(S(\mathbf{x}_i))$
 その2、二項ロジット正規モデル:
  $Y_i | S(\mathbf{x}_i) \sim Binomial(t_i, \mu_i/t_i)$
  $\mu_i = t_i \exp(S(\mathbf{x}_i)) / (1+\exp(S(\mathbf{x}_i)))$

 事前分布は次のようにご用意しております...[めんどくさいので略]

 プログラミング上の工夫と並列処理...[パス]
 アルゴリズムの工夫...[パス]

 プログラム例... [パス。実際に使うときに読めばいいや]

 他のパッケージとの比較。

  • 同様のモデルはgeoRglmパッケージでも組めるけど、MCMCの際のパラメータ化のやり方が違っており...[よくわからんので中略... 何だか知らんが、うちらのほうがイケているよとのこと]。ただし、計算時間はこっちのほうがかかるかも。[←よく考えてみたらオイオイって話ですね。C++かつ効率的アルゴリズムで速いんちゃうんかと]
  • INLAパッケージと比べると...そもそもINLAは潜在ガウスモデルを共分散関数として特徴づけるのではなくて、ガウス・マルコフ確率場(GMRF)として特徴づける。ベイズ推論の際、MCMCは周辺分布を確率的に近似しようとするが、INLAはラプラス近似とnumerical quadrature rules[←なにそれ]を用いて決定論的に近似する。
     GMRFを使う利点は計算が速いこと。欠点は、定常性とか等方性といった特徴を持つ共分散関数を持つGRMFを構築するのが難しいこと。しかしこれには部分的解決策が提案されている。Matern族の一種である、ある族に属する共分散関数を持つガウス確率場は、ガウシアン・ホワイトノイズを持つある種のstochastic partial differential equations(SPDE)の解であることが示されているのだ。これらの解はある種のGMRFによって与えられる、ランダム係数を持つある種の基底関数へと拡張できることが示されている。その結果、Matern族の共分散関数を使ってモデリングつつ、それと対応する確率場をGMRFで表現して計算をすることができる。これがINLAで採用されているアプローチである。[←ああなるほど...って、何言ってんだかさっぱり分かんないよこの野郎]
     geoCountと比べてみると、計算は確かに速い。結果は似てたり違っていたり。[←詳細略]

 限界。

  • データがでかいと遅くなる。たまに収束しないこともある。どちらの場合もINLAがおすすめ。
  • ポワソン分布のときは対数リンク、二項の場合はロジットリンクしか選べない。実はidentityリンクのほうがいいという話もある[←へー]。
  • 事前分布についても改善の余地がある。

 云々。

2017/10/27追記: INLAと比較しているくだりについて、メモを追加。

論文:データ解析 - 読了:Jing & De Oliveira (2015) RのgeoCountパッケージ

 空間データのモデリングではバリオグラムというのを描くけれど、たいてい従属変数は量で、誤差に正規性が仮定されている。もし従属変数が二値とかカウントとかだったら、モデルはGLMとして組めるにしても、バリオグラムはいったいどうなるのか。
 これ、ここんところの疑問の一つだったのだが、このたびやっと説明を見つけたのでメモしておく。すいません、純粋に自分のための覚え書きです。

 $S(x)$を平均0, 分散$\sigma^2$の定常ガウス過程とする。観察$Y_i$は互いに独立で平均$\mu_i=g(\alpha+S(x_i))$, 分散$v_i = v(\mu_i)$とする。$g(\cdot)$はリンク関数$h(\cdot)$の逆関数である。

 $Y$の素直なバリオグラムについて考えよう。横軸は$u=||x_i - x_j||$、縦軸は
 $\tau_Y(u)=E[(1/2)(Y_i-Y_j)^2]$
これを展開するとこうなる。
 $\tau_Y(u)$
 $= (1/2)E_S[{g(\alpha+S(x_i)) - g(\alpha+S(x_j))}]$
 $+ E_s[v(g(\alpha+S(x_i)))]$
第2項に$i$だけ出てきているのは、$S(x_i)$の周辺分布が場所を問わず一定だから。よくよくみればこの項は定数なので、以後$\bar{\tau}^2$と書く。これ、$S(\cdot)$の条件付き分布の分散を平均したもので、ガウシアンモデルで言うところのナゲットである。

 さて。第一項に出てくる$g(\alpha+S(x_i))$が、$g(\alpha) + S(x_i) g'(\alpha)$とテイラー近似できるのはわかるね。[←わかんないよ!こっちは文系なんだから!]
 すると結局こうなる。
 $\tau_Y(u) \approx g'(\alpha)^2 \gamma_S(u) + \bar{\tau}^2$
つまり、バリオグラムの縦軸は、潜在ガウス過程$S(\cdot)$だけを取り出したバリオグラム$\gamma_S(u)$に比例し、そこにさきほどの平均ナゲット効果を乗せたものになる。
 というわけで、上記のような通常のバリオグラムをナゲットのサイズ把握のために使うのはOK。しかし、GLMは$S(x)$と$Y$のあいだに非線形的関係を仮定しているわけだから、診断のために使うのはおかしい。[←「診断のため」ってどういうことだろう。バリオグラム・モデルの診断には使えないということだろうか]

 バリオグラム・モデルをどうやって推定するか。
 まずは一般的に、GLMの尤度関数について考えよう。ランダム効果$S$の下での$Y$の条件付き分布を$f_i(y_i; S, \theta)$とすると、パラメータ$\theta$の条件付き尤度は
 $L(\theta | S) | \prod_i f_i(u_i|S, \theta)$
混合モデルの尤度は、$S$の同時分布を$g(S; \phi)$として
 $L(\theta, \phi) = \int_S \prod_i f_i(u_i|S, \theta) g(s|\phi) ds$
さて、こんな風に多重積分が掛かっている尤度が計算できちゃうのは、ひとえに$S_i$が互いに独立だからである。しかし空間データでは$S_i = S(x_i)$は互いに独立でないので、こんなもん計算できっこない。実をいうとこの積分を近似する方法も提案されているんだけど、$S$の分散が大きいときには役に立たない。

 というわけで... モンテカルロ最尤法、階層尤度、GEEなどの方法が提案されておる、という話になるのだが、力尽きたのでこの辺にしておこう。

 以上, Diggle & Ribeiro (2010) "Model-based Geostatistics", Chapter 4-5より。勤務先の経費とはいえ2万円以上した本だったので、元を取らねば、と読んだ次第。
 実をいうと、ポワソン対数線形モデルを使ったクリギングの際に、いい感じのバリオグラム・パラメータを視認で決める方法について知りたかったんだけど、それはよくわかんなかった。先の章に書いてあるのかも。 

雑記:データ解析 - GLM空間モデルのバリオグラムとはなんぞや

<< 読了:Finley, Banerjee, Gelfand (2015) RのspBayesパッケージ
 
validate this page / CSS