この文書は、 Rao & Molina (2015) “Small Area Estimation”, 2nd Ed. (以下RM本)を読みながら、小地域推定の具体的な手続きについて、 Rで実習してみた記録である。
いうまでもなく、私による私のための私のノートであり、すべての誤りは私に帰属いたします。
なお、RM本の読書ノートは理論編として別途まとめてある。
更新履歴:
- 2018/06/25 理論編から分離。Fay-Herriotモデルの説明を追加。 まだまだこれからだ…
- 2018/06/27 Fay-Herriotモデル(milkデータ)を, BayesSAEパッケージ, nlmeパッケージでも試した。
- 2018/06/28 ちょっとだけ更新。空間FHモデルに進む勇気がまだ出ない…
- 2018/06/29 Fay-Herriotモデルのシミュレーションを追加。いろいろと 発見がありました。
- 2018/07/01 空間Fay-Herriotモデルの実習を追加。
- 2018/07/02 空間Fay-Herriotモデルの実習に、BayesSAEパッケージとspdepパッケージを 追加。RM本のささやかな誤りのせいで、えらい目に…
- 2018/07/05 誤字などの修正。
- 2018/07/08 ようやく空間Fay-Herriotモデルの実習を終了。長かった。
- 2018/07/09 Fay-Herriotモデルの実習にrstanを追加。
- 2018/07/10 空間Fay-Herriotモデルの実習にrstanを追加。
- 2018/07/12 誤字など訂正。いったんこれで終了とします。
- 2018/07/18 誤字など訂正。
- 2018/07/19 誤字など訂正。
なにかの事柄について標本調査をやったとき、標本サイズは全体としては十分に大きいんだけど、 標本がなんらかのたくさんのドメイン(たとえば地域とか)に分かれていて、 ドメイン別にみると標本サイズが足りない、 困ったな、ということがある。
標本サイズが足りないんなら、いさぎよく推定をあきらめればいいと 思うんだけど、なかなかそうもいっていられない。 標本サイズが小さいドメインについても、なんとかしてそれらしい推定値を求めたい。 こういう問題を小地域推定という。
小地域推定は主に公的統計の分野で発展してきたんだけど、標本調査を利用する他の分野 (たとえば市場調査)の関係者にとっても、学ぶ価値があるように思う。
というわけで、以下ではRを使い、小地域推定の実習をいたします。
小地域推定の歴史は長く、いろんなアプローチがあるが、それらは以下の3つに分類できる。
本稿ではモデル・ベース間接推定、つまり、小地域モデルによる小地域推定を扱う。
先に述べたように、 地域モデルによる小地域推定では、 (1)データの背後になんらかの小地域モデルを想定し、 (2)そのモデルをデータにあてはめて、モデルのパラメータを推定し、 (3)それを通じて、小地域の特性について推定しようとするわけである。
(1)の小地域モデルにはいろいろあるが、いずれの場合にも、 (2)(3)における推定方法(推定量のタイプ)は、大きく分けて次の3種類である。
以下ではEBLUP推定とHB推定を試す。
EBLUP推定では、モデルのパラメータを次の2段階にわけて推定する。
HB推定では、データとモデルに基づき、小地域の特性がどうなっているかを 確率分布の形で推定する(これを「事後分布」という)。事後分布を求める方法には、 大きく分けて次の2つがある。
以下の実習では次のパッケージを使用する。
### x <- Sys.setlocale("LC_ALL","English")
###
library(dplyr)
library(tidyr)
library(ggplot2)
library(gridExtra) ### grid.arrange()
library(R.utils) ### witTimeout()
library(GGally) ### ggpairs()
library(igraph) ### plot network graph
library(tagcloud) ### smoothPalette()
library(raster) ### required by {rstan}
# library(RColorBrewer)
# library(scales)
# library(coda)
###
### for small area estimation
library(sae)
library(hbsae)
library(BayesSAE)
library(nlme)
library(spdep)
library(rstan)
###
### for rstan
rstan_options(auto_write = TRUE)
options(mc.cores=parallel::detectCores())
sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
Matrix products: default
locale:
[1] LC_COLLATE=Japanese_Japan.932 LC_CTYPE=Japanese_Japan.932
[3] LC_MONETARY=Japanese_Japan.932 LC_NUMERIC=C
[5] LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rstan_2.17.3 StanHeaders_2.17.2 spdep_0.7-7
[4] spData_0.2.9.0 nlme_3.1-137 BayesSAE_1.0-2
[7] lattice_0.20-35 coda_0.19-1 Formula_1.2-3
[10] hbsae_1.0 sae_1.2 lme4_1.1-17
[13] Matrix_1.2-14 MASS_7.3-50 raster_2.6-7
[16] sp_1.3-1 tagcloud_0.6 Rcpp_0.12.17
[19] igraph_1.2.1 GGally_1.4.0 R.utils_2.6.0
[22] R.oo_1.22.0 R.methodsS3_1.7.1 gridExtra_2.3
[25] ggplot2_3.0.0 tidyr_0.8.1 dplyr_0.7.6
loaded via a namespace (and not attached):
[1] deldir_0.1-15 gtools_3.8.1 assertthat_0.2.0
[4] rprojroot_1.3-2 digest_0.6.15 R6_2.2.2
[7] plyr_1.8.4 backports_1.1.2 stats4_3.5.1
[10] evaluate_0.11 pillar_1.3.0 rlang_0.2.1
[13] lazyeval_0.2.1 gdata_2.18.0 rstudioapi_0.7
[16] minqa_1.2.4 nloptr_1.0.4 gmodels_2.18.1
[19] rmarkdown_1.10 splines_3.5.1 stringr_1.3.1
[22] munsell_0.5.0 compiler_3.5.1 pkgconfig_2.0.1
[25] base64enc_0.1-3 htmltools_0.3.6 tidyselect_0.2.4
[28] expm_0.999-2 tibble_1.4.2 reshape_0.8.7
[31] crayon_1.3.4 withr_2.1.2 grid_3.5.1
[34] jsonlite_1.5 arm_1.10-1 gtable_0.2.0
[37] magrittr_1.5 scales_0.5.0 stringi_1.2.3
[40] LearnBayes_2.15.1 bindrcpp_0.2.2 boot_1.3-20
[43] RColorBrewer_1.1-2 tools_3.5.1 glue_1.3.0
[46] purrr_0.2.5 abind_1.4-5 parallel_3.5.1
[49] yaml_2.1.19 inline_0.3.15 colorspace_1.3-2
[52] knitr_1.20 bindr_0.1.1
また、チャート作成のために、あらかじめ以下の汎用的な関数を定義しておく。
sub_PlotTwoEstimate.1 <- function(
agEst1, agMSE1=NULL, agEst2, agMSE2=NULL, asLabel, sVarLabel1, sVarLabel2
){
# purpose: plot scatter plots to compare two type of estimation
# arg:
# agEst1: estimates of type 1
# agMSE1: MSEs of agEst1 (or NULL)
# agEst2: estimates of type 2
# agMSE2: MSE2 of agEst2 (or NULL)
# asLabel: labels
# sVarLabel1: a label of estimation type 1
# sVarLable1: a label of estimation type 2
bPlotMSE <- !is.null(agMSE1)
# traps
stopifnot(length(agEst1) == length(agEst2))
stopifnot(length(agEst1) == length(asLabel))
if (bPlotMSE){
stopifnot(length(agEst1) == length(agMSE1))
stopifnot(length(agEst1) == length(agMSE2))
}
# plot of estimates
dfPlotEst <- data.frame(
gEst1 = agEst1,
gEst2 = agEst2,
sLabel = asLabel
)
g1 <- ggplot(data=dfPlotEst, aes(x = gEst1, y = gEst2, label=sLabel))
g1 <- g1 + geom_text(size=3)
g1 <- g1 + labs(
title = expression("Estimate"),
x = sVarLabel1,
y = sVarLabel2
)
g1 <- g1 + geom_abline(intercept = 0, slope = 1)
g1 <- g1 + theme_bw()
# print(g1)
# plot of MSEs
if (bPlotMSE){
dfPlotMSE <- data.frame(
gMSE1 = agMSE1,
gMSE2 = agMSE2,
sLabel = asLabel
)
g2 <- ggplot(data=dfPlotMSE, aes(x = gMSE1, y = gMSE2, label=sLabel))
g2 <- g2 + geom_text(size=3)
g2 <- g2 + labs(
title = expression("MSE"),
x = sVarLabel1,
y = sVarLabel2
)
g2 <- g2 + geom_abline(intercept = 0, slope = 1)
g2 <- g2 + theme_bw()
}
# plot
if (bPlotMSE){
grid.arrange(g1, g2, ncol = 2)
} else {
plot(g1)
}
}
gg_color_hue <- function(n) {
# purpose: Emulate ggplot2 default color palette
# arg:
# n: number of colors
# return:
# a vector of colors
# notes;
# see https://stackoverflow.com/questions/8197559/emulate-ggplot2-default-color-palette
hues = seq(15, 375, length = n + 1)
hcl(h = hues, l = 65, c = 100)[1:n]
}
以下では、小地域推定に用いるパッケージをご紹介します。
RM本の著者のひとりMolinaさんが開発に関わっているパッケージ。
このパッケージが提供する関数は以下の通り。
saeパッケージで小地域モデルを推定する場合、推定量はEBLUPである。分散成分の 推定量はいくつかの選択肢から選ぶことができる。
小地域推定のパッケージとしては古株。 開発者の方のプレゼン資料によれば、この方はオランダ統計局にお勤めで、このパッケージはオランダの労働力調査で実際に使われているんだそうです。そういわれるとぐっと信頼感が増しますね。
このパッケージが提供する関数は以下の通り。
小地域モデルを推定する場合、推定量としてEBLUPとHBが用意されている。EBLUP推定の場合、 分散成分の推定量はいくつかの選択肢から選ぶことができる。 HB推定は数値積分を使用する。
Fay-Herriotモデル(後述する)のHB推定に特化したパッケージ。
このパッケージが提供する関数は BayesSAE()である。シンプルでよろしい。
推定量はHBのみ。MCMCを使用する。
小地域推定に限らず、線形・非線形の混合モデルのためのパッケージ。 Rのコアパッケージである。
小地域モデルは基本的には一種の混合モデルなので、 こういう汎用的なパッケージでも推定できる…はずである。
推定量はEBLUP。分散成分の推定量はいくつかの選択肢から選べる。
小地域推定というより、空間データ分析の基本パッケージ。多種多様な機能を 提供している。
空間Fay-Harrietモデルは要するにSAR(同時自己回帰)モデルなので、 こういう汎用的なパッケージでも推定できるはずだ…と思ったんです。
小地域の特性を推測する場合、推定量はEBLUPだということになると 思うのですが… すいません、よくわからないっす。
(本節はRM本4.2節, 6.1節, 10.3節のまとめである。)
Fay-Herriotモデルとは、小地域レベルのデータを使って 小地域特性を推定するためのモデルである。
地域\(i\)(\(=1,\ldots,m\))における、なにかの変数の平均(たとえば一人当たり所得)を\(\bar{Y}_i\)とする。
\(\bar{Y}_i\)そのものについてはわからないが、 その地域の標本から、直接推定量\(\hat{\bar{Y}}_i\)(典型的には標本平均) が手に入っているとしよう。ただし、標本サイズが十分でないので、 この値があてにならない… というのが問題の発端である。
なんらかの関数\(g(\cdot)\)があって、\(\theta_i = g(\bar{Y}_i)\)であるとする。 典型的には\(\theta_i = \bar{Y}_i\)である。 同様に、\(\hat{\theta}_i = g(\hat{\bar{Y}}_i)\)とする。
以上のデータとともに、 地域の補足データ\(\mathbf{z}_i\)(長さ\(p\)のベクトル)も手に入っているとする。 つまり、その地域の\(\theta_i\)を予測するのに役立つような、なんらかの変数群についての 値である。
さて、各地域の\(\bar{Y}_i\)について推測したい。 そこで以下のように想定する。なお以下では、変数\(x\)が平均0, 分散\(\sigma^2\)のなんらかの分布に 従うということを、\(x \sim (0, \sigma^2)\)と略記する。
まず、 \[ \theta_i = \mathbf{z}_i^T \mathbf{\beta} + v_i, \ \ v_i \sim (0, \sigma_v^2)\]
と考える。 つまり、関心の対象である地域特性(\(\theta_i\))は、その地域の補足データ\(\mathbf{z}_i\) で決まっている部分と、それ以外の部分(いわば「地域効果」)\(v_i\)の和だと考えるわけである。 さらに、地域効果はランダムに決まり、その平均は0で、分散は\(\sigma_v^2\)だと考える。
次に、 \[\hat{\theta}_i = \theta_i + e_i, \ \ e_i|\theta_i \sim (0, \psi_i)\]
と考える。つまり、手元にある直接推定量\(\hat{\theta}_i\)は、地域特性の真の値 (\(\theta_i\))を反映してはいるが、標本調査につきものである標本抽出誤差\(e_i\)も 含んでしまっている、と考えるわけである。さらに、各地域の標本抽出誤差は ランダムに決まり、その平均は0で、分散は\(\psi_i\)だと考える。 \(\psi_i\)については既知とする。
このモデルに基づき、\(\theta_i\)について推定することができる。 これがFay-Herriotモデルである。
では、Fay-Herriotモデルのパラメータ (\(\mathbf{\beta}\), \(\sigma_v^2\))をどうやって推定するか。
EBLUP推定の場合は、まず\(\sigma_v^2\)を推定し、次に\(\mathbf{\beta}\)や \(\theta_i\)を推定する。
HB推定の場合は、Fay-Herriotモデルを次の階層ベイズモデルに書き換える。
\[ \begin{aligned} {\rm (i)} \ \ & \hat{\theta}_i | \theta_i, \mathbf{\beta}, \sigma^2_v \mathop{\sim}^{iid} N(\theta_i, \psi_i), \ \ i = 1, \ldots, m \\ {\rm (ii)} \ \ & \theta_i | \mathbf{\beta}, \sigma^2_v \mathop{\sim}^{iid} N(\mathbf{z}_i^T, \sigma_v^2), \ \ i = 1, \ldots, m \\ {\rm (iii)} \ \ & f(\mathbf{\beta}, \sigma^2_v) = f(\mathbf{\beta}) f(\sigma^2_v) \propto f(\sigma^2_v) \end{aligned} \] このモデルを推定するためには、\(\sigma^2_v\)の事前分布\(f(\sigma^2_v)\)を指定する必要がある。
代表的な推定方法として、次の2つがある。
まずはかんたんなデータについて、基本的なFay-Herriotモデルを推定してみよう。
saeパッケージについているサンプルデータ milk を使う。 このデータは、43地域について以下の変数を持っている。
データの出典となっている論文 Arora & Lahiri (1997) をあたってみたところ、これは1989年のデータで、地域は “publication area”だそうである(恥ずかしながら意味がわかりません)。 なお、MajorAreaは現論文にはない。あとから付けたものらしい。
### a.dfMilk: dataset
data(milk, package="sae")
a.dfMilk <- milk %>%
mutate(
SmallArea = as.factor(SmallArea),
MajorArea = as.factor(MajorArea),
var = SD^2
)
rm(milk)
print(a.dfMilk %>% dplyr::select(-var))
各地域の標本サイズは95から633。小地域じゃないじゃん!これで十分じゃん! と私なんかは思うわけですが、そのへんは分野によるんでしょうね。
各地域の支出をチャートにしてみよう。
g <- ggplot(
data = a.dfMilk,
aes(
x = SmallArea,
y = yi,
ymin = yi - SD,
ymax = yi + SD,
color = MajorArea
)
)
g <- g + geom_pointrange()
g <- g + labs(
x = "Small Area",
y = expression(mean %+-% SE),
color = "Major Area"
)
g <- g + theme_bw()
print(g)
上位地域番号によって支出が異なることがわかる。共変量にするとよさそうだ。
(本節はRM本6.5節の再現である。)
まずはsaeパッケージで推定してみよう。
saeパッケージでは、 基本的なFay-Herriotモデルの関数として、 eblupFH()とmseFH()が用意されている。eblupFH()は推定値しか返さないが、 mseFH()は推定量のMSEの推定値も返す、という 違いがあるらしい。 以下ではmseFH()を使うことにする。
まず、共変量なしで。
### a.oModel_saeREML_nox: saeREML model with no covariate
a.oModel_saeREML_nox <- mseFH(
formula = a.dfMilk$yi ~ 1,
vardir = a.dfMilk$SD ^ 2,
method = "REML"
)
cat("refvar:", a.oModel_saeREML_nox$est$fit$refvar, "\n")
refvar: 0.05431129
\(\sigma_v^2\)は0.05431と推定された。
では、共変量として上位地域番号をいれてみよう。
### a.oModel_saeREML: saeREML model
a.oModel_saeREML <- mseFH(
formula = a.dfMilk$yi ~ a.dfMilk$MajorArea,
vardir = a.dfMilk$SD ^ 2,
method = "REML"
)
cat("refvar:", a.oModel_saeREML$est$fit$refvar, "\n")
refvar: 0.01855022
\(\sigma_v^2\)は0.01855と推定された。ぐっと下がっている。
\(\mathbf{\beta}\)は…
print(a.oModel_saeREML$est$fit$estcoef)
試みに、地域効果\(v_i\)を無視した回帰モデルを推定してみよう。
print(summary(lm(yi ~ MajorArea, data = a.dfMilk))$coefficients)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.9854286 0.06937886 14.203585 5.314727e-17
MajorArea2 0.1750000 0.09811653 1.783593 8.227312e-02
MajorArea3 0.2175714 0.08874974 2.451516 1.880886e-02
MajorArea4 -0.2390952 0.08176377 -2.924220 5.725864e-03
固定効果の係数の推定値が異なる。つまり、直接推定量\(\hat{\theta}_i\) を目的変数として、 Fay-Herriotモデル \[\hat{\theta}_i = \mathbf{z}_i^T \mathbf{\beta} + v_i + e_i\] を推定した場合と、 地域効果を無視した回帰モデル \[\hat{\theta_i} = \mathbf{z}^T_i \mathbf{\beta}' + e'_i \] を推定した場合、 \(\mathbf{\beta}\)の推定値と\(\mathbf{\beta}'\)の推定値は異なるのだ。このことは、 RM本6.1.1節の\(\hat{\mathbf{\beta}}\)の式をみてもわかる。 ここ、わたしずーっと誤解してました…反省…
本題に戻って、\(\theta_i\)のEBLUP推定値を見比べてみる。
### a.dfResult_saeREML: result from models by saeREML
a.dfResult_saeREML <- a.dfMilk %>%
mutate(
gEBLUP1 = as.vector(a.oModel_saeREML_nox$est$eblup),
gEBLUP2 = as.vector(a.oModel_saeREML$est$eblup),
gMSE1 = a.oModel_saeREML_nox$mse,
gMSE2 = a.oModel_saeREML$mse,
gCV_0 = CV,
gCV_1 = sqrt(gMSE1)/gEBLUP1,
gCV_2 = sqrt(gMSE2)/gEBLUP2
)
g <- ggplot(
data = a.dfResult_saeREML,
aes(
x = SmallArea,
y = yi,
ymin = yi - SD ,
ymax = yi + SD
)
)
g <- g + geom_pointrange(aes(color = MajorArea))
g <- g + geom_path(aes(y = gEBLUP1, group = 0))
g <- g + geom_path(aes(y = gEBLUP2, group = 0, color = MajorArea))
g <- g + labs(
x = "Small Area",
y = expression(mean %+-% SE),
color = "Major Area"
)
g <- g + theme_bw()
print(g)
折れ線のうち、黒いのは共変量を使わなかった場合、色がついているのは 使った場合。前者は全体平均に、後者はMajor rea別の 平均に値が寄っている。
今度はMSEの推定値を比べてみよう。
平均が大きいとそのMSEも大きくなりそうなので、 変動係数で比較する。
dfPlot <- a.dfResult_saeREML %>%
arrange(-ni) %>%
mutate(
fSmallArea = factor(
seq_along(SmallArea),
labels = as.character(SmallArea)
)
) %>%
dplyr::select(fSmallArea, MajorArea, ni, starts_with("gCV_")) %>%
gather(sVar, gValue, starts_with("gCV_")) %>%
separate(sVar, c("sVar1", "sVar2")) %>%
mutate(
fMethod = factor(
as.integer(sVar2),
levels = 0:2,
labels = c("Direct", "EBLUP(w/o covariate)", "EBLUP(w/ covariate)")
)
)
g <- ggplot(
data = dfPlot,
aes(
x = fSmallArea,
y = gValue,
group = fMethod,
color = fMethod
)
)
g <- g + geom_path()
g <- g + labs(
x = "Small Area (sorted by decreasing sample size)",
y = "Coef.Variation",
color = "Method"
)
g <- g + theme_bw()
print(g)
今度は横軸を標本サイズ順に並べている。
直接推定量、共変量なしのEBLUP, 共変量をいれたEBLUPと進むにつれて 変動係数が小さくなっている。 また、標本サイズが大きい地域では、変動係数があまり減っていない。なるほどね。
なお、以上は\(\sigma_v^2\)をREML推定した場合であった。 後述する理由により、ここでML推定も試しておく。
### a.oModel_saeML: model by saeML
a.oModel_saeML <- mseFH(
formula = a.dfMilk$yi ~ a.dfMilk$MajorArea,
vardir = a.dfMilk$SD ^ 2,
method = "ML"
)
print(a.oModel_saeML$est$fit$refvar)
[1] 0.01551755
\(\sigma_v^2\)は0.01552と推定された。ちょっぴり低めですね。
sub_PlotTwoEstimate.1(
agEst1 = as.vector(a.oModel_saeREML$est$eblup),
agMSE1 = a.oModel_saeREML$mse,
agEst2 = as.vector(a.oModel_saeML$est$eblup),
agMSE2 = a.oModel_saeML$mse,
asLabel = a.dfMilk$SmallArea,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "sae (ML)"
)
\(\theta_i\)の推定値も、推定量のMSEの推定値も、ほとんどかわらない。
共変量をいれたモデルの推定を、今度はhbsaeパッケージで試してみる。
hbsaeパッケージは推定量として以下を御用意している。
まずはEBLUP推定を試してみよう。saeパッケージの出力と揃えるため、REMLを選ぶ。
### a.oModel_hbsaeREML: model by hbsaeREML
a.oModel_hbsaeREML <- fSAE.Area(
a.dfMilk$yi,
a.dfMilk$SD^2,
model.matrix(yi ~ MajorArea, data=a.dfMilk),
method = "REML"
)
REML estimate of variance ratio: 0.01855
sub_PlotTwoEstimate.1(
agEst1 = as.vector(a.oModel_saeREML$est$eblup),
agMSE1 = a.oModel_saeREML$mse,
agEst2 = EST(a.oModel_hbsaeREML),
agMSE2 = SE(a.oModel_hbsaeREML)^2,
asLabel = a.dfMilk$SmallArea,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "hbsae (REML)"
)
\(\sigma_v^2\)は, hbsaeのfSAE.Area()ではvariance ratioと呼ばれるらしい。 なんだか変な呼び方だが、 きっとfSAE.Area()の内側でユニットレベルモデルの関数をコールしている から、ユニットレベルモデルの用語がそのまま残っているのだろう。 推定値は0.01855、saeパッケージの推定値とほぼ同じだ。
\(\theta_i\)の推定値もほぼ同じ。しかし、推定量のMSEの推定値は hbsaeのほうが少し低めになっている。なぜかわからない。
こんどは階層ベイズ推定。 hbsaeパッケージにおけるHB推定は、MCMCではなく数値積分を使っている。
### a.oModel_hbsaeHB: model by hbsaeHB
set.seed(12345)
a.oModel_hbsaeHB <- fSAE.Area(
a.dfMilk$yi,
a.dfMilk$SD^2,
model.matrix(yi ~ MajorArea, data=a.dfMilk)
)
REML estimate of variance ratio: 0.01855
numerical integration of f(x) (normalization constant): 2.231 with absolute error < 5.1e-10
numerical integration of x*f(x): 0.05055 with absolute error < 6.6e-11
posterior mean for variance ratio: 0.02266
sub_PlotTwoEstimate.1(
agEst1 = as.vector(a.oModel_saeREML$est$eblup),
agMSE1 = a.oModel_saeREML$mse,
agEst2 = EST(a.oModel_hbsaeHB),
agMSE2 = SE(a.oModel_hbsaeHB)^2,
asLabel = a.dfMilk$SmallArea,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "hbsae (HB)"
)
\(\sigma_v^2\)の推定値(事後平均)は0.02266。少し高くなった。 \(\theta_i\)の推定値は、saeパッケージとほぼ同じだが、MSEの推定値はちょっとずれている。
こんどはBayesSAEパッケージで試してみよう。 MCMCによるHB推定である。
事前分布の指定がデフォルトでどうなっているのか、 マニュアルを読んでもよくわからなかったので、 ソースコードを調べて突き止めた。以下のコードでprior=に与えている リスト これがデフォルトである(従ってprior=オプションを消しても 結果は同じ)。\(\beta\)は無情報事前分布, \(\sigma_v^2\)は \((0, 100)\)の一様分布。
set.seed(123)
a.oModel_BayesSAE <- BayesSAE::BayesSAE(
yi ~ MajorArea|var,
data=a.dfMilk,
mcmc=20000,
prior=list(
beta.type = "non_in",
beta.prior = 0,
sigv.type = "unif",
sigv.prior = list(eps2=1e-02)
)
)
print(summary(a.oModel_BayesSAE))
Call:
BayesSAE::BayesSAE(formula = yi ~ MajorArea | var, prior = list(beta.type = "non_in",
beta.prior = 0, sigv.type = "unif", sigv.prior = list(eps2 = 0.01)),
mcmc = 20000, data = a.dfMilk)
Basic Fay-Herriort Model
Sampling model: yi ~ theta
Linking model: theta ~ MajorArea + var + NA
Rao-Blackwellization of theta's based on the simulation:
1 2 3 4 5 6 7 8 9 10
1.0257 1.0489 1.0698 0.7530 0.8406 0.9745 1.0672 1.0983 1.2321 1.2044
11 12 13 14 15 16 17 18 19 20
0.7759 1.2266 1.2175 0.9772 1.1861 1.1535 1.2282 1.2932 1.2388 1.2382
21 22 23 24 25 26 27 28 29 30
1.0854 1.1923 1.1175 1.2255 1.1939 0.7646 0.7669 0.7356 0.7716 0.6108
31 32 33 34 35 36 37 38 39 40
0.7751 0.8035 0.7745 0.6088 0.6994 0.7611 0.5250 0.7446 0.7559 0.7721
41 42 43
0.7486 0.8078 0.6790
Coefficients in the linking model:
Mean SD 2.5% 97.5%
beta.1 0.97049 0.07142 0.82950 1.110
beta.2 0.13513 0.10713 -0.06956 0.351
beta.3 0.22345 0.09577 0.03532 0.416
beta.4 -0.24317 0.08397 -0.40953 -0.076
Variance of residual in the linking model:
Mean SD 2.5% 97.5
[1,] 0.022241 0.009261 0.008402 0.044
DIC: -28.35111
\(\sigma_v^2\)の推定値(事後平均)は0.022241。hbsaeのHB推定に近い。
sub_PlotTwoEstimate.1(
agEst1 = as.vector(a.oModel_saeREML$est$eblup),
agMSE1 = a.oModel_saeREML$mse,
agEst2 = a.oModel_BayesSAE$HB,
agMSE2 = summary(MCMC(a.oModel_BayesSAE))$statistics[1:43,2]^2,
asLabel = a.dfMilk$SmallArea,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "BayesSAE"
)
\(\theta_i\)の推定値、ならびにMSEの推定値も、hbsaeのHB推定に近いようだ。
混合モデルの汎用パッケージであるnlmeパッケージで試してみよう。 正規混合モデルの関数lme()を使う。\(\theta_i\)の推定量はEBLUP。 \(\sigma^2_v\)の推定量はデフォルトではREMLになる。
a.oModel_lmeREML <- lme(
fixed = yi ~ MajorArea,
random = ~ 1 | SmallArea,
control = lmeControl(sigma = 1, apVar = T),
weights = varFixed(~var),
data = a.dfMilk
)
print(a.oModel_lmeREML)
Linear mixed-effects model fit by REML
Data: a.dfMilk
Log-restricted-likelihood: 10.34588
Fixed: yi ~ MajorArea
(Intercept) MajorArea2 MajorArea3 MajorArea4
0.9680768 0.1316132 0.2269008 -0.2415905
Random effects:
Formula: ~1 | SmallArea
(Intercept) Residual
StdDev: 0.1332918 1
Variance function:
Structure: fixed weights
Formula: ~var
Number of Observations: 43
Number of Groups: 43
lme()関数のweights=varFixed(~ var)で、「グループ内の分散は変数varである」 と指定していることになるのだそうだ(ここでグループとは小地域のことであろう)。 また、control=オプションのうち、sigma=1とは「残差分散を1に固定する」という意味、 apVar = Tとは「分散共分散行列を求める」という意味らしい。 こうすると残差分散は var*1 になる、ということなのかな。
へーこうやって書くのか、勉強になるなあ、と思ったのだが、 よくよくみると、結果がちょっぴりsaeパッケージと異なる。 ここでは\(\sigma_v\)の推定値が 0.1332918 となっている。二乗すると0.0177667。 saeパッケージでは0.01855であった。あれ?ちょっと低いぞ?
実はこの話、nlmeパッケージのバグとして報告されている。 saeパッケージとぴったり合わないのはおかしいのだそうで、これはどうやら nlme側のバグらしい。うわー。まさかそんなことがー。
なお、ML推定を指定した場合の出力は正しいそうだ。
a.oModel_lmeML <- lme(
fixed = yi ~ MajorArea,
random = ~ 1 | SmallArea,
control = lmeControl(sigma = 1, apVar = T),
weights = varFixed(~var),
method = "ML",
data = a.dfMilk
)
print(a.oModel_lmeML)
Linear mixed-effects model fit by maximum likelihood
Data: a.dfMilk
Log-likelihood: 12.77117
Fixed: yi ~ MajorArea
(Intercept) MajorArea2 MajorArea3 MajorArea4
0.9677986 0.1278755 0.2266909 -0.2425804
Random effects:
Formula: ~1 | SmallArea
(Intercept) Residual
StdDev: 0.1245693 1
Variance function:
Structure: fixed weights
Formula: ~var
Number of Observations: 43
Number of Groups: 43
0.1245693^2 = 0.01551751, saeパッケージでのML推定とほぼ一致する。
sub_PlotTwoEstimate.1(
agEst1 = as.vector(a.oModel_saeREML$est$eblup),
agMSE1 = NULL,
agEst2 = predict(a.oModel_lmeML),
agMSE2 = NULL,
asLabel = a.dfMilk$SmallArea,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "nlme (ML)"
)
推定値はsaeパッケージとたいして変わらない。なお、推定量のMSEの推定値のほうは、 nlmeパッケージだけでは出力できず、なにか工夫しないといけない模様である。
最後に、ベイジアン・モデリングのための汎用パッケージrstanで試してみよう。
stanコードはこちら。 わたくし、普段はほぼ使わないので、正直なところ全然自信がないんですが…
// Fay-Herriot model
// cf. https://rdrr.io/github/ssoro411/bayesae/src/R/Models_fast.R
data {
int<lower=0> M;
int<lower=0> P;
vector[M] Y;
vector<lower=0>[M] SD;
matrix[M,P] X;
}
parameters {
real<lower=0> sigma;
vector[P] beta;
vector[M] theta;
}
transformed parameters {
vector[M] mu;
mu = X*beta;
}
model {
theta ~ normal(mu, sigma);
Y ~ normal(theta, SD);
}
では実行します…
# a.oStanModel <- stan_model(file="./FH1_vectorized.stan")
# save(a.oStanModel, file="./a.oStanModel.RData")
#
# load(file = "./a.oStanModel.RData")
#
# mnModelMatrix <- model.matrix(~ MajorArea, data=a.dfMilk)
# lData <- list(
# M = nrow(a.dfMilk),
# P = ncol(mnModelMatrix),
# Y = a.dfMilk$yi,
# SD = a.dfMilk$SD,
# X = mnModelMatrix
# )
# a.oStanFit <- sampling(
# a.oStanModel,
# data = lData,
# seed = 1234,
# control = list(adapt_delta = 0.95)
# )
# save(a.oStanFit, file="./a.oStanFit.RData")
load(file = "./a.oStanFit.RData")
stan_trace(a.oStanFit, c("sigma", "theta[28]"), inc_warmup = T)
stan_rhat((a.oStanFit))
\(\sigma_v\)と\(\theta_{28}\)のトレースプロット(地域28は推定値のSEがいちばん大きい地域)、 ならびに\(\hat{R}\)のヒストグラムを示す。 余裕で収束している模様。
地域効果の推定値ならびにMSEの推定と、saeパッケージと比較してみると…
lResult <- rstan::extract(a.oStanFit)
agEst <- colMeans(lResult$theta)
agMSE <- apply(lResult$theta, 2, var)
sub_PlotTwoEstimate.1(
agEst1 = as.vector(a.oModel_saeREML$est$eblup),
agMSE1 = a.oModel_saeREML$mse,
agEst2 = agEst,
agMSE2 = agMSE,
asLabel = a.dfMilk$SmallArea,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "rstan"
)
推定値のほうはほぼ同じ。MSEは少しずれていて、hbsaeパッケージ(HB推定)や、 BayesSAEパッケージと似ている。
rm(list=grep("^a\\.", ls(), value=T))
こうしてみてみると、いろいろ疑問が湧いてくる。わたくしと致しましては、 特に次の3点に関心がある。
そこでシミュレーションしてみることにしました。
いま、地域\(i\)(\(=1,\ldots,m\))のひとりひとりの住民\(j\)(\(=1,\ldots,N_i\))が、 量的な反応変数\(y_{ij}\)を持っている。えーっと、幸福度ってことにしましょうか。
地域の母集団サイズ\(N_i\)はいずれも十分に大きい。 各地域の幸福度の母平均\(\bar{Y}_i = \frac{1}{N_i} \sum_j y_{ij}\)について推測したい。 そこで、各地域\(i\)からサイズ\(n_i\)の標本を単純無作為抽出し、 標本平均\(\hat{\bar{Y}}_i = \frac{1}{n_i} \sum_j y_{ij}\), ならびにその分散の推定値 \(\hat{Var}(\hat{\bar{Y}})_i = \frac{1}{n_i(n_i-1)} \sum_j (y_{ij} - \hat{\bar{Y}}_i)\) を求めた。
ところが、地域ごとの標本サイズ\(n_i\)は必ずしも十分でない。 標本平均\(\hat{\bar{Y}}_i\)を\(\bar{Y}_i\)の推定量とするのは、いささか心許ない。 さらに、各地域は\(p+1\)個のブロックに分類されており、 ブロックによって幸福度平均\(\bar{Y}_i\)の高さが異なる、 と私は信じている。
そこで私はこう考える。幸福度の標本平均\(\hat{\bar{Y}}_i\)だけでなく、 その地域がどのブロックに属するかという情報もつかって、 各地域の幸福度の母平均\(\bar{Y}_i\)を 推測しよう!
というわけで、Fay-Erriotモデルの登場である。
以下、記号が煩雑になるので、 RM本にあわせて\(\theta_i = \bar{Y}_i\), \(\hat{\theta}_i = \hat{\bar{Y}}_i\)と書く。 また、\(\hat{\psi}_i = \hat{Var}(\hat{\bar{Y}}_i)\)と書く。さらに、 ある地域がどこに分類されるかを、長さ\(p\)の二値ベクトル\(\mathbf{z}_i\)で表す。
おさらいすると、Fay-Herriotモデルは以下の通り。 \[ \hat{\theta}_i = \mathbf{z}^T_i \mathbf{\beta} + v_i + e_i, \ \ v_i \sim (0, \sigma^2_v), \ \ e_i \sim (0, \psi_i) \]
以下では\(v_i\)を地域効果、\(e_i\)を標本抽出誤差と呼ぼう。
さて、これから分析用のデータを生成します。
データ生成モデルは次のとおり。\(i=1,\ldots,m\)について、 \[ y_{ij} = 50 + \mathbf{z}^T_i \mathbf{0} + v_i + \epsilon_{ij}, \ \ v_i \mathop{\sim}^{iid} N(0, \sigma^2_v), \ \ \epsilon_{ij} \mathop{\sim}^{iid} N(0,100) \]
つまり、私は知らないのだが、実は
母集団サイズは大きく、\(\epsilon_{ij}\)の母集団平均は 0とみなせるので、結局、知りたい特性\(\theta_i\)はここでは\(v_i\)に等しい。
地域の標本サイズ\(n_i\)は、地域1から地域10まで順に\(10, 20, ..., 100\)とし、 地域11以降もこれを繰り返す。
共変量\(\mathbf{z}_i\)(長さ\(p\))は、地域のブロック(\(p+1\)水準)を表すものとする。 ブロック番号は、地域1から地域\(p+1\)までは順に \(1,...,p+1\)とし、 残りの地域については等確率でランダムに決める。
…というデータを生成し分析する、一連の関数を作りました。
### for rstan ...
# b.oStanModel <- stan_model(file="./FH1_vectorized.stan")
# save(b.oStanModel, file="./b.oStanModel.RData")
load(file = "./b.oStanModel.RData")
b.sub_GenerateDataset1.1 <- function(
nNumArea, nNumGroup, gVarV, gVarE=100, anPossibleN=seq(10, 100, by = 10)
){
## purpose: generate a dataset for FH model
## arg:
## nNumArea: number of areas
## nNumGroup: number of levels of a categorical covariate
## gVarV: variance of area effects
## anPossibleN: vector of possible sample size
## return:
## a data frame, area level data
## $nArea: Area ID
## $nSize: sample size
## $fGroup: categorical covariate with (nNumGroup) levels
## $gEst_Direct: sample means of a response variable
## $gMSE_Direct: sample est. of variance of sample means
## $gV: area effect
out <- data.frame(
nArea = 1:nNumArea,
nSize = rep(anPossibleN, nNumArea %/% length(anPossibleN) + 1)[1:nNumArea],
fGroup = factor(
c(1:nNumGroup, sample(1:nNumGroup, nNumArea - nNumGroup, replace = T)),
levels = 1:nNumGroup
),
gV = 50 + rnorm(n = nNumArea, mean = 0, sd = sqrt(gVarV))
)
lError <- lapply(
out$nSize,
function(nSize){
rnorm(n = nSize, mean = 0, sd = sqrt(gVarE))
}
)
out$gEst_Direct <- out$gV + sapply(lError, mean)
out$gMSE_Direct <- sapply(lError, function(x){ var(x) / length(x) })
out$nTotalSize <- sum(out$nSize)
return(out)
}
b.sub_GenerateMultiDataset1.1 <- function(nNumDataset=1, ...){
## purpose: generate multiple datasets for FH model
## arg:
## nNumDataset: number of data sets
## ...: pass to sub_GenerateDataset1.1()
## return:
## a data frame, area level data
## $nDataset: Dataset ID
## $... return from sub_GenerateDataset1.1()
out <- data.frame(nDataset = 1:nNumDataset) %>%
group_by(nDataset) %>%
dplyr::do(b.sub_GenerateDataset1.1(...)) %>%
ungroup()
return(out)
}
b.sub_EstimateDatasetByMethod1.1 <- function(
dfData, sMethod, nTimeOut = 1800, bVerbose=FALSE
){
## purpose: estimate FH model for a dataset by a method
## arg:
## dfData: a data set
## sMethod: {"saeREML", "saeREMLNOMSE", "hbsaeREML", "hbsaeHB",
## "bayessae", "nlme", "rstan"}
## nTimeOut: time limit
## bVerbose: verbose flag (logical)
## return:
## a list
## notes:
## - nTimeOut is effective only for hbsaeREML and hbsaeHB
## - rstan requires b.oStanModel in the environment
## trap: check sMethod
stopifnot(
sMethod %in% c(
"saeREML", "saeREMLNOMSE", "hbsaeREML", "hbsaeHB",
"bayessae", "nlme", "rstan"
)
)
if (bVerbose) cat("[sub_EstimateDataset1.1]", sMethod, "... \n")
t <- proc.time()
if (sMethod == "saeREML") {
oModel <- try(
mseFH(
formula = dfData$gEst_Direct ~ dfData$fGroup,
vardir = dfData$gMSE_Direct,
method = "REML"
)
)
if (!inherits(oModel, "try-error")) {
agEst <- as.vector(oModel$est$eblup)
agMSE <- oModel$mse
}
}
if (sMethod == "saeREMLNOMSE") {
oModel <- try(
eblupFH(
formula = dfData$gEst_Direct ~ dfData$fGroup,
vardir = dfData$gMSE_Direct,
method = "REML"
)
)
if (!inherits(oModel, "try-error")) {
agEst <- as.vector(oModel$eblup)
agMSE <- rep(NA, nrow(dfData))
}
}
if (sMethod == "hbsaeREML") {
oModel <- withTimeout(
try(
fSAE.Area(
dfData$gEst_Direct,
dfData$gMSE_Direct,
model.matrix(gEst_Direct ~ fGroup, data = dfData),
method = "REML",
silent = TRUE
)
),
timeout = nTimeOut
)
if (!inherits(oModel, "try-error")) {
agEst <- EST(oModel)
agMSE <- SE(oModel)^2
}
}
if (sMethod == "hbsaeHB") {
oModel <- withTimeout(
try(
fSAE.Area(
dfData$gEst_Direct,
dfData$gMSE_Direct,
model.matrix(gEst_Direct ~ fGroup, data = dfData),
method = "HB",
silent = TRUE
)
),
timeout = nTimeOut
)
if (!inherits(oModel, "try-error")) {
agEst <- EST(oModel)
agMSE <- SE(oModel)^2
}
}
if (sMethod == "bayessae") {
oModel <- try(
BayesSAE(
gEst_Direct ~ fGroup | gMSE_Direct,
data = dfData,
mcmc = 5000
)
)
if (!inherits(oModel, "try-error")) {
agEst <- oModel$HB
agMSE <- summary(MCMC(oModel))$statistics[seq_along(dfData$nArea),2]^2
}
}
if (sMethod == "nlme") {
oModel <- try(
lme(
fixed = gEst_Direct ~ fGroup,
random = ~ 1 | nArea,
control = lmeControl(sigma = 1, apVar = T),
weights = varFixed(~gMSE_Direct),
method = "ML",
data = dfData
)
)
if (!inherits(oModel, "try-error")) {
agEst <- fitted(oModel)
agMSE <- rep(NA, nrow(dfData))
}
}
if (sMethod == "rstan") {
## trap: b.oStanModel is found in the enviroment
stopifnot(exists("b.oStanModel"))
mnModelMatrix <- model.matrix(~ fGroup, data = dfData)
lData <- list(
M = nrow(dfData),
P = ncol(mnModelMatrix),
Y = dfData$gEst_Direct,
SD = sqrt(dfData$gMSE_Direct),
X = mnModelMatrix
)
oModel <- try(
sampling(
b.oStanModel,
data = lData,
seed = 1234,
iter = 400,
control = list(adapt_delta = 0.95),
open_progress = FALSE,
show_messages = FALSE
)
)
if ( !inherits(oModel, "try-error") ) {
if ( any(summary(oModel)$summary[,"Rhat"] > 1.1) ) {
oModel <- try(stop("not converged!"))
}
}
if ( !inherits(oModel, "try-error") ) {
lResult <- rstan::extract(oModel)
agEst <- colMeans(lResult$theta)
agMSE <- apply(lResult$theta, 2, var)
}
}
## common process ...
if (!inherits(oModel, "try-error")) {
out <- list(
sMethod = sMethod,
bStatus = TRUE,
agEst = agEst,
agMSE = agMSE,
gTime = (proc.time() - t)[3]
)
if (bVerbose) cat("[sub_EstimateDataset1.1] Success.", out$gTime, "sec.\n")
} else {
out <- list(
sMethod = sMethod,
bStatus = FALSE,
agEst = rep(NA, nrow(dfData)),
agMSE = rep(NA, nrow(dfData)),
gTime = (proc.time() - t)[3]
)
if (bVerbose) cat("[sub_EstimateDataset1.1] Fail.", out$gTime, "sec.\n")
}
return(out)
}
b.sub_EstimateDataset1.1 <- function(
dfData,
asMethod=c("saeREML", "hbsaeREML", "hbsaeHB", "bayessae", "nlme", "rstan"),
...
){
## purpose: estimate FH model for a dataset by multiple methods
## arg:
## dfData: a dataset
## asMethod
## ...
## return:
## a data frame
lOut <- lapply(
asMethod,
function(sMethod){
lResult <- b.sub_EstimateDatasetByMethod1.1(dfData, sMethod = sMethod, ...)
out <- data.frame(
nArea = dfData$nArea,
sMethod = rep(lResult$sMethod, nrow(dfData)),
bStatus = rep(lResult$bStatus, nrow(dfData)),
gEst = lResult$agEst,
gMSE = lResult$agMSE,
gTime = rep(lResult$gTime, nrow(dfData)),
stringsAsFactors = FALSE
)
return(out)
}
)
out <- bind_rows(lOut) %>%
gather(sVar, gValue, c(bStatus, gEst, gMSE, gTime)) %>%
mutate(sVar = paste0(sVar, "_", sMethod)) %>%
dplyr::select(nArea, sVar, gValue) %>%
spread(sVar, gValue)
out <- full_join(dfData, out, by = "nArea")
return(out)
}
b.sub_EstimateMultiDataset1.1 <- function(dfData, ...){
## purpose: estimate FH models for multiple datasets by multiple methods
## arg:
## dfData: multiple datasets
## return:
## a data frame of multiple results
out <- dfData %>%
group_by(nDataset) %>%
dplyr::do(b.sub_EstimateDataset1.1(., ...)) %>%
ungroup()
return(out)
}
b.sub_EstimatorAsFactor.1 <- function(asEstimator){
## purpose: convert estimator names into factor
## args:
## asEstimator: a character vector of estimator names
## return:
## a factor vector
factor(
asEstimator,
levels = c(
"True", "Direct", "saeREML", "saeREMLNOMSE",
"hbsaeREML", "hbsaeHB", "bayessae", "nlme", "rstan"
),
labels = c(
"True", "Direct", "sae(REML)", "sae(REML,noMSE)",
"hbsae(REML)", "hbsae(HB)", "BayesSAE", "nlme", "rstan"
)
)
}
b.sub_EstimatorPalette.1 <- function(afEstimator){
## purpose: get consisntent color pallette for estimators
## args:
## afEstimator: a factor vector of estimator
## return:
## a vector of color palette
gg_color_hue(9)[sort(unique(as.integer(afEstimator)))]
}
ためしに、上記のデータ生成関数とデータ分析関数を使ってみよう。
まずデータ生成。10地域、4ブロック、地域効果の分散は50とする。 地域効果の分散は、地域内の個人差の半分の大きさだ。 結構大きな地域差があるんですね。
## b.dfDemoData: data for demo
set.seed(100)
b.dfDemoData <- b.sub_GenerateDataset1.1(nNumArea=10L, nNumGroup=4L, gVarV=50)
print(b.dfDemoData)
では、このデータについて推定し、次の8個を比較する。
dfResult <- b.sub_EstimateDataset1.1(b.dfDemoData) %>%
dplyr::select(nArea, nSize, gV, starts_with("gEst_"))
dfPlot <- dfResult %>%
dplyr::select(nArea, gV, starts_with("gEst_")) %>%
rename(gEst_True = gV) %>%
gather(sVar, gValue , starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sVar2")) %>%
mutate(fEstimator = b.sub_EstimatorAsFactor.1(sVar2))
g <- ggplot(
data=dfPlot,
aes(x=as.factor(nArea), y=gValue, group=fEstimator, color=fEstimator)
)
g <- g + geom_path(size=1)
g <- g + labs(
x = "area",
y = expression(paste(theta, " or ", hat(theta))),
color = "Estimator"
)
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
縦軸は真値ないし推定値である。
たとえば地域1の正解は56.3点。標本平均は59.2点, 正解よりもかなり高めになってしまった (地域1の標本サイズは10しかない)。 しかるに、saeパッケージの推定値は57.5点。標本平均よりも正解に近くなっている。
もっとも、この結果はたまたまかもしれない。 同じことを何度も繰り返してみないと、なんともいえない。
おまたせしました。なんども繰り返してみます。
先程と同じく、10地域、4ブロック、地域効果の分散は50とする。 500回繰り返してみた。40分ほどかかりました。
# set.seed(123456)
# dfData <- b.sub_GenerateMultiDataset1.1(
# nNumDataset = 500, nNumArea = 10L, nNumGroup = 4L, gVarV = 50
# )
# b.dfSimulation1 <- b.sub_EstimateMultiDataset1.1(dfData)
# save(b.dfSimulation1, file="./b_dfSimulation1.RData")
load("./b_dfSimulation1.RData")
500回中1回、rstanが収束しなかった試行 (いずれかのパラメータの\(\hat{R}\)が1.10を上回わった試行)があった。 これは反復数が足りなかったせいだと思うので、 公平を期するため、以下では1試行まるごと除外する。つまり、他の推定量による 推定結果も除外する。
推定値はどのくらいあてになるか、調べてみよう。 おそらく、地域の標本サイズ\(n_i\)が小さいとき、推定値は正解からずれてしまうだろうから、 \(n_i\)別に見る必要がありそうだ。
dfPlot <- b.dfSimulation1 %>%
filter(!is.na(gEst_rstan)) %>%
dplyr::select(nSize, gV, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
mutate(
gDiff = gEst - gV,
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator)
)
g <- ggplot(data=dfPlot, aes(x = nSize, y=gDiff, color=fEstimator))
g <- g + geom_point(alpha=1/5)
g <- g + labs(
x = "Sample Size",
y = expression(hat(theta)-theta),
color = "Estimator"
)
g <- g + geom_hline(yintercept = 0)
g <- g + scale_x_continuous(breaks=seq(10, 100, by = 20))
g <- g + scale_color_manual(
guide=FALSE,
values = b.sub_EstimatorPalette.1(dfPlot$fEstimator)
)
g <- g + facet_grid(~ fEstimator)
g <- g + theme_bw()
print(g)
縦軸は、推定値と正解との差を表している。 やはり、標本サイズが小さいときに、正解からのずれが大きくなる。
推定方法 x 標本サイズ別に、推定値と正解との差を集計してみよう。
dfPlot <- b.dfSimulation1 %>%
filter(!is.na(gEst_rstan)) %>%
dplyr::select(nSize, gV, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(sEstimator, nSize) %>%
summarize(
nFreq = n(),
gRMSE = sqrt(mean((gEst - gV)^2))
) %>%
ungroup() %>%
mutate(fEstimator = b.sub_EstimatorAsFactor.1(sEstimator))
g <- ggplot(data=dfPlot, aes(x = nSize, y=gRMSE, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + labs(x = "Sample Size", y = "observed RMSE", color="Estimator")
g <- g + scale_x_continuous(breaks=seq(10, 100, by = 10))
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
縦軸は「推定値と正解との差の二乗の平均の平方根」(RMSE)を示す。平たく言っちゃうと、 推定した地域平均は、真の地域平均からだいたい何点くらいずれるか、を示している。
標本サイズが10のときには、だいたい3点くらいずれてしまう。 標本サイズが大きくなると、ずれは小さくなる。
ここで関心があるのは推定方法の間の差だ。 上のチャートだと差を読み取れないので描き直そう。
### observed RMSE / theoretical RMSE of Direct
dfPlot <- b.dfSimulation1 %>%
filter(!is.na(gEst_rstan)) %>%
dplyr::select(nSize, gV, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(sEstimator, nSize) %>%
summarize(
gRMSE = sqrt(mean((gEst - gV)^2))
) %>%
ungroup() %>%
mutate(
gRatio = gRMSE / sqrt(100/nSize),
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator)
)
g1 <- ggplot(data=dfPlot, aes(x = nSize, y=gRatio, group=fEstimator, color=fEstimator))
g1 <- g1 + geom_path(size=1)
g1 <- g1 + geom_hline(yintercept = 1)
g1 <- g1 + labs(
x = "Sample Size",
y = "(observed RMSE)/(theoretical RMSE of Direct)",
color = "Estimator"
)
g1 <- g1 + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g1 <- g1 + scale_x_continuous(breaks=seq(10, 100, by = 10))
g1 <- g1 + theme_bw()
print(g1)
### (observed RMSE) / (observed RMSE of Direct)
# dfTemp <- b.dfSimulation1 %>%
# filter(!is.na(gEst_rstan)) %>%
# dplyr::select(nSize, gV, starts_with("gEst_")) %>%
# gather(sVar, gEst, starts_with("gEst_")) %>%
# separate(sVar, c("sVar1", "sEstimator")) %>%
# group_by(sEstimator, nSize) %>%
# summarize(
# gRMSE = sqrt(mean((gEst - gV)^2))
# ) %>%
# ungroup()
# dfTemp.1 <- dfTemp %>% filter(sEstimator != "Direct")
# dfTemp.2 <- dfTemp %>%
# filter(sEstimator == "Direct") %>%
# rename(gRMSE_Direct = gRMSE) %>%
# dplyr::select(-sEstimator)
# dfPlot <- full_join(dfTemp.1, dfTemp.2, by = "nSize") %>%
# mutate(
# gRatio = gRMSE / gRMSE_Direct,
# fEstimator = b.sub_EstimatorAsFactor.1(sEstimator)
# )
# g1 <- ggplot(data=dfPlot, aes(x = nSize, y=gRatio, group=fEstimator, color=fEstimator))
# g1 <- g1 + geom_path(size=1)
# g1 <- g1 + geom_hline(yintercept = 1)
# g1 <- g1 + labs(
# x = "Sample Size",
# y = "(RMSE)/(RMSE of Direct)",
# color = "Estimator"
# )
# g1 <- g1 + scale_x_continuous(breaks=seq(10, 100, by = 10))
# g1 <- g1 + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
# g1 <- g1 + theme_bw()
# print(g1)
g2 <- ggplot(
data=dfPlot %>% filter(nSize == 10),
aes(x = fEstimator, y = gRatio, color=fEstimator)
)
g2 <- g2 + geom_point(size=5)
g2 <- g2 + geom_hline(yintercept = 1)
g2 <- g2 + labs(
x = "Estimator",
y = "(observed RMSE)/(theoretical RMSE of Direct), size=10"
)
g2 <- g2 + scale_color_manual(
guide = FALSE,
values = b.sub_EstimatorPalette.1(dfPlot$fEstimator)
)
g2 <- g2 + theme_bw()
print(g2)
# grid.arrange(g1, g2, ncol = 2)
さきほど求めたRMSEはデータから得たものだが、 標本抽出誤差は実は100なので、これを既知とするならば、 標本平均のRMSEは\(\sqrt{100/n_i}\)である。 上の2枚の図の縦軸は、この理屈上の値を基準に、それぞれの推定方法による推定値の RMSEを比で表している(つまり、長い目でみれば標本平均(Direct)の折れ線は 縦軸1の位置の水平線になるはずで、そうなっていないのは運の問題である)。 下にいくほど、推定値と真の値のずれが小さい。
上の図は全地域。下の図は標本サイズ10の地域のみ拡大している。
以下の点がみてとれる。
こんどはMSEの推定値について調べてみよう。
推定量のMSEとは、「仮にある推定値を得て、その値を真の値と比べることができるなら、 差の二乗の期待値がどうなるか」という値である。 この実験では真の値が分かっているのだから、 推定値と真値の差の二乗そのものが手に入る。この実績を、 推定されたMSEと比べてみようではないか。
dfPlot <- b.dfSimulation1 %>%
filter(!is.na(gEst_rstan)) %>%
dplyr::select(nDataset, nArea, nSize, gV, starts_with("gEst_"), starts_with("gMSE_")) %>%
gather(sVar, gValue, c(starts_with("gEst_"), starts_with("gMSE_"))) %>%
separate(sVar, c("sVar", "sEstimator")) %>%
spread(sVar, gValue) %>%
filter(sEstimator != "nlme") %>%
mutate(
gDiff = gMSE - (gEst - gV)^2,
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator)
)
g <- ggplot(data=dfPlot, aes(x = nSize, y = gDiff, color = fEstimator))
g <- g + geom_point(alpha=1/5)
g <- g + labs(
x = "Sample Size",
y = "(Estimated MSE)-(Actual Sq.Err)"
)
g <- g + geom_hline(yintercept = 0)
g <- g + scale_x_continuous(breaks=seq(10, 100, by = 20))
g <- g + facet_grid(~ fEstimator)
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
縦軸は、「推定されたMSEから、推定値と正解との差の二乗を引いた値」である。 縦軸が0だということは、 推定値と正解との間にはずれがあったかもしれないが、そのずれの大きさは予想通りだった、 ということを表している。正の値は「思ったよりずれなかった」つまり MSEが過大に推定されていたことを表し、 負の値は「思ったよりずれた」つまりMSEが過小に推定されていたことを表す。 なお、nlmeパッケージがないのは、すでに述べたように、 nlmeパッケージがEBLUP推定量のMSEの推定を提供していないからである。
標本サイズが小さい時、MSEの推定は 過小だったり過大だったりしやすくなることが見て取れる。
では、推定方法 x 標本サイズ別に、MSEが過小・過大だった程度について調べてみよう。
dfPlot <- b.dfSimulation1 %>%
filter(!is.na(gEst_rstan)) %>%
dplyr::select(nSize, gV, starts_with("gEst_"), starts_with("gMSE_")) %>%
gather(sVar, gValue, c(starts_with("gEst_"), starts_with("gMSE_"))) %>%
separate(sVar, c("sVar", "sEstimator")) %>%
spread(sVar, gValue) %>%
filter(sEstimator != "nlme") %>%
mutate(
gDiff = gMSE - (gEst - gV)^2
) %>%
group_by(sEstimator, nSize) %>%
summarize( gMeanDiff = mean(gDiff) ) %>%
ungroup() %>%
mutate( fEstimator = b.sub_EstimatorAsFactor.1(sEstimator) )
g <- ggplot(data=dfPlot, aes(x = nSize, y=gMeanDiff, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + scale_x_continuous(breaks=seq(10, 100, by = 10))
g <- g + geom_hline(yintercept = 0)
g <- g + labs(
x = "Sample Size",
y = "mean( (Estimated MSE)-(Actual Sq.Err) )",
color = "Estimator"
)
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
縦軸は、「推定されたMSEから、推定値と正解との差の二乗を引いた値」の平均である。 正の値はMSEの推定が概して過大であること、負の値は概して過小であることを表す。
全体的なでこぼこは運の問題として、 ここで注目したいのは推定方法間の差である。よくみると、標本サイズが10のとき、 EBLUP推定のMSEの推定が過小であることが見て取れる。 とくにhbsaeパッケージのEBLUP推定は過小となる程度が大きく、標本サイズが20, 30のときにも やや過小になっている。
EBLUP推定のMSEの推定が過小になるのは、 EBLUP推定では(1)いったん\(\hat{\sigma^2_v}\)の点推定値を求め、 (2)これを使って\(\theta_i\)を推定する、という方法を取っているからだと思う。 つまり、本当は\(\hat{\sigma^2_v}\)にも分散があるのだが、それを無視しているわけで、 その分MSEが過小評価されるのだ… という理解で、あってますかね?
さきほどは地域の数を10に固定していた。これが多くなったらどうなるか?
さきほどと同じく、4ブロック、地域効果の分散は50とし、 地域の数を10, 100, 200, 300と動かしながら、 500回づつ繰り返してみよう。
なお、このシミュレーションはすごく時間がかかるので、nlmeパッケージと rstanパッケージは省略した。
# set.seed(123)
# lOut <- lapply(
# c(10, 100, 200, 300),
# function(nNumArea){
# cat("nNumArea = ", nNumArea, "...\n")
# dfData <- b.sub_GenerateMultiDataset1.1(
# nNumDataset = 500, nNumArea = nNumArea, nNumGroup = 4L, gVarV = 50
# )
# dfResult <- b.sub_EstimateMultiDataset1.1(
# dfData,
# asMethod = c("saeREML", "hbsaeREML", "hbsaeHB", "bayessae")
# )
# dfResult$nNumArea <- nNumArea
# cat("\n")
# return(dfResult)
# }
# )
# b.dfSimulation2 <- bind_rows(lOut) %>%
# filter(nArea == 1)
# save(b.dfSimulation2, file="./b_dfSimulation2.RData")
load("./b_dfSimulation2.RData")
なお、以下では標本サイズが10である地域だけに注目する。 また、条件間で分析件数を揃えるために、地域番号1だけに限定して分析する。 各条件において500件の推定結果が使えるわけである。
地域数 x 推定方法別に、推定値と正解とのずれを集計してみよう。
dfPlot <- b.dfSimulation2 %>%
dplyr::select(nNumArea, gV, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(sEstimator, nNumArea) %>%
summarize(
nFreq = n(),
gRMSE = sqrt(mean((gEst - gV)^2))
) %>%
ungroup() %>%
mutate(
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator) # ,
# fNumArea = factor(nNumArea, levels=c(10, 100, 200, 300))
)
g <- ggplot(data=dfPlot, aes(x = nNumArea, y=gRMSE, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + scale_x_continuous(breaks = c(0, 10, 100, 200, 300))
g <- g + labs(x = "Num. Area", y = "observed RMSE", color="Estimator")
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
さきほど同様、縦軸は「推定値と正解との差の二乗の平均の平方根」(RMSE)を示している。 さきほど確認したように、地域数10の場合には、だいたい3点くらいずれる。 そのずれは、HB推定, EBLUP推定, 標本平均の順に小さい。
さて、上の図からは、地域数が増えるとHB推定とEBLUP推定におけるずれがさらに小さくなる ことがわかる。また、 HB推定とEBLUP推定の間の差は消失し、標本平均との差が広がる。
なるほど、小地域モデルによる推定は、いわば他の地域から 情報を借りてくるようなものだから、他の地域の数が大きいほうがより有利になるのだろう。
こんどは、地域数 x 推定方法別に、MSEの推定が 過小・過大だった程度について調べてみよう。
dfPlot <- b.dfSimulation2 %>%
dplyr::select(nNumArea, gV, starts_with("gEst_"), starts_with("gMSE_")) %>%
gather(sVar, gValue, c(starts_with("gEst_"), starts_with("gMSE_"))) %>%
separate(sVar, c("sVar", "sEstimator")) %>%
spread(sVar, gValue) %>%
mutate(
gTrueSE = (gEst - gV)^2,
gDiff = gMSE - gTrueSE
) %>%
group_by(sEstimator, nNumArea) %>%
summarize(
gMeanDiff = mean(gDiff)
) %>%
ungroup() %>%
mutate(
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator) #,
# fNumArea = factor(nNumArea, levels=c(10, 100, 200, 300))
)
g <- ggplot(
data=dfPlot,
aes(x = nNumArea, y=gMeanDiff, group=fEstimator, color=fEstimator)
)
g <- g + geom_path(size = 1)
g <- g + geom_hline(yintercept = 0)
g <- g + scale_x_continuous(breaks = c(0, 10, 100, 200, 300))
g <- g + labs(
x = "Num. Area",
y = "mean ((Estimated MSE) - (Actual Sq.Err))",
color = "Estimator"
)
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
例によって、縦軸は「推定されたMSEから推定値と正解との差の二乗を引いた値」の平均である。 正の値はMSEの推定が概して過大であること、負の値は概して過小であることを表す。
すでに、地域数10で標本サイズが10のとき、EBLUP推定(特にhbsaeパッケージ)で MSEが過小評価されるということがわかった。 しかし上の図からは、それは地域数が少ないときの話であって、 地域数が多ければ気にしなくてよいということがわかる。
ここでちょっと変な実験を試してみる。
ここまでは、地域レベルの共変量として、地域のブロック(4水準のカテゴリ変数)を投入していた。 いいかえると、3つの二値変数を共変量として投入していた。
実は、地域のブロックは幸福度と無関係である。つまり、 真の係数\(\mathbf{\beta}\)は\(\mathbf{0}\)である。でも分析者は そのことに気づいていないわけである。
この共変量の数が減ったら、また増えたら、何が起きるだろうか。
変なことに関心を持つと思われるかもしれないけれど、個人的には ちょっと切実な疑問なんです。だって、共変量というものは、往々にして 意味があるのかないのかわからないものじゃないですか。
というわけで、意味のない共変量をうっかり投入したときの弊害というのはどんなものなのか、 それは推定量によって違うのか、ひとつ調べてみたいものだと思ったわけです。
地域数は10とし、ブロックの数を2から5まで動かしながら、 500回繰り返す。 nlmeパッケージとrstanパッケージは省略した。
# set.seed(123)
# lOut <- lapply(
# c(2, 3, 4, 5),
# function(nNumGroup){
# cat("nNumGroup = ", nNumGroup, "...\n")
# dfData <- b.sub_GenerateMultiDataset1.1(
# nNumDataset = 500, nNumArea = 10, nNumGroup = nNumGroup, gVarV = 50
# )
# dfResult <- b.sub_EstimateMultiDataset1.1(
# dfData,
# asMethod = c("saeREML", "hbsaeREML", "hbsaeHB", "bayessae")
# )
# dfResult$nNumGroup <- nNumGroup
# cat("\n")
# return(dfResult)
# }
# )
# lOut <- lapply(
# lOut,
# function(dfIn){
# dfIn$fGroup <- as.integer(dfIn$fGroup)
# return(dfIn)
# }
# )
# b.dfSimulation3 <- bind_rows(lOut) %>%
# filter(nSize == 10)
# save(b.dfSimulation3, file="./b_dfSimulation3.RData")
load("./b_dfSimulation3.RData")
ここでも、標本サイズが10である地域だけに注目する。 各条件において500件の推定結果が使える。
ブロックの数 x 推定方法別に、推定値と正解とのずれを集計してみよう。
dfPlot <- b.dfSimulation3 %>%
dplyr::select(nNumGroup, gV, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(sEstimator, nNumGroup) %>%
summarize(
nFreq = n(),
gRMSE = sqrt(mean((gEst - gV)^2))
) %>%
ungroup() %>%
mutate(
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator),
fNumGroup = factor(nNumGroup, levels=2:5)
)
g <- ggplot(
data=dfPlot,
aes(x = fNumGroup, y=gRMSE, group=fEstimator, color=fEstimator)
)
g <- g + geom_path(size=1)
g <- g + labs(
x = "Num. Blocks",
y = "observed RMSE",
color = "Estimator"
)
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
ブロックの数が少ない、つまり共変量が少ないほうが、 推定値は真値に近づく。えっ、こんなに違うの…
ブロックの数 x 推定方法別に、MSEの推定が過小・過大だった程度について調べてみる。
dfPlot <- b.dfSimulation3 %>%
dplyr::select(nNumGroup, gV, starts_with("gEst_"), starts_with("gMSE_")) %>%
gather(sVar, gValue, c(starts_with("gEst_"), starts_with("gMSE_"))) %>%
separate(sVar, c("sVar", "sEstimator")) %>%
spread(sVar, gValue) %>%
mutate(
gTrueSE = (gEst - gV)^2,
gDiff = gMSE - gTrueSE
) %>%
group_by(sEstimator, nNumGroup) %>%
summarize(
gMeanDiff = mean(gDiff)
) %>%
ungroup() %>%
mutate(
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator),
fNumGroup = factor(nNumGroup, levels=2:5)
)
g <- ggplot(data=dfPlot, aes(x = fNumGroup, y=gMeanDiff, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + geom_hline(yintercept = 0)
g <- g + labs(
x = "Num. Blocks",
y = "mean ((Estimated MSE) - (Actual Sq.Err))",
color = "Estimator"
)
g <- g + scale_color_manual(values = b.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
ブロック数5のときにちょっと高くなっている理由がよくわかんないんだけど、 ブロック数を2から3に増やすと、MSEが過小評価されていることがわかる。 無駄な共変量は、MSEの過小評価にもつながるようだ…
最後に、地域数と計算にかかる時間との関連について調べてみる。
100地域, 1,000地域, 10,000地域, 100,000地域の5通りについて、 5回ずつためしてみた。
ここでは、以下の推定方法を試す。
なお、
# set.seed(12345)
# lOut <- lapply(
# c(100,1000,10000,100000),
# function(nNumArea){
# cat("nNumArea = ", nNumArea, "...\n")
#
# cat("generate data ...\n")
# dfData <- b.sub_GenerateMultiDataset1.1(
# nNumDataset = 5, nNumArea = nNumArea, nNumGroup = 4L, gVarV = 50
# )
#
# asMethod <- c(
# "saeREML", "saeREMLNOMSE", "hbsaeREML", "hbsaeHB", "bayessae", "nlme"
# )
# ## sae takes too long with nNumArea >= 10000
# if (nNumArea >= 10000){
# asMethod <- asMethod[!(asMethod %in% c("saeREML", "saeREMLNOMSE"))]
# }
# ## BayesSAE raises out-of-memory error with nNumArea >= 100000
# if (nNumArea >= 100000){
# asMethod <- asMethod[asMethod != "bayessae"]
# }
#
# dfResult <- b.sub_EstimateMultiDataset1.1(
# dfData, asMethod = asMethod, bVerbose = TRUE
# ) %>%
# filter(nArea == 1) %>%
# mutate(nNumArea = nNumArea)
# cat("\n")
# return(dfResult)
# }
# )
# b.dfSimulation4 <- bind_rows(lOut)
# save(b.dfSimulation4, file="./b_dfSimulation4.RData")
load(file="./b_dfSimulation4.RData")
一回あたりの計算時間は下図となった。
dfPlot1 <- b.dfSimulation4 %>%
dplyr::select(nNumArea, nDataset, starts_with("gTime_"), starts_with("bStatus_")) %>%
gather(sVar, gValue, c(starts_with("gTime_"), starts_with("bStatus_"))) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
spread(sVar1, gValue) %>%
filter(bStatus == 1) %>%
mutate(
fEstimator = b.sub_EstimatorAsFactor.1(sEstimator),
fNumArea = factor(
nNumArea,
levels = c(100,1000,10000,100000),
labels = c("100","1000","10000","100000")
)
)
dfPlot2 <- dfPlot1 %>%
group_by(fNumArea, fNumArea, fEstimator) %>%
summarize(gTime = mean(gTime)) %>%
ungroup()
g <- ggplot(data=dfPlot2, aes(x = fNumArea, y=gTime, fill=fEstimator))
g <- g + geom_col()
g <- g + geom_point(data=dfPlot1, aes(x = fNumArea, y = gTime), size=1, alpha=1/10)
g <- g + labs(y = "Mean time (sec)", x = "Num.Area", fill="Estimator")
g <- g + scale_fill_manual(guide=FALSE, values = b.sub_EstimatorPalette.1(dfPlot1$fEstimator))
g <- g + facet_grid(~ fEstimator)
g <- g + theme_bw()
g <- g + theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust = 0.3))
print(g)
MCMCで推定しているBayesSAEパッケージが遅くなるのは わかるんだけど、不思議なのはsaeパッケージである。 MSEを算出してもしなくても、 地域数がある程度大きくなると急激に時間を食うようになる。
地域数が10万を超えるようなケースでは、hbsaeパッケージでEBLUP推定するか、 MSEは出せないけどnlmeパッケージのlme()を使うのが速そうだ。 lme()の推定精度が悪かった(3.3.4)ことを考えると、 hbsaeパッケージによるEBLUP推定が最善手であろう。
ここだけの話、このシミュレーションが今回の最大の発見でありました。 これまでsaeパッケージだけで悪戦苦闘していたのはなんだったのか…
このシミュレーションでわかったことをまとめておく。
rm(list=grep("^b\\.", ls(), value=T))
(この項はRM本4.4.4節, 8.4節に基づく。)
おさらいすると、Fay-Herriotモデルは以下の通りであった。 \[ \hat{\theta}_i = \mathbf{z}^T_i \mathbf{\beta} + v_i + e_i, \ \ v_i \sim (0, \sigma^2_v), \ \ e_i \sim (0, \psi_i) \] ところで、上の\(v_i\)の式はこう書き換えることができる。 \(\mathbf{v} = (v_1, \ldots, v_m)^T\)として、 \[ \mathbf{v} \sim (0, \mathbf{\Gamma}), \ \ \mathbf{\Gamma} = \sigma_v^2 \mathbf{I} \]
\(\mathbf{\Gamma}\)は\(m \times m\)の行列で、\(i\)行\(j\)列目の要素が、 地域\(i\)と\(j\)の共分散を表す。Fay-Herriotモデルの 場合は、すべての対角成分に\(\sigma_v^2\) が入っており、すべての非対角成分に\(0\)がはいる(地域効果\(v_i\)は地域間で 独立だと考えているから)。
このように、Fay-Herriotモデルでは、 地域効果\(v_i\)は地域間で独立だと考えていた。 でも実際には、 地域の間には隣接関係があり、隣り合う地域のあいだでは地域効果が似ている ことが多い。こうした空間相関を考慮すれば、 地域特性をもっと正確に推測できるかもしれない。
そこで、\(\mathbf{\Gamma}\)の非対角成分が0でないと考える。 地域効果\(\mathbf{v}\)が、地域間の空間的な隣接関係とどのように関連しているかを モデル化し、このモデルに基づいて、\(\mathbf{\Gamma}\)がどういう形式になっている かを規定する。これが空間Fay-Herriotモデルである。
推定の方法には、Fay-Herriotモデルと同様、EBLUP推定, HB推定などがある。
では、共分散行列\(\mathbf{\Gamma}\)をどのように規定するか。 代表的なアプローチが2つある。
いきなり\(\mathbf{v}\)の同時分布をモデル化するアプローチ。
次のように考える。 \[ \mathbf{v} = \phi \mathbf{W} \mathbf{v} + \mathbf{u}, \ \ \mathbf{u} \sim N(\mathbf{0}, \sigma_u^2 \mathbf{I}) \] つまりこういうことだ。オリジナルのFay-Herriotモデルでは、 地域\(i\)の地域効果\(v\)は、平均0, 分散\(\sigma_v^2\)のなんらかの分布に従う ランダムな効果だと想定していた。 いっぽうこのモデルでは、地域\(i\)の地域効果\(v_i\)は次のふたつの要素の 和だと想定する。
この\(w_{ij}\)を行列の形に並べたのが\(\mathbf{W}\)である。この\(\mathbf{W}\)は 地域間の隣接関係を表現している行列で、分析者が指定する。
以上の想定はこういう風に書き換えることができる。 \[ \mathbf{v} = (\mathbf{I}-\phi \mathbf{W})^{-1} \mathbf{u}, \ \ \mathbf{u} \sim N(\mathbf{0}, \sigma_u^2 \mathbf{I}) \] ここから、\(\mathbf{v}\)の分布は \[ \mathbf{v} \sim N(0, \mathbf{\Gamma}), \ \ \mathbf{\Gamma} = \sigma_v^2[(\mathbf{I}-\phi \mathbf{W})(\mathbf{I}-\phi \mathbf{W})^T]^{-1} \] となる。空間相関を考慮しないFay-Herriotモデルでは\(\sigma_v^2\)を推定する必要があったが、 ここではそのかわりに\(\sigma_u^2\)と\(\phi\)を推定する必要がある。
ここで鍵になるのは、隣接行列\(\mathbf{W}\)の定義である。 通常は次のように定義する。
このとき、\(\phi \in (0,1)\)となる。このときの\(\phi\)を 空間自己相関パラメータという。
\(\mathbf{v}\)の同時分布ではなく、 他のエリア\(\{v_l: l \neq i\}\)の下での、\(v_i\)の条件つき分布を モデル化する、というアプローチである。
各地域\(i\)は既知の定数\(b_i\)を持っているとする(典型的には1とする)。 次のように考える。地域\(i\)の「近接」地域の集合を\(A_i\)として \[ b_i v_i | \{v_l: l \neq i\} \sim N \left( \rho \sum_{l \in A_l} q_{il} b_l v_l, b_i^2 \sigma_v^2 \right) \] ここで\(\{q_{il}\}\)は既知の定数で、\(q_{il}b_l^2 = q_{li} b_i^2\)を満たす。 典型的には、\(l \in A_i\)のとき\(q_{il} = q_{li} = 1\)とする。
つまりこういうことだ。オリジナルのFay-Herriotモデルでは、 地域\(i\)の地域効果\(v_i\)は、平均0, 分散\(\sigma_v^2\)のなんらかの分布に従う ランダムな効果だと想定していた。 いっぽうこのモデルでは、なぜだか知らないが、\(v_i\)は隣接地域の 地域効果が固定されたときに、次のように正規分布する、と想定する。
ここから、\(\mathbf{\Gamma}\)を規定する次式が得られる。 \(b_1^2, \ldots, b_m^2\)を持つ対角行列を\(\mathbf{B}\)とし、 すべての\(q_{il}\)を持つ\(m \times m\)行列(近接じゃないところと対角は0にする)を \(\mathbf{Q}\)として \[ \mathbf{B}^{1/2} \mathbf{v} \sim N_m(\mathbf{0}, \mathbf{\Gamma}) \] \[ \Gamma = \sigma^2_v(\mathbf{I} - \rho \mathbf{Q})^{-1} \mathbf{B} \] 空間相関を考慮しないFay-Herriotモデルでは\(\sigma_v^2\)を推定する必要があったが、 ここでは\(\sigma_v^2\)と\(\rho\)を推定する必要がある。
前項まででは、空間Fay-Herriotモデルにおいて空間相関を規定するための2つの モデル(SARモデルとCARモデル)を取り上げた。
なんだか難しいので、すごく単純な具体例で考えてみよう。
いま地域が4つあり、地域1と2, 2と3, 3と4が隣接しているとする。 SARモデルとCARモデルでは、空間相関はどのように規定されるか。
まず隣接行列について。 SARモデルの場合、隣接行列\(\mathbf{W}\)は次のように定義されることが多い。 各行の和が1になっている点に注意。
e.mgProx_unstd <- matrix(0, nrow=4, ncol=4)
e.mgProx_unstd[1,2] <- 1
e.mgProx_unstd[2,3] <- 1
e.mgProx_unstd[3,4] <- 1
e.mgProx_unstd[2,1] <- 1
e.mgProx_unstd[3,2] <- 1
e.mgProx_unstd[4,3] <- 1
e.mgProx_std <- sweep(e.mgProx_unstd, 1, rowSums(e.mgProx_unstd), "/")
print(e.mgProx_std)
[,1] [,2] [,3] [,4]
[1,] 0.0 1.0 0.0 0.0
[2,] 0.5 0.0 0.5 0.0
[3,] 0.0 0.5 0.0 0.5
[4,] 0.0 0.0 1.0 0.0
いっぽうCARモデルの場合、隣接行列\(\mathbf{Q}\)は次のように定義されることが多い。 隣接する地域間のみ1, ほかを0として
print(e.mgProx_unstd)
[,1] [,2] [,3] [,4]
[1,] 0 1 0 0
[2,] 1 0 1 0
[3,] 0 1 0 1
[4,] 0 0 1 0
次に地域効果の定義について。
SARモデルの場合は下式のように定義される。
\[ v_1 = \phi v_2 + u_1, \ \ u_1 \sim N(0, \sigma_u^2) \] \[ v_2 = \phi (0.5 v_1 + 0.5v_3) + u_2, \ \ u_2 \sim N(0, \sigma_u^2) \] \[ v_3 = \phi (0.5v_2 + 0.5v_4) + u_3, \ \ u_3 \sim N(0, \sigma_u^2) \] \[ v_4 = \phi v_3 + u_4, \ \ u_4 \sim N(0, \sigma_u^2) \]
ここから、\(\mathbf{v}\)の分布は \[ \mathbf{v} \sim N(0, \mathbf{\Gamma}), \ \ \mathbf{\Gamma} = \sigma_v^2[(\mathbf{I}-\phi \mathbf{W})(\mathbf{I}-\phi \mathbf{W})^T]^{-1} \]
となる。
いっぽうCARモデルの場合は、すべての\(b_i\)を1とすれば
\[ v_1 | v_2 \sim N(\rho v_2, \sigma_v^2) \] \[ v_2 | \{ v_1, v_3\} \sim N(\rho (v_1+v_3), \sigma_v^2) \] \[ v_3 | \{ v_2, v_4\} \sim N(\rho (v_2+v_4), \sigma_v^2) \] \[ v_4 | v_3 \sim N(\rho v_3, \sigma_v^2) \]
これを同時分布に書き直すと \[ \mathbf{v} \sim N(\mathbf{0}, \mathbf{\Gamma}), \ \ \mathbf{\Gamma} = \sigma^2_v(\mathbf{I} - \rho \mathbf{Q})^{-1} \]
となる。
では、共分散行列\(\mathbf{\Gamma}\)をくらべてみよう。 SARモデルの場合、\(\sigma_u^2=1, \phi=0.5\)のとき
mbUnit <- diag(rep(1,4))
e.mgCov2 <- solve( (mbUnit - 0.5 * e.mgProx_std) %*% t(mbUnit - 0.5 * e.mgProx_std) )
print(e.mgCov2)
[,1] [,2] [,3] [,4]
[1,] 1.4419753 1.145679 0.454321 0.1580247
[2,] 1.1456790 2.093827 1.106173 0.4543210
[3,] 0.4543210 1.106173 2.093827 1.1456790
[4,] 0.1580247 0.454321 1.145679 1.4419753
CARモデルの場合、\(\sigma_v^2=1, \rho=0.5\)のとき
e.mgCov1 <- solve(diag(rep(1,4)) - 0.5 * e.mgProx_unstd)
print(e.mgCov1)
[,1] [,2] [,3] [,4]
[1,] 1.6 1.2 0.8 0.4
[2,] 1.2 2.4 1.6 0.8
[3,] 0.8 1.6 2.4 1.2
[4,] 0.4 0.8 1.2 1.6
となる。
相関行列に直してみると, SARモデルは
mgCor <- sweep(e.mgCov2, 1, sqrt(diag(e.mgCov2)), "/")
mgCor <- sweep(mgCor, 2, sqrt(diag(e.mgCov2)), "/")
print(mgCor)
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.6593464 0.2614649 0.1095890
[2,] 0.6593464 1.0000000 0.5283019 0.2614649
[3,] 0.2614649 0.5283019 1.0000000 0.6593464
[4,] 0.1095890 0.2614649 0.6593464 1.0000000
CARモデルは
mgCor <- sweep(e.mgCov1, 1, sqrt(diag(e.mgCov1)), "/")
mgCor <- sweep(mgCor, 2, sqrt(diag(e.mgCov1)), "/")
print(mgCor)
[,1] [,2] [,3] [,4]
[1,] 1.0000000 0.6123724 0.4082483 0.2500000
[2,] 0.6123724 1.0000000 0.6666667 0.4082483
[3,] 0.4082483 0.6666667 1.0000000 0.6123724
[4,] 0.2500000 0.4082483 0.6123724 1.0000000
このように、SARモデルとCARモデルは、定式化は全く異なるものの、 空間相関についての規定は、結局それほど変わらない。 いずれのモデルでも、
rm(list=grep("^e\\.", ls(), value=T))
かんたんなデータについて、空間Fay-Herriotモデルを推定してみよう。
saeパッケージについているサンプルデータ grapes を使う。 イタリアはトスカーナ地方の224自治体についてのデータで、以下の変数を 持っている。
### c.dfGrapes
data(grapes)
c.dfGrapes <- grapes %>%
mutate(
nAreaID = seq_along(grapehect),
SE = sqrt(var)
)
rm(grapes)
### c.mgProx_std: row-standardized neighbor matrix
data(grapesprox)
c.mgProx_std <- as.matrix(grapesprox)
rm(grapesprox)
colnames(c.mgProx_std) <- as.character(1:ncol(c.mgProx_std))
rownames(c.mgProx_std) <- as.character(1:nrow(c.mgProx_std))
## print(c.mgProx_std)
### c.mbProx_unstd: unstandardized neighbor matrix
c.mbProx_unstd <- c.mgProx_std
c.mbProx_unstd[c.mbProx_unstd > 0] <- 1
## print(c.mbProx_unstd)
summary(c.dfGrapes %>% dplyr::select(-c(nAreaID, SE)))
grapehect area workdays
Min. : 0.2478 Min. : 35.81 Min. : 22.15
1st Qu.: 38.9887 1st Qu.: 205.24 1st Qu.:104.94
Median : 60.3268 Median : 443.42 Median :138.25
Mean : 69.5096 Mean : 639.16 Mean :149.61
3rd Qu.: 85.7415 3rd Qu.: 830.28 3rd Qu.:178.59
Max. :342.9919 Max. :2890.13 Max. :497.68
var
Min. : 0.00
1st Qu.: 25.97
Median : 211.52
Mean : 2588.58
3rd Qu.: 1105.37
Max. :102745.66
データを散布図行列で観察してみよう。 なお、varはみにくいので平方根をとって標準誤差(SE)にする。
ggpairs(
c.dfGrapes %>% dplyr::select(grapehect, area, workdays, SE),
progress = FALSE,
aes(alpha = 1/10)
)
平均面積の直接推定量grapehectは、総面積areaと弱い相関があり、 謎の地域特性workdaysと比較的強い相関がある。 またSEとも正の相関がある(それはそうなるだろう。広い畑が多い自治体の ほうが、面積のばらつきも大きくなるだろうから)。
このデータに付随して、地域間隣接行列grapesproxも提供されている。 左上の8行8列だけ、ちらっとみてみると…
print(c.mgProx_std[1:8,1:8])
1 2 3 4 5 6 7
1 0.00 0.3333333 0.3333333 0.0000000 0.0000000 0.0000000 0.0000000
2 0.25 0.0000000 0.0000000 0.2500000 0.0000000 0.0000000 0.2500000
3 0.50 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
4 0.00 0.3333333 0.0000000 0.0000000 0.3333333 0.0000000 0.3333333
5 0.00 0.0000000 0.0000000 0.1428571 0.0000000 0.1428571 0.1428571
6 0.00 0.0000000 0.0000000 0.0000000 0.5000000 0.0000000 0.0000000
7 0.00 0.2000000 0.0000000 0.2000000 0.2000000 0.0000000 0.0000000
8 0.20 0.2000000 0.2000000 0.0000000 0.0000000 0.0000000 0.2000000
8
1 0.3333333
2 0.2500000
3 0.5000000
4 0.0000000
5 0.0000000
6 0.0000000
7 0.2000000
8 0.0000000
こんな感じ。各行の合計が1になっている。
隣接行列を眺めていても全然面白くない。地図をみたいんだけど、あいにく 提供されていない。仕方がないので、地域をノード、隣接をリンクにした グラフを描いてみた。
sub_ShowAreaGraph1.1 <- function(
mgProx, agSizeVar, agColorVar, sPal="OrRd"
){
g <- graph.adjacency(mgProx, mode="undirected")
V(g)$size <- agSizeVar
V(g)$color <- smoothPalette(agColorVar, pal=sPal)
set.seed(122)
par(mar=c(0,0,0,0))
plot(
g,
vertex.label.color = "black",
vertex.label.family = "mono",
vertex.label.cex = 0.5,
margin = -0.15
)
}
sub_ShowAreaGraph1.1(
c.mbProx_unstd, c.dfGrapes$area/120, c.dfGrapes$grapehect
)
円の大きさはブドウ畑総面積、 色の濃さはブドウ畑平均面積(色が濃い地域は平均面積が広い)。 なお、位置はソフトが適当に決めており、あまり意味がない。
グラフの左側に孤立した地域群があるけど、これは 島かなにかかなあ、とイタリアの地図を調べてみた。
グラフと見比べてみると、グラフ左側の孤立した地域群は、ポルト・フェライロという 島なんじゃないかと思うが、他の箇所については全く対応がわからない。うーむ、残念。 でもまあ、 色の濃い地域が隣接しているということがわかったのでよしとしよう。
地理には全く疎く、海外旅行もしないつまらない男なので、全然見当がつかないんですけど、 トスカーナ地方というのは、北はこの地図の左上のラ・スペツィアというあたりから、 南は地図の下側のオルベテッロというあたりまで、イタリア半島の西側の一帯を占める 州なのだそうだ。
イタリア政府観光局様によれば、 「糸杉が連なる美しい丘陵地帯は、世界遺産のヴァル・ドルチャをはじめ、 世界的に有名な大変美しいトスカーナの自然景観で、 街々を移動する間も旅行者の目を楽しませてくれ」るのだそうである。 州都はフィレンツェ。あ、それ知ってる!須賀敦子の本に出てくるよね!
写真を眺めていたら、人生どこで間違えたんだろうかと、 だんだん辛くなってきたので、話をもとに戻して…
さきほどのグラフによれば、ブドウ畑の平均面積は、どうやら隣接している 地域間で似ているようであった。 こういう性質を空間自己相関という。
空間自己相関の強さはどのくらいだろうか。 RM本の内容からは離れるが、いくつか指標を求めてみよう。
空間自己相関の指標には大きく分けて2種類ある。
その1、グローバルな指標。つまり、全ての地域を通じた空間自己相関の 指標。その代表選手が、MoranのI統計量である。 地域\(i\)の値を\(x_i\), 地域\(i\)と\(j\)のあいだの 隣接関係を表す重みを\(w_{ij}\)として、 \[ I = \frac{m}{\sum_i \sum_k w_{ij}} \frac{ \sum_i \sum_j w_{ij} (x_i-\bar{x})(x_j - \bar{x}) }{ \sum_i (x_i - \bar{x})^2 } \] この指標は、空間自己相関がないときに\(-1/(m-1)\)に近づく。
計算してみよう。\(w_{ij}\)を4.1で述べた重みとすると、
moran.test(c.dfGrapes$grapehect, mat2listw(c.mgProx_std, style="W"))
Moran I test under randomisation
data: c.dfGrapes$grapehect
weights: mat2listw(c.mgProx_std, style = "W")
Moran I statistic standard deviate = 6.2334, p-value = 2.282e-10
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.233579488 -0.003663004 0.001448550
I=0.23、弱い正の空間自己相関がある。
その2, ローカルな指標。隣接地域と似ているのはどの地域かを調べる ために使う。その代表選手はMoranのローカルI統計量で、 これはグローバルな\(I\)を各地域に分解する。つまり、 各地域について値が得られ、その合計はグローバルな\(I\)に一致する。 \[ I_i = \frac{m}{\sum_i \sum_k w_{ij}} \frac{ (x_i-\bar{x})\sum_j w_{ij} (x_j - \bar{x}) }{ \sum_i (x_i - \bar{x})^2 } \] 計算してみよう。
mgLocalMoran <- localmoran(
c.dfGrapes$grapehect,
mat2listw(c.mgProx_std, style="W")
)
sub_ShowAreaGraph1.1(
c.mbProx_unstd, c.dfGrapes$area/120, mgLocalMoran[,1], sPal = "Blues"
)
隣接地域との空間相関が高い地域を濃い青にしてみた。中央に、 青色が濃い地域が固まっている。さっきの図に戻ってみると、 これらはブドウ畑平均面積が広い地域である。
(本節はRM本8.11節の再現である。)
ブドウ畑平均面積をsaeパッケージで推定してみよう。
saeパッケージでは、SARモデルに基づく 空間Fay-Herriotモデルを推定できる。 地域特性はEBLUP推定。 \(\sigma_u^2\)と\(\phi\)の推定量はデフォルトではREMLとなる。
saeパッケージは空間Fay-Herriotモデルの関数として次の4つの関数を提供している。 いずれも地域特性の推定は同じで、MSEの扱いが異なる。
ここではmseSFH()を使うことにする。
まず、切片項を抜いたモデルを推定する。 ブドウ畑総面積areaとブドウ畑平均面積の散布図をみてみると、 回帰直線は原点を通りそうだから、という理由であろう。 RM本のモデルもそうなってますし(p.266), saeパッケージのmseSFH()のマニュアルのコード例でもそうなってます。
oModel <- mseSFH(
formula = grapehect ~ area + workdays - 1 ,
vardir = var,
proxmat = c.mgProx_std,
data = c.dfGrapes
)
cat(
"refvar:", oModel$'est'$fit$refvar,
"\nspatialcorr:", oModel$'est'$fit$spatialcorr, "\n"
)
refvar: 69.74899
spatialcorr: 0.6142697
\(\sigma_u^2\)は69.7, \(\phi\)は0.61と推定された。 … あれ? RM本p.266とまったく同じコードなのに、値がかなりちがう (RM本では71.189, 0.583)。なぜだ???
これは私のパソコンがおかしいわけではなくて、 FAO統計部の傘下にあるらしき組織による テクニカル・レポート での解説でも、私の出力と同じ値になっている(p.41)。うーむ…
とかなり思い悩んだんだけど、途中で気が付きました。切片項を抜くのをやめてみると
### c.oModel_saeREML
c.oModel_saeREML <- mseSFH(
formula = grapehect ~ area + workdays ,
vardir = var,
proxmat = c.mgProx_std,
data = c.dfGrapes
)
cat(
"refvar:", c.oModel_saeREML$'est'$fit$refvar,
"\nspatialcorr:", c.oModel_saeREML$'est'$fit$spatialcorr, "\n"
)
refvar: 71.1893
spatialcorr: 0.5826043
print(c.oModel_saeREML$'est'$fit$estcoef)
\(\sigma_u^2\)は71.189, \(\phi\)は0.583。 RM本の出力と一致しました。やれやれ。
このように、RM本はコード例と出力が矛盾しているわけだ。 では、どっちを採用するべきか。ここでは切片項を入れるほうを採用する。 「平均面積は総面積とworkdaysの線形和になる」というのは なんだか変なモデルだと思うので。
さて、グラフの色を、このモデルによる推定値に描き変えると…
sub_ShowAreaGraph1.1(c.mbProx_unstd, c.dfGrapes$area/120, c.oModel_saeREML$'est'$eblup)
色の濃さがちょっと変わった。たとえば、中心の色の濃いノード(地域173)の左上にある、 地域187の色がちょっと薄くなっている。
さて、地域特性の推定値は、空間相関を考慮するかどうかによっても変わるし、 どんな共変量を使うかによって変わる。次の5種類を比べてみよう。
### c.oModel_saeREML_{nox, area, areaworkday}
c.oModel_saeREML_nox <- mseFH(
formula = grapehect ~ 1,
vardir = var,
data = c.dfGrapes
)
c.oModel_saeREML_area <- mseFH(
formula = grapehect ~ area ,
vardir = var,
data = c.dfGrapes
)
c.oModel_saeREML_areaworkday <- mseFH(
formula = grapehect ~ area + workdays,
vardir = var,
data = c.dfGrapes
)
モデルの適合度を調べてみると…
lModel <- list(
c.oModel_saeREML_nox,
c.oModel_saeREML_area,
c.oModel_saeREML_areaworkday,
c.oModel_saeREML
)
dfTemp <- data.frame(
Estimator = c("none", "area", "area+workdays", "area+workdays+spatialcorr"),
AIC = sapply(lModel, function(oModel) oModel$est$fit$goodness[2]),
BIC = sapply(lModel, function(oModel) oModel$est$fit$goodness[3]),
stringsAsFactors = FALSE
)
print(dfTemp)
AIC, BICいずれも、空間Fay-Herriotモデルが最小となった。
推定結果は下図の通り。
c.dfResult_saeREML <- c.dfGrapes %>%
mutate(
gEst_1 = as.vector(c.oModel_saeREML_nox$est$eblup),
gEst_2 = as.vector(c.oModel_saeREML_area$est$eblup),
gEst_3 = as.vector(c.oModel_saeREML_areaworkday$est$eblup),
gEst_4 = as.vector(c.oModel_saeREML$est$eblup),
gMSE_3 = c.oModel_saeREML_areaworkday$mse,
gMSE_4 = c.oModel_saeREML$mse
) %>%
mutate(nArea = seq_along(grapehect))
dfPlot <- c.dfResult_saeREML %>%
arrange(-var) %>%
top_n(30, var) %>% # attention!!!!
mutate(
fArea = factor(seq_along(nArea), labels = as.character(nArea))
) %>%
rename(gEst_0 = grapehect) %>%
gather(sVar, gValue, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sVar2")) %>%
mutate(
fModel = factor(
as.integer(sVar2),
levels=0:4,
labels=c("direct", "none", "area", "area+workdays", "area+workdays+spatialcorr")
)
)
g <- ggplot(data=dfPlot, aes(x=fArea, y=gValue, group=fModel, color=fModel))
g <- g + geom_path(size=1)
g <- g + labs(x = "AreaID", y = expression(hat(theta)), color="Model")
g <- g + theme_bw()
print(g)
標本抽出分散が大きい地域(つまり、おそらくは標本サイズが小さかった地域) を30個選んで示す。 横軸は地域で、標本抽出分散の大きさで並べている。地域187は左から4番目にある。
推定値は、共変量としてworkdaysを投入したことによって 大きく変わり、空間相関を考慮してもしなくてもあまり変わっていない。 これはこの30地域に限ったことではなくて…
sub_PlotTwoEstimate.1(
agEst1 = as.vector(c.dfResult_saeREML$gEst_3),
agMSE1 = as.vector(c.dfResult_saeREML$gMSE_3),
agEst2 = as.vector(c.dfResult_saeREML$gEst_4),
agMSE2 = as.vector(c.dfResult_saeREML$gMSE_4),
asLabel = c.dfResult_saeREML$nArea,
sVarLabel1 = "area+workdays",
sVarLabel2 = "area+workdays+spatialcoor"
)
左側のパネルは各地域の推定値を表している。 横軸は空間相関を考慮していない場合、縦軸は考慮した場合。たいして変わっていない。
では、わざわざ空間Fay-Herriotモデルなんて使わなくてよかったんじゃないか? というと、これがそうでもない。右側は各地域の推定値のMSEの推定値を表している。 空間相関を考慮した場合のほうがMSEの推定値が下がっている。つまり、 もともとブドウ畑平均面積は隣接する自治体で似ているという傾向があるので、 それを考慮することで、推定値の精度を高めに見積もることができたわけだ。
こんどはBayesSAEパッケージ。
BayesSAEパッケージはCARモデルによる空間Fay-Herriotモデルを推定できる… はずである。推定量はHB。推定方法はMCMC。
では推定してみよう…と試みたところで、 ドツボにはまりました。
延々時間をかけた末、以下のことがわかった。
Rの再起動を繰り返しながら何度も試したなかで、運よく結果が得られたケースを 示す。mcmc=500。
# dfTemp <- expand.grid(
# nRow = 1:nrow(c.mbProx_unstd),
# nCol = 1:ncol(c.mbProx_unstd)
# )
# dfTemp$gValue <- c.mbProx_unstd[as.matrix(dfTemp)]
# mbTemp <- dfTemp %>%
# filter(gValue > 0, nRow < nCol) %>%
# dplyr::select(nRow, nCol) %>%
# as.matrix(.)
#
# set.seed(13)
# c.oModel_BayesSAE <- BayesSAE(
# grapehect ~ area + workdays | var,
# data = c.dfGrapes,
# spatial = TRUE,
# prox = mbTemp,
# mcmc = 500,
# burnin = 250
# )
# save(c.oModel_BayesSAE, file="./c_oModel_BayesSAE.RData")
load(file="./c_oModel_BayesSAE.RData")
summary(c.oModel_BayesSAE)
Call:
BayesSAE(formula = grapehect ~ area + workdays | var, spatial = TRUE,
prox = mbTemp, mcmc = 500, burnin = 250, data = c.dfGrapes)
CAR Area-Level Fay-Herriot Model
Sampling model: grapehect ~ theta
Linking model: theta ~ area + workdays with spatial random effect
Rao-Blackwellization of theta's based on the simulation:
1 2 3 4 5 6 7 8
31.1089 69.6696 73.8351 63.1305 38.9171 78.5383 49.6683 41.2826
9 10 11 12 13 14 15 16
111.1351 10.2331 79.9256 78.4572 60.6078 52.3993 59.1869 69.0589
17 18 19 20 21 22 23 24
103.4756 42.1818 34.9507 73.1143 71.6234 82.8210 60.2036 65.6269
25 26 27 28 29 30 31 32
69.5191 59.0232 76.5645 15.8110 2.0183 130.3570 34.4952 21.6494
33 34 35 36 37 38 39 40
15.0136 63.6233 60.2669 91.5382 54.9043 41.7187 82.4428 57.9240
41 42 43 44 45 46 47 48
0.6297 97.5886 77.6175 38.9341 16.9809 139.6731 48.8367 52.5079
49 50 51 52 53 54 55 56
29.1266 78.4874 65.3626 37.9817 59.7702 39.2380 19.6406 95.0978
57 58 59 60 61 62 63 64
74.3208 58.6386 57.5934 41.4696 96.3495 66.7592 80.4993 52.3589
65 66 67 68 69 70 71 72
38.8288 62.9432 45.9401 42.2684 43.9611 49.1812 77.7751 62.4326
73 74 75 76 77 78 79 80
35.1562 74.2700 53.7352 86.4640 101.7038 52.5921 44.5070 32.9756
81 82 83 84 85 86 87 88
144.3330 77.3188 88.7806 73.9564 52.1338 51.2183 24.3714 57.3625
89 90 91 92 93 94 95 96
36.2671 26.7954 29.0365 38.4890 162.4122 31.4012 48.4079 95.9136
97 98 99 100 101 102 103 104
46.4321 46.8414 68.3132 71.2937 70.9937 50.5978 106.7134 59.2189
105 106 107 108 109 110 111 112
37.9505 54.3702 29.9997 45.5565 48.4020 83.2776 29.5474 98.5622
113 114 115 116 117 118 119 120
94.9086 75.5132 80.4144 63.5484 63.1694 60.2385 51.2838 66.7885
121 122 123 124 125 126 127 128
39.0232 35.9737 72.3765 60.0282 52.7013 44.1830 26.9104 63.9960
129 130 131 132 133 134 135 136
115.4749 56.6792 118.7125 55.8009 61.2699 127.7067 38.8180 44.6875
137 138 139 140 141 142 143 144
55.9713 51.1405 63.5967 54.4470 151.9690 66.3482 70.7935 44.9629
145 146 147 148 149 150 151 152
67.6876 19.3465 73.1860 34.1474 81.6416 25.1924 54.2336 73.8958
153 154 155 156 157 158 159 160
71.9959 77.8823 39.0749 72.3258 33.4203 111.5003 43.4258 78.4384
161 162 163 164 165 166 167 168
58.2105 69.9291 89.3437 144.1329 141.5826 59.0028 32.9067 55.4793
169 170 171 172 173 174 175 176
74.1697 108.5298 76.2338 69.2812 218.3891 121.7793 19.9590 51.8617
177 178 179 180 181 182 183 184
62.9212 34.1808 33.3153 51.0139 52.6428 52.4443 88.4668 60.0182
185 186 187 188 189 190 191 192
72.9722 67.1534 168.9997 78.2465 87.9143 107.0522 86.4187 22.8224
193 194 195 196 197 198 199 200
28.5026 48.0452 68.5629 66.8568 18.0972 53.7315 84.2729 102.6104
201 202 203 204 205 206 207 208
36.2346 26.8451 52.4201 88.5696 55.4012 25.8781 50.3355 43.7183
209 210 211 212 213 214 215 216
97.2438 60.0605 76.1216 51.1254 68.3506 119.0295 51.7719 29.9704
217 218 219 220 221 222 223 224
37.8990 52.2895 96.8565 99.2800 49.6618 54.4060 39.7284 96.7855
225 226 227 228 229 230 231 232
75.5608 118.1376 57.9896 225.7792 81.1546 89.0848 91.3760 60.7030
233 234 235 236 237 238 239 240
47.3624 67.2553 42.3490 108.5470 55.3341 50.5188 85.9824 75.5469
241 242 243 244 245 246 247 248
49.0319 57.2949 34.9818 29.5782 32.9417 54.5701 32.2340 63.4536
249 250 251 252 253 254 255 256
16.0221 43.8043 40.6643 87.8827 46.4739 61.5904 114.4208 168.8226
257 258 259 260 261 262 263 264
58.3548 107.0767 89.7515 43.8791 93.2973 67.8464 80.2149 42.8581
265 266 267 268 269 270 271 272
53.6011 157.5936 119.6216 124.7967 111.2727 46.2415 52.2796 85.2828
273 274
88.1354 25.7205
Coefficients in the linking model:
Mean SD 2.5% 97.5%
beta.1 -3.947019 3.551080 -8.655680 2.854
beta.2 -0.011854 0.002403 -0.016726 -0.007
beta.3 0.516320 0.021010 0.478431 0.557
Variance of residual in the linking model:
Mean SD 2.5% 97.5
[1,] 404.3 654.7 156.7 2092
DIC: 2231.764
sub_PlotTwoEstimate.1(
agEst1 = as.vector(c.oModel_saeREML$est$eblup),
agMSE1 = c.oModel_saeREML$mse,
agEst2 = c.oModel_BayesSAE$HB,
agMSE2 = summary(MCMC(c.oModel_BayesSAE))$statistics[1:nrow(c.dfGrapes),2]^2,
asLabel = c.dfGrapes$nAreaID,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "BayesSAE"
)
… まあとにかく、パッケージを世に送り出す開発者の方には感謝しておりますが、 BayesSAEパッケージで空間Fay-Herriotモデルを推定するのは、現状かなり怖いな、 と思う次第である。
SARモデルについて、今度はrstanでやってみよう。 SARモデルはベイズ推定との相性があまりよくないようなのだが、 推定できないわけじゃないらしい。
まずstanコード。以下ではテキストとして示すが、 本物はファイルになっている。 正直、全然自信がないのですが…
// Spatial Fay-Herriot model (SAR)
// cf. https://rdrr.io/github/ssoro411/bayesae/src/R/Models_fast.R
data {
int<lower=1> M;
int<lower=0> P;
vector[M] Y;
vector<lower=0>[M] SD;
matrix<lower=0, upper=1>[M,M] W;
matrix[M,P] X;
}
transformed data {
matrix[M,M] unit;
unit = diag_matrix(rep_vector(1,M));
}
parameters {
real<lower=-0.99999, upper=0.99999> rho;
real<lower=0> sigma_sq;
vector[P] beta;
vector[M] theta;
}
model {
theta ~ multi_normal_prec(X*beta, (1/sigma_sq)*((unit-rho*(W))*(unit-rho*(W'))) );
// slower
// theta ~ multi_normal_prec( X*beta, (1/sigma_sq)*crossprod(unit-rho*(W)) ) ;
Y ~ normal(theta, SD);
}
Rコードは以下。ふだんはMplusの忠実な信者で、stanにあまり慣れてないもので、 これが速いのか遅いのかわからないのだが、PCを6時間くらい 唸らせる羽目になった。
# c.oStanModel <- stan_model(file="./SFH1.stan")
# save(c.oStanModel, file="./c.oStanModel.RData")
# load(file = "./c.oStanModel.RData")
# mgModelMatrix <- model.matrix(~ area + workdays, data=c.dfGrapes)
# lData <- list(
# M = nrow(c.dfGrapes),
# P = ncol(mgModelMatrix),
# Y = c.dfGrapes$grapehect,
# SD = c.dfGrapes$SE,
# W = c.mgProx_std,
# X = mgModelMatrix
# )
#
# c.oStanFit <- sampling(
# oStanModel_vec,
# data = lData,
# seed = 1234,
# iter = 500,
# verbose = TRUE,
# control = list(adapt_delta = 0.95),
# refresh = 10
# )
# save(c.oStanFit, file="./c.oStanFit.RData")
load(file="./c.oStanFit.RData")
stan_trace(c.oStanFit, c("rho", "sigma_sq", "theta[116]"), inc_warmup = T)
stan_rhat(c.oStanFit)
\(\rho, \sigma_u^2, \theta_{116},\)のトレースプロットと、 \(\hat{R}\)のヒストグラムを描いてみた。 なお、地域116は直接推定量の分散が最も高い。MCMCにおいてパラメータの\(\hat{R}\)が 一番高いのも\(\theta_{116}\)である(1.099)。
まあ… ぎりぎり収束しているということにしたいが、 warmup(現状ではデフォルトの250)がちょっと足りないような気がするので、 以下ではさらに50反復を捨て、最後の200反復だけ使うことにする。
地域特性の推定値をsaeパッケージと比較してみると…
mgResult <- rstan::extract(c.oStanFit, permuted = FALSE) # 3 dims
mgTheta <- mgResult[51:250, , grepl("^theta", dimnames(mgResult)$parameters)]
mgTheta <- do.call(rbind, lapply(1:4, function(nID) mgTheta[,nID,]))
agEst <- colMeans(mgTheta)
agMSE <- apply(mgTheta, 2, var)
sub_PlotTwoEstimate.1(
agEst1 = as.vector(c.oModel_saeREML$est$eblup),
agMSE1 = c.oModel_saeREML$mse,
agEst2 = agEst,
agMSE2 = agMSE,
asLabel = c.dfGrapes$nAreaID,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "rstan"
)
推定値はほぼ同じ、MSEは高めに推定されている。
saeパッケージで推定した空間Fay-Herriotモデルとは、要は 各地域の標本抽出誤差分散\(\psi_i\)が既知なSARモデルかCARモデルなのだから、 汎用的なパッケージで推定できるはずだ。
…と考え、SARモデルのほうについて、 果敢にもspdepパッケージでの推定にチャレンジしてみた。 正直、正しいコードなのか全く自信がない。
特に自信がないのが、標本抽出誤差分散\(\psi_i\)の指定の仕方である。 errorsarlm()にはweights=という引数がある。マニュアルにいわく、 「オプション。フィッティングの過程で用いられるウェイトのベクトル。 非NULLのウェイトを指定することで、観察によって分散が異なることを 指示することができる。この場合、ウェイトの値は分散と反比例する値とする」 とある。それで、きっとweightsには\(1/\psi_i\)を指定するんだろうと 思ったわけなのですが…
c.oModel_spdep <- errorsarlm(
grapehect ~ area + workdays ,
data = c.dfGrapes,
listw = mat2listw(c.mgProx_std, style="W"),
weights = 1/var
)
print(summary(c.oModel_spdep))
Call:
errorsarlm(formula = grapehect ~ area + workdays, data = c.dfGrapes,
listw = mat2listw(c.mgProx_std, style = "W"), weights = 1/var)
Residuals:
Min 1Q Median 3Q Max
-84.0203 -2.6305 11.4841 24.5922 225.4653
Type: error
Coefficients: (asymptotic standard errors)
Estimate Std. Error z value Pr(>|z|)
(Intercept) -25.3884753 1.0907410 -23.276 <2e-16
area 0.0015826 0.0015890 0.996 0.3192
workdays 0.4030366 0.0151736 26.562 <2e-16
Lambda: 0.58134, LR test value: 257.4, p-value: < 2.22e-16
Asymptotic standard error: 0.061051
z-value: 9.5223, p-value: < 2.22e-16
Wald statistic: 90.673, p-value: < 2.22e-16
Log likelihood: -1640.073 for error model
ML residual variance (sigma squared): 56.904, (sigma: 7.5435)
Number of observations: 274
Number of parameters estimated: 5
AIC: 3290.1, (AIC for lm: 3545.5)
\(\sigma_u^2\)は56.904(seaパッケージでは71.189), \(\phi\)は0.581(saeパッケージでは0.583)と推定された。
\(\sigma_u^2\)の推定値が全然ちがう。\(\mathbf{\beta}\)の推定値も かなり違う。暗雲垂れ込めてきたけれど、\(\theta_i\)の推定値を比べてみると…
sub_PlotTwoEstimate.1(
agEst1 = as.vector(c.oModel_saeREML$est$eblup),
agEst2 = as.vector(predict(c.oModel_spdep, pred.type="TS")),
asLabel = c.dfGrapes$nAreaID,
sVarLabel1 = "sae (REML)",
sVarLabel2 = "spdep"
)
案の定、かなり違う。なぜだ??
二つの可能性が考えられる。
あれこれ考えていて、実は後者なのかも、という気がしてきた。 その理由は次の通り。
空間Fay-Herriotモデルは、 \(\hat{\mathbf{\theta}} = (\theta_1, \ldots, \theta_m)^T\)とし、 \(\psi_1, \ldots, \psi_m\)を持つ対角行列を\(\mathbf{\Psi}\)と書くとして \[ \hat{\mathbf{\theta}} = \mathbf{Z} \mathbf{\beta} + \mathbf{v} + \mathbf{e}, \ \ \mathbf{e} \sim (0, \mathbf{\Psi}) \] \[ \mathbf{v} = \phi \mathbf{W} \mathbf{v} + \mathbf{u}, \ \ \mathbf{u} \sim N(\mathbf{0}, \sigma_u^2 \mathbf{I}) \] となると思う。いっぽうerrorsarlm()のモデルは、マニュアルによれば
y = X beta + u, u = lambda W u + e
である。記法を揃えると、これは \[ \hat{\mathbf{\theta}} = \mathbf{Z} \mathbf{\beta} + \mathbf{v} \] \[ \mathbf{v} = \phi \mathbf{W} \mathbf{v} + \mathbf{u}, \ \ \mathbf{u} \sim N(\mathbf{0}, \sigma_u^2 \mathbf{I}) \] ではないかと思うのである。weights=引数は 推定の際になんらか使われる重みづけに過ぎず、モデルには標本抽出誤差が含まれて いないのではなかろうか?
うーむ。残念ながら、よくわからない。 どなたか教えてくださらないでしょうか。
rm(list=grep("^c\\.", ls(), value=T))
こうしてみてみると、いろいろ疑問が湧いてくる。わたくしと致しましては、 特に次の3点に関心がある。
こんばんは、シミュレーションの時間です。
\(c\)行\(c\)列のグリッドで区切られた\(m=c*c\)個の地域があるとする。 いま、地域\(i\)(\(=1,\ldots,m\))のひとりひとりの住民\(j\)(\(=1,\ldots,N_i\))が、 量的な反応変数\(y_{ij}\)を持っている。えーっと、幸福度ってことにしましょうか。
地域の母集団サイズ\(N_i\)はいずれも十分に大きい。 各地域の幸福度の母平均\(\bar{Y}_i = (1/N_i) \sum_j y_{ij}\)について推測したい。
そこで、各地域\(i\)からサイズ\(n_i\)の標本を単純無作為抽出し、 標本平均\(\hat{\bar{Y}}_i = \frac{\sum_j y_{ij}}{n_i}\), ならびにその分散の推定値 \(\hat{Var}(\hat{\bar{Y}}_i) = \frac{\sum_j (y_{ij} - \hat{\bar{Y}}_i)}{n_i(n_i-1)}\) を求めた。
ところが、地域ごとの標本サイズ\(n_i\)は必ずしも十分でない。 また、辺を接する地域間では、幸福度平均\(\bar{Y}_i\)が似ている、 と私は信じている(なんかその、境界を超えて漂ってくるんじゃないですかね、 幸せのオーラみたいなものが)。 だから、幸福度の標本平均\(\hat{\bar{Y}}_i\)だけでなく、 その地域が他のどの地域と接しているかという情報もつかって、 各地域の幸福度の母平均\(\bar{Y}_i\)を 推測できるはずである。
そこで空間Fay-Erriotモデルの登場である。
以下、記号が煩雑になるので、 RM本にあわせて\(\theta_i = \bar{Y}_i\), \(\hat{\theta}_i = \hat{\bar{Y}}_i\) と書く。 また、 \(\psi_i = Var(\hat{\bar{Y}}_i)\)と書く。
おさらいすると、空間Fay-Herriotモデル(SARモデル)は以下の通り。 上の問題設定にあわせ、共変量を省いて書く。 \(\hat{\mathbf{\theta}} = (\theta_1, \ldots, \theta_m)^T\)とし、 \(\psi_1, \ldots, \psi_m\)を持つ対角行列を\(\mathbf{\Psi}\)と書く。
\[ \hat{\mathbf{\theta}} = \mathbf{v} + \mathbf{e}, \ \ \mathbf{e} \sim (\mathbf{0}, \mathbf{\Psi}) \] \[ \mathbf{v} = (\mathbf{I}-\phi \mathbf{W})^{-1} \mathbf{u}, \ \ \mathbf{u} \sim N(\mathbf{0}, \sigma_u^2 \mathbf{I}) \]
データ生成モデルは次のとおり。SARモデルに基づく。
地域 \(i=1,\ldots,m\) のユニット \(j=1,\ldots,n_j\) について、 \[ y_{ij} = 50 + v_i + \epsilon_{ij}, \ \ \epsilon_{ij} \mathop{\sim}^{iid} N(0,100) \] \(\mathbf{v} = (v_1, \ldots, v_m)^T\)について、 \[ \mathbf{v} = (\mathbf{I} - \phi \mathbf{W})^{-1} \mathbf{u}, \ \ \mathbf{u} \sim N(\mathbf{0}, \sigma_u^2 \mathbf{I}) \] \(\mathbf{W}\)は、辺を接する地域を1とした隣接行列を行標準化したものとする。
つまり、私は知らないのだが、実は
母集団サイズは大きく、\(\epsilon_{ij}\)の母集団平均は 0とみなせるので、結局、知りたい特性\(\theta_i\)はここでは\(v_i\)に等しい。
地域の標本サイズ\(n_i\)は、\(10, 20, \ldots, 100\)を長さ\(m\)に達するまで繰り返し、 これをランダムに並び替えたものとする。
…というデータを生成するための、一連の関数を作ってみました。
d.sub_MakeNeighborMatrix.1 <- function(nNumGridRow, nNumGridCol){
# purpose: make a neighbor matrix of meshes in specified grid
# args:
# nNumGridRow, nNumGridCol: size of grid
# return:
# a neighbor matrix (row-standardized)
# notes:
# e.g. nNumGridRow=2, nNumGridCol=3
# meshID in grid:
# 1, 2, 3
# 4, 5, 6
# neighbor matrix (unstandardized):
# (meshID=1) 0, 1, 0, 1, 0, 0
# (meshID=2) 1, 0, 1, 0, 1, 0
# (meshID=3) 0, 1, 0, 0, 0, 1
# (meshID=4) 1, 0, 0, 0, 1, 0
# (meshID=5) 0, 1, 0, 1, 0, 1
# (meshID=6) 0, 0, 1, 0, 1, 0
dfValue <- expand.grid(
# nGridRow1, nGridCol1: location of mesh1 in the grid
# nGridRow2, nGridCol2: location of mesh2 in the grid
nGridRow1 = 1:nNumGridRow,
nGridCol1 = 1:nNumGridCol,
nGridRow2 = 1:nNumGridRow,
nGridCol2 = 1:nNumGridCol
) %>%
mutate(
# nLoc1: meshID of mesh1
# nLoc2: meshID of mesh2
nLoc1 = (nGridRow1-1)*nNumGridCol + nGridCol1,
nLoc2 = (nGridRow2-1)*nNumGridCol + nGridCol2,
# bValue: neighbor or not
bValue = if_else(
(
(nGridRow1 == nGridRow2)
& (abs(nGridCol1 - nGridCol2) == 1)
)|(
(abs(nGridRow1 - nGridRow2)==1)
& (nGridCol1 == nGridCol2)
),
1, 0
)
)
dfName <- dfValue %>%
dplyr::select(nGridRow1, nGridCol1, nLoc1) %>%
distinct() %>%
mutate(sName = paste0("r", nGridRow1, "c", nGridCol1)) %>%
arrange(nLoc1)
## print(dfName)
out <- matrix(0, nrow=nNumGridRow*nNumGridCol, ncol=nNumGridRow*nNumGridCol)
out[as.matrix(dfValue[c("nLoc1", "nLoc2")])] <- dfValue$bValue
out <- sweep(out, MARGIN = 1, FUN = "/", STATS = rowSums(out))
rownames(out) <- dfName$sName
colnames(out) <- dfName$sName
# print(out)
return(out)
}
d.sub_GenerateDataset2.1 <- function(
mgNB, gVarU, gPhi, gVarE=100, anPossibleN=seq(10, 100, by = 10)
){
## purpose: generate a dataset for spatial FH model
## arg:
## mgNB: neighbor matrix
## gVarU: variance of area effects
## gPhi: spatial correlation
## gVarE: variance in each area
## anPossibleN: vector of possible sample size
## return:
## a data frame, area level data
## $nArea: Area ID
## $sName: Grid Location (e.g. "r2c3" means row2-col3)
## $nSize: sample size
## $gU: area effect without spatial correlation
## $gV: area effect with spatial correlation
## $gTrue: true value
## $gEst_Direct: sample means of a response variable
## $gMSE_Direct: sample est. of variance of sample means
nNumArea <- nrow(mgNB)
# anSize: sample size
anSize <- rep(anPossibleN, nNumArea %/% length(anPossibleN) + 1)[1:nNumArea]
anSize <- sample(anSize, nNumArea, replace=FALSE)
## print(anSize)
# agU: \sim N(\mathbf{0}, \sigma_u^2 \mathbf{I})
agU <- rnorm(n = nNumArea, mean = 0, sd = sqrt(gVarU))
# agV: (\mathbf{I} - \phi \mathbf{W})^{-1} \mathbf{u}
mgCoef <- solve(diag(x=1, nrow=nNumArea) - gPhi * mgNB)
agV <- as.vector(mgCoef %*% agU)
asName <- rownames(mgNB)
out <- data.frame(
nArea = 1:nNumArea,
sName = asName,
nGridRow = as.integer(sub("^r([0-9]+)c([0-9]+)$", "\\1", asName)),
nGridCol = as.integer(sub("^r([0-9]+)c([0-9]+)$", "\\2", asName)),
nSize = anSize,
gU = agU,
gV = agV,
gTrue = 50 + agV,
stringsAsFactors = FALSE
)
lError <- lapply(
out$nSize,
function(nSize){
rnorm(n = nSize, mean = 0, sd = sqrt(gVarE))
}
)
out$gEst_Direct <- out$gTrue + sapply(lError, mean)
out$gMSE_Direct <- sapply(lError, function(x){ var(x) / length(x) })
## print(out)
return(out)
}
d.sub_GenerateMultiDataset2.1 <- function(nNumDataset=1, mgNB, gVarU, gPhi, ...){
## purpose: generate multiple datasets for FH model
## arg:
## nNumDataset: number of data sets
## mgNB: neighbor matrix
## gVarU: variance of area effects
## gPhi: spatial correlation,
## ...: pass to sub_GenerateDataset2.1() (gVarE, anPossibleN)
##
## return:
## a dataset
## $nDataset: Dataset ID
## $... generated by sub_GenerateDataset2.1()
out <- data.frame(nDataset = 1:nNumDataset) %>%
group_by(nDataset) %>%
dplyr::do(d.sub_GenerateDataset2.1(mgNB = mgNB, gVarU = gVarU, gPhi = gPhi)) %>%
ungroup()
return(out)
}
d.sub_EstimateDatasetByMethod2.1 <- function(
dfData, mgNB, sMethod, bVerbose=FALSE
){
## purpose: estimate spatial FH model for a dataset by a method
## arg:
## dfData: a data set
## mgNB: neighbor matrix
## sMethod: {"saeFH", "saeSFH", "spdepWGT", "spdepNOWGT"}
## bVerbose
## return:
## a list
## trap: check sMethod
stopifnot(sMethod %in% c("saeFH", "saeSFH", "spdepWGT", "spdepNOWGT"))
if (bVerbose) cat("[sub_EstimateDataset2.1]", sMethod, "... \n")
t <- proc.time()
if (sMethod == "saeFH") {
oModel <- try(
mseFH(
formula = dfData$gEst_Direct ~ 1,
vardir = dfData$gMSE_Direct,
method = "REML"
)
)
if (!inherits(oModel, "try-error")) {
agEst <- as.vector(oModel$est$eblup)
agMSE <- oModel$mse
gSCorr <- NA
gAIC <- oModel$est$fit$goodness[2]
gBIC <- oModel$est$fit$goodness[3]
}
}
if (sMethod == "saeSFH") {
oModel <- try(
mseSFH(
formula = dfData$gEst_Direct ~ 1,
vardir = dfData$gMSE_Direct,
proxmat = mgNB,
method = "REML"
)
)
if (!inherits(oModel, "try-error")) {
agEst <- as.vector(oModel$est$eblup)
agMSE <- oModel$mse
gSCorr <- oModel$est$fit$spatialcorr
gAIC <- oModel$est$fit$goodness[2]
gBIC <- oModel$est$fit$goodness[3]
}
}
if (sMethod == "spdepWGT") {
oModel <- try(
errorsarlm(
gEst_Direct ~ 1 ,
data = dfData,
listw = mat2listw(mgNB, style="W"),
weights = 1/gMSE_Direct
)
)
## print(oModel)
if (!inherits(oModel, "try-error")) {
agEst <- as.vector(predict(oModel, pred.type="TS"))
agMSE <- rep(NA, nrow(dfData))
gSCorr <- oModel$lambda
gAIC <- NA
gBIC <- NA
}
}
if (sMethod == "spdepNOWGT") {
oModel <- try(
errorsarlm(
gEst_Direct ~ 1 ,
data = dfData,
listw = mat2listw(mgNB, style="W")
)
)
## print(oModel)
if (!inherits(oModel, "try-error")) {
agEst <- as.vector(predict(oModel, pred.type="TS"))
agMSE <- rep(NA, nrow(dfData))
gSCorr <- oModel$lambda
gAIC <- NA
gBIC <- NA
}
}
## common process ...
if (!inherits(oModel, "try-error")) {
out <- list(
sMethod = sMethod,
bStatus = TRUE,
agEst = agEst,
agMSE = agMSE,
gSCorr = gSCorr,
gAIC = gAIC,
gBIC = gBIC,
gTime = (proc.time() - t)[3]
)
if (bVerbose) cat("[sub_EstimateDataset2.1] Success.", out$gTime, "sec.\n")
} else {
out <- list(
sMethod = sMethod,
bStatus = FALSE,
agEst = rep(NA, nrow(dfData)),
agMSE = rep(NA, nrow(dfData)),
gSCorr = NA,
gAIC = NA,
gBIC = NA,
gTime = (proc.time() - t)[3]
)
if (bVerbose) cat("[sub_EstimateDataset2.1] Fail.", out$gTime, "sec.\n")
}
return(out)
}
d.sub_EstimateDataset2.1 <- function(
dfData, mgNB, asMethod=c("saeFH", "saeSFH", "spdepWGT", "spdepNOWGT"), ...
){
## purpose: estimate spatial FH model for a dataset by multiple methods
## arg:
## dfData: a data set
## mgNB: neighbor matrix
## asMethod: {"saeFH", "saeSFH", "spdepWGT", "spdepNOWGT"}
## ...: pass to sub_EstimateDatasetByMethod2.1() (bVerbose)
## return:
## a data frame
lOut <- lapply(
asMethod,
function(sMethod){
lResult <- d.sub_EstimateDatasetByMethod2.1(dfData, mgNB, sMethod, ...)
out <- data.frame(
nArea = dfData$nArea,
sName = dfData$sName,
sMethod = rep(lResult$sMethod, nrow(dfData)),
bStatus = rep(lResult$bStatus, nrow(dfData)),
gEst = lResult$agEst,
gMSE = lResult$agMSE,
gSCorr = rep(lResult$gSCorr, nrow(dfData)),
gAIC = rep(lResult$gAIC, nrow(dfData)),
gBIC = rep(lResult$gBIC, nrow(dfData)),
gTime = rep(lResult$gTime, nrow(dfData)),
stringsAsFactors = FALSE
)
return(out)
}
)
out <- bind_rows(lOut) %>%
gather(sVar, gValue, c(bStatus, gEst, gMSE, gSCorr, gAIC, gBIC, gTime)) %>%
mutate(sVar = paste0(sVar, "_", sMethod)) %>%
dplyr::select(nArea, sVar, gValue) %>%
spread(sVar, gValue)
out <- full_join(dfData, out, by = "nArea")
return(out)
}
d.sub_EstimateMultiDataset2.1 <- function(dfMultiData, mgNB, ...){
## purpose: estimate FH models for multiple datasets by multiple methods
## arg:
## dfMultiData: multiple data
## mgNB: neighbor matrix
## ...: pass to sub_EstimateDataset2.1() (asMethod, bVerbose)
## return:
## a data frame of multiple results
out <- dfMultiData %>%
group_by(nDataset) %>%
dplyr::do(d.sub_EstimateDataset2.1(., mgNB = mgNB, ...)) %>%
ungroup()
return(out)
}
d.sub_EstimatorAsFactor.1 <- function(asEstimator){
## purpose: convert estimator names into factor
## args:
## asEstimator: a character vector of estimator names
## return:
## a factor vector
factor(
asEstimator,
levels = c(
"True", "Direct", "saeFH", "saeSFH",
"spdepWGT", "spdepNOWGT"
),
labels = c(
"True", "Direct", "sae(w/o sc)", "sae(w/ sc)",
"spdep(wgt)", "spdep(unwgt)"
)
)
}
d.sub_EstimatorPalette.1 <- function(afEstimator){
## purpose: get consisntent color pallette for estimators
## args:
## afEstimator: a factor vector of estimator
## return:
## a vector of color palette
gg_color_hue(6)[sort(unique(as.integer(afEstimator)))]
}
ためしに、上記のデータ生成関数とデータ分析関数を使ってみよう。
地域数を\(10 \times 10 = 100\), 地域効果の分散\(\sigma_u^2\)を50, 空間自己相関パラメータ\(\phi\)を0.7とする。
set.seed(100)
d.mgNeighbor <- d.sub_MakeNeighborMatrix.1(10, 10)
d.dfDemoData <- d.sub_GenerateDataset2.1(mgNB = d.mgNeighbor, gVarU=50, gPhi=0.7)
各地域の標本サイズと、幸福度の標本平均を描いてみよう。
dfPlot <- d.dfDemoData %>%
dplyr::select(nArea, nGridRow, nGridCol, nSize, gEst_Direct)
g <- ggplot(data=dfPlot, aes(x=nGridCol, y=-nGridRow, size = nSize, label=nArea))
g <- g + geom_count(aes(color=gEst_Direct), shape=15)
g <- g + geom_text(size=4)
g <- g + scale_size_area(max_size = 12)
g <- g + scale_x_discrete()
g <- g + scale_y_discrete()
g <- g + scale_color_gradient(low="yellow", high="red")
g <- g + labs(x=NULL, y=NULL, size="SampleSize", color="Direct")
g <- g + theme_bw()
print(g)
正方形は各地域、数字は地域番号, サイズは標本サイズ、色の濃さは標本平均の高さを表す。 心なしか、隣接している地域の色は似ているような気がしませんか?
MorranのIを求めてみると…
moran.test(d.dfDemoData$gEst_Direct, mat2listw(d.mgNeighbor, style="W"))
Moran I test under randomisation
data: d.dfDemoData$gEst_Direct
weights: mat2listw(d.mgNeighbor, style = "W")
Moran I statistic standard deviate = 4.5121, p-value = 3.209e-06
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.323555742 -0.010101010 0.005468176
I=0.32。有意ではあるが、そんなに大きくない。これは標本平均に標本抽出誤差が 載っているからだろうと思い、ためしに真値\(v_i\)について MorranのIを求めてみたら…
moran.test(d.dfDemoData$gTrue, mat2listw(d.mgNeighbor, style="W"))
Moran I test under randomisation
data: d.dfDemoData$gTrue
weights: mat2listw(d.mgNeighbor, style = "W")
Moran I statistic standard deviate = 4.864, p-value = 5.753e-07
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.35029214 -0.01010101 0.00548998
I=0.35, たいしてかわらない。そんなもんなのか…
では、推定してみましょう。 次の6個を比較する。
rstanも試したいところだが、計算時間が長すぎるので断念した。
dfResult <- d.sub_EstimateDataset2.1(d.dfDemoData, d.mgNeighbor, bVerbose=F)
dfPlot <- dfResult %>%
filter(nSize == 10) %>%
dplyr::select(nArea, nSize, gTrue, starts_with("gEst_")) %>%
arrange(-gTrue) %>%
mutate(fArea = factor(nArea, levels=nArea, labels=as.character(nArea))) %>%
rename(gEst_True = gTrue) %>%
gather(sVar, gValue, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
mutate(fEstimator = d.sub_EstimatorAsFactor.1(sEstimator))
g <- ggplot(data=dfPlot, aes(x=fArea, y=gValue, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + labs(
x = "area (of sample size == 10)",
y = expression(paste(theta, " or ", hat(theta))),
color = "Estimator"
)
g <- g + scale_color_manual(values = d.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
上の図は、標本サイズ10の地域のみについて、真値\(\theta_i = v_i\)ならびにその推定値を 示したものである。真値が高い順に並べてある。
spdepによる推定値は、真値よりもずっと全体平均に近い推定値となっている。 はっきりいってぼろ負けである。
saeによる推定値は直接推定(標本平均)に近い。また、空間相関を考慮しようがしまいが、 結果はたいして変わらない。
… もっとも、この結果はたまたまかもしれない。 同じことを何度も繰り返してみないと、なんともいえない。
お待たせしました。繰り返します。
さきほどと同じく、 地域数を\(10 \times 10 = 100\), 地域効果の分散\(\sigma_u^2\)を50, 地域効果の分散\(\sigma_u^2\)を50, 空間自己相関パラメータ\(\phi\)を0.7とする。 500回繰り返す。
# set.seed(100)
# d.dfSimData1 <- d.sub_GenerateMultiDataset2.1(
# nNumDataset = 500,
# mgNB = d.mgNeighbor,
# gVarU = 50,
# gPhi = 0.7
# )
# d.dfSimulation1 <- d.sub_EstimateMultiDataset2.1(
# dfMultiData = d.dfSimData1,
# mgNB = d.mgNeighbor
# )
# save(d.dfSimulation1, file="./d_dfSimulation1.RData")
load("./d_dfSimulation1.RData")
推定値はどのくらいあてになるか。
dfPlot <- d.dfSimulation1 %>%
dplyr::select(nSize, gTrue, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
mutate(
gDiff = gEst - gTrue,
fEstimator = d.sub_EstimatorAsFactor.1(sEstimator)
)
g <- ggplot(data=dfPlot, aes(x = nSize, y=gDiff, color=fEstimator))
g <- g + geom_point(alpha=1/5)
g <- g + labs(
x = "Sample Size",
y = expression(hat(theta)-theta),
color = "Estimator"
)
g <- g + geom_hline(yintercept = 0)
g <- g + scale_x_continuous(breaks = seq(10, 100, by = 20))
g <- g + scale_color_manual(
guide = FALSE,
values = d.sub_EstimatorPalette.1(dfPlot$fEstimator)
)
g <- g + facet_grid(~ fEstimator)
g <- g + theme_bw()
print(g)
縦軸は、推定値と真値との差を表している。
まずみてとれるのは、次の点である: spdep!! おまえはダメな子!!
ま、spdepのerrorsarlm()というのは空間Fay-Herriotモデルとは違うのだろう。 ないし、私の使い方が間違っているのだろう。 ともあれ、以下ではspdepの結果は省略する。
推定値と正解との差を集計すると、
dfPlot <- d.dfSimulation1 %>%
dplyr::select(nSize, gTrue, gEst_Direct, gEst_saeFH, gEst_saeSFH) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(sEstimator, nSize) %>%
summarize(
nFreq = n(),
gRMSE = sqrt(mean((gEst - gTrue)^2))
) %>%
ungroup() %>%
mutate(
fEstimator = d.sub_EstimatorAsFactor.1(sEstimator)
)
g <- ggplot(data=dfPlot, aes(x = nSize, y=gRMSE, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + scale_x_continuous(breaks = seq(10, 100, by = 10))
g <- g + labs(x = "Sample Size", y = "observed RMSE", color="Estimator")
g <- g + scale_color_manual(values = d.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
縦軸は「推定値と正解との差の二乗の平均の平方根」(RMSE)を示す。平たく言っちゃうと、 推定した地域平均は、真の地域平均からだいたい何点くらいずれるか、を示している。
標本サイズが10のときには、だいたい3点くらいずれてしまう。 標本サイズが大きくなると、ずれは小さくなる。
ここで関心があるのは推定量の間の差だ。拡大します。
dfPlot <- d.dfSimulation1 %>%
dplyr::select(nSize, gTrue, gEst_Direct, gEst_saeFH, gEst_saeSFH) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(sEstimator, nSize) %>%
summarize(
gRMSE = sqrt(mean((gEst - gTrue)^2))
) %>%
ungroup() %>%
mutate(
gRatio = gRMSE / sqrt(100/nSize),
fEstimator = d.sub_EstimatorAsFactor.1(sEstimator)
)
g <- ggplot(data=dfPlot, aes(x = nSize, y=gRatio, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + geom_hline(yintercept = 1)
g <- g + scale_x_continuous(breaks = seq(10, 100, by = 10))
g <- g + labs(
x = "Sample Size",
y = "(observed RMSE)/(theoretical RMSE of Direct)",
color = "Estimator"
)
g <- g + scale_color_manual(values = d.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
標本抽出誤差は実は100なので、 標本平均の真のRMSEは\(\sqrt{100/n_i}\)である。 これを基準に、それぞれのパッケージによる推定値の RMSEを比で表している。下にいくほど、推定値と真の値のずれが小さい。
このように、空間相関を考慮することで、推定値と真の値のずれが小さくなる ことがわかる。標本サイズが小さいときに効果が大きい。
MSEについてもみておこう。
dfPlot <- d.dfSimulation1 %>%
dplyr::select(
nSize, gTrue,
gEst_Direct, gEst_saeFH, gEst_saeSFH,
gMSE_Direct, gMSE_saeFH, gMSE_saeSFH
) %>%
gather(sVar, gValue, c(starts_with("gEst_"), starts_with("gMSE_"))) %>%
separate(sVar, c("sVar", "sEstimator")) %>%
spread(sVar, gValue) %>%
mutate(
gTrueSE = (gEst - gTrue)^2,
gDiff = gMSE - gTrueSE
) %>%
group_by(sEstimator, nSize) %>%
summarize(
gMeanDiff = mean(gDiff)
) %>%
ungroup() %>%
mutate(
fEstimator = d.sub_EstimatorAsFactor.1(sEstimator)
)
g <- ggplot(data=dfPlot, aes(x = nSize, y=gMeanDiff, group=fEstimator, color=fEstimator))
g <- g + geom_path(size=1)
g <- g + geom_hline(yintercept = 0)
g <- g + scale_x_continuous(breaks = seq(10, 100, by = 10))
g <- g + labs(
x = "Sample Size",
y = "mean ((Estimated MSE)-(Actual Sq.Err))",
color = "Estimator"
)
g <- g + scale_color_manual(values = d.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
縦軸は、「推定されたMSEから、推定値と正解との差の二乗を引いた値」の平均である。 つまり、正の値はMSEの推定が概して過大であること、 負の値は概して過小であることを表す。
3.3.4で確認したように、 標本サイズが小さい時、BLUP推定のMSEの推定は過小になる。 ここでわかるのは、空間相関を考慮するとさらに過小になるということである。
つまりこういうことだ。 標本サイズが小さいとき、空間相関を正しく考慮することで、真値と推定値との ずれは小さくなる(推定量の真のMSEは下がる)。いっぽう、真値と推定値との ずれに対する推定(推定量のMSEの推定)は、もともと過度に楽観的だったのが、 さらに楽観的になってしまうわけだ。へー!
今度は、空間相関の強さと推定の精度の関係について調べてみたい。
地域効果の分散\(\sigma_u^2\)を50, 空間自己相関パラメータ\(\phi\)を0.7とする。 空間自己相関パラメータ\(\phi\)を0, 0.1, 0.2, …, 0.8と変えながら、 データ生成と分析を500試行繰り返した。 推定方法はsaeパッケージのみ、2種類とする。
# set.seed(123)
# lOut <- lapply(
# seq(from=0, to=0.8, by = 0.1),
# function(gPhi){
# cat("gPhi = ", gPhi, "...")
# out <- d.sub_GenerateMultiDataset2.1(
# nNumDataset = 500,
# mgNB = d.mgNeighbor,
# gVarU = 50,
# gPhi = gPhi
# )
# out$gPhi <- gPhi
# cat("\n")
# return(out)
# }
# )
# d.dfSimData2 <- bind_rows(lOut)
#
# d.dfSimulation2 <- d.dfSimData2 %>%
# group_by(gPhi) %>%
# dplyr::do(
# d.sub_EstimateMultiDataset2.1(
# dfMultiData = .,
# mgNB = d.mgNeighbor,
# asMethod = c("saeFH", "saeSFH")
# )
# ) %>%
# ungroup()
# save(d.dfSimulation2, file="./d_dfSimulation2.RData")
load("./d_dfSimulation2.RData")
まずは推定の精度から。ここでは、 標本サイズ10の地域のみに注目する。各試行について10地域を利用できる。
dfPlot <- d.dfSimulation2 %>%
filter(nSize == 10) %>%
dplyr::select(gPhi, gTrue, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(gPhi, sEstimator) %>%
summarize(
gRMSE = mean((gTrue-gEst)^2)
) %>%
ungroup() %>%
mutate(
fEstimator = d.sub_EstimatorAsFactor.1(sEstimator)
)
## print(dfPlot)
g <- ggplot(data=dfPlot, aes(x=gPhi, y = gRMSE, color=fEstimator, group=fEstimator))
g <- g + geom_path(size=1)
g <- g + labs(x = expression(phi), y = "observed RMSE", color="Estimator")
g <- g + scale_color_manual(values = d.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + theme_bw()
print(g)
縦軸は「推定値と正解との差の二乗の平均の平方根」(RMSE)を示す。平たく言っちゃうと、 推定した地域平均は、真の地域平均からだいたい何点くらいずれるか、を示している。
空間自己相関パラメータ\(\phi\)が0.2以下のときは、 空間相関を考慮しないほうが、推定値と真値のずれがわずかに小さい。 いっぽう0.3以上になると、空間相関を考慮することによって、推定値の精度が高くなる。
では、\(\phi\)の大きさを知ることはできるのか。 こんどは\(\phi\)の推定値について調べてみよう。
dfPlot1 <- d.dfSimulation2 %>%
filter(nArea == 1) %>%
dplyr::select(gPhi, gSCorr_saeSFH)
dfPlot2 <- dfPlot1 %>%
group_by(gPhi) %>%
summarize(gMean = mean(gSCorr_saeSFH)) %>%
ungroup()
g <- ggplot(data=dfPlot2, aes(x = gPhi, y = gMean))
g <- g + geom_abline(intercept = 0, slope = 1)
g <- g + geom_violin(data=dfPlot1, aes(y=gSCorr_saeSFH, group = gPhi))
# g <- g + geom_point(data=dfPlot1, aes(y=gSCorr_saeSFH), size=1, alpha=1/10)
g <- g + geom_point(size=5)
g <- g + labs(x = expression(phi), y = expression(paste("estimated ", phi)))
g <- g + theme_bw()
print(g)
横軸は空間自己相関パラメータの真値\(\phi\)、縦軸はsaeパッケージ(空間Fay-Herriotモデル)に よるその推定値。 分布(バイオリンプロット)と平均(黒丸)を示している。
\(\phi\)の推定は、平均すると正しいが、かなりばらつきが大きい。\(\phi\)が低い場合には、 さらにばらつきが大きくなる。 つまり、標本から\(\phi\)を知るのは難しいわけだ。
具体的な分析の場面で、空間相関を考慮するかどうかを検討する際には、 空間Fay-Herriotモデルを推定して\(\phi\)の推定値を調べるよりも、 空間相関を考慮するモデルと考慮しないモデルのAICを比較することが 多いのではないかと思う。試してみよう。
dfPlot <- d.dfSimulation2 %>%
filter(nArea == 1) %>%
mutate(bWin = if_else(gAIC_saeSFH < gAIC_saeFH, 1, 0)) %>%
group_by(gPhi) %>%
summarize(gWin = mean(bWin)) %>%
ungroup()
g <- ggplot(data=dfPlot, aes(x = gPhi, y = gWin))
g <- g + geom_col()
g <- g + labs(x = expression(phi), y = "Ratio of selecting spatial model")
g <- g + theme_bw()
print(g)
縦軸は、空間相関を考慮しないFay-Herriotモデルと空間Fay-Herriotモデルを比較したとき、 後者のほうがAICが小さくなった試行の割合である。
この図からわかるように、真の空間相関がゼロでも、 空間Fay-HerriotモデルのAICのほうが低くなってしまうことも ある。また、真の空間相関が小さいとき、空間Fay-HerriotモデルのAICのほうが 必ず低くなるとも限らない。 要するに、空間相関がない、ないし小さいとき、そのことを正しく知るのは 結構難しい。
ついでに、空間相関の有無についての判断を誤り、うっかり変なモデルを採用しちゃったら なにが起きるか試してみよう。
dfPlot <- d.dfSimulation2 %>%
filter(nSize == 10, gPhi < 0.4 ) %>%
mutate(bWin = if_else(gAIC_saeSFH < gAIC_saeFH, 1, 0)) %>%
dplyr::select(gPhi, gTrue, bWin, starts_with("gEst_")) %>%
gather(sVar, gEst, starts_with("gEst_")) %>%
separate(sVar, c("sVar1", "sEstimator")) %>%
group_by(gPhi, sEstimator, bWin) %>%
summarize(
gRMSE = mean((gTrue-gEst)^2)
) %>%
ungroup() %>%
mutate(
fWin = factor(bWin, levels=c(0,1), labels=c("w/o sc", "w/ sc")),
fEstimator = d.sub_EstimatorAsFactor.1(sEstimator)
)
# print(dfPlot)
g <- ggplot(data=dfPlot, aes(x=fWin, y = gRMSE, color=fEstimator, group=fEstimator))
g <- g + geom_path(size=1)
g <- g + labs(x = "selected model", y = "observed RMSE", color="Estimator")
g <- g + scale_color_manual(values = d.sub_EstimatorPalette.1(dfPlot$fEstimator))
g <- g + facet_grid(~gPhi)
g <- g + theme_bw()
print(g)
パネルは真の空間自己相関パラメータ。パネル内の左側は、 AICによって空間相関を考慮しないモデルが支持された試行、 右側は、AICによって空間相関を考慮するモデルが支持された試行。縦軸は 「推定値と正解との差の二乗の平均の平方根」(RMSE)で、小さいほうが 推定値の精度が高い。
さきほどの分析では、 真の空間自己相関パラメータが0.2までであれば空間相関を考慮しないほうが いいし、それより大きければ空間相関を考慮したほうがいい、ということが わかった。この図からは、手元のデータから得られたAICがどうであれ (つまり、パネル内の左側の状況であれ右側の状況であれ)、 同じことがいえるのだ、と確認できる。
このシミュレーションの状況に関する限り、 空間相関が弱い(\(\phi\)が0.2以下だ)という確信があるならば、 たとえAICによって空間Fay-Herriotモデルが支持されようが、 断固として、ただのFay-Herriotモデルを使うべきなのである。
このシミュレーションでわかったことをまとめておく。
rm(list=grep("^d\\.", ls(), value=T))
いったんこれで終了とします。続きはまた後日…