library(diagram)
という pacakge があると教えていただいた
(RGM による example 図一覧)．
なかなか便利そう．
ここではある GLSM (generalized linear spatial model) の推定のために library(geoRglm) の glsm.mcmc() 関数 と likfit.glsm() 関数が一度ずつ使われていますが，そも そもこのような使いかたで良いのかどうか疑問です． というのも，これらは ・glsm.mcmc(): あるパラメーターを与えられた GLSM の ものとでのガウス確率場の MCMC によるシミュレイション ・likfit.glsm(): あるガウス確率場の実現値のセットを与えら れたときの GLSM のパラメーターの最尤推定 という関数なので，それぞれ一度ずつ実行すればよいと いったものではないように思います．うーむ， 今まで見せていただいた例 (他の人たち) も一度ずつしか使ってなかったような ……
おそらく以下のような繰り返し計算 (これは一種の EM アル
ゴリズムだと思うのですが) をする必要があるでしょう．
(1) 係数・分散構造の vector である beta, cov.pars にてきとうな
値をいれる
(2) glsm.mcmc() によるガウス確率場のサンプリング
(3) likfit.glsm() による GLSM のパラメーター推定
(4) この最尤推定値で beta と cov.pars をおきかえて (2) にもどる
といったくりかえしです．
実際に，下記のような R コードを準備して
library(geoRglm)
D <- read.csv("GeoP.csv")
geoD <- as.geodata(
D,
coords.col = c(2, 3),
data.col = 4,
covar.col = c(5, 6, 7)
)
beta <- c(-15, 0.03, 0.02, 0.01)
cov.pars <- c(0.6, 0.1) # c(sigma^2, phi)
file <- "kubo.results.txt"
sink(file); sink() # to clear file
for (i in 1:5) {
mcmc1 <- mcmc.control(
S.scale = 0.6,
thin = 20,
n.iter = 10000,
burn.in = 1000
)
model1 <- list(
cov.pars = cov.pars,
beta = beta,
trend = ~ Covar1 + Covar2 + Covar3,
family = "poisson"
)
outmcmc1 <- glsm.mcmc(geoD, model = model1, mcmc.input = mcmc1)
mcmcobj1 <- prepare.likfit.glsm(outmcmc1)
lik1 <- likfit.glsm(mcmcobj1, ini.phi = cov.pars[2])
beta <- lik1$beta
cov.pars <- lik1$cov.pars
sink(file, append = TRUE)
cat("\n# Iteration", i, "\n")
print(lik1)
sink()
}
5 回ほど (2)-(4) の過程を繰りかえしてみると，このような
結果がえられました:
# Iteration 1
likfit.glsm: estimated model parameters:
beta0 beta1 beta2 beta3 sigmasq phi tausq.rel
"-18.2325" " 0.0871" " 0.0512" " 0.0308" " 7.8854" " 18.6653" " 0.6107"
likfit.glsm : maximised log-likelihood = 4010
# Iteration 2
likfit.glsm: estimated model parameters:
beta0 beta1 beta2 beta3 sigmasq phi tausq.rel
"-7.8130" " 0.0735" " 0.0459" " 0.0239" " 3.1119" " 6.0664" " 0.0488"
likfit.glsm : maximised log-likelihood = 93.6
# Iteration 3
likfit.glsm: estimated model parameters:
beta0 beta1 beta2 beta3 sigmasq phi tausq.rel
"-6.2595" " 0.0598" " 0.0375" " 0.0185" " 2.9844" "10.6527" " 0.0648"
likfit.glsm : maximised log-likelihood = 52.2
# Iteration 4
likfit.glsm: estimated model parameters:
beta0 beta1 beta2 beta3 sigmasq phi tausq.rel
"-5.3985" " 0.0530" " 0.0333" " 0.0151" " 3.0246" "13.2617" " 0.0172"
likfit.glsm : maximised log-likelihood = 28.3
# Iteration 5
likfit.glsm: estimated model parameters:
beta0 beta1 beta2 beta3 sigmasq phi tausq.rel
"-4.8941" " 0.0472" " 0.0305" " 0.0133" " 2.6217" "13.2617" " 0.0096"
likfit.glsm : maximised log-likelihood = 18.0
試行回数とともにパラメーターが変わっていきますが，
次第にあるところに落ちつきつつあるようです．
なかなか時間がかかるので (とくに likfit.glsm() に
よる最尤推定) 5 回のくりかえしで終了しましたが，
10-20 回ぐらいくりかえせば最尤推定値らしきセット
が得られるのではないでしょうか．
いぜんとして最大化対数尤度はゼロより大ですが，
このあたりは尤度の計算方法が glm() のそれとは異なる
ということで単純に比較できないかもしれません．geoR
系の文献やソースコードを読んだりするのはたいへんそう
なので，outmcmc1 や lik1 に格納されている推定値を使っ
て，自分で対数尤度を計算したほうが良いかもしれません．
不完全ながら，とりあえずのご報告です．
どんなもんでしょうね?
library(geoRglm)
の理解が少しだけススんだ
……
こういう進歩があったりするので，
ついついデータ解析こんさるに時間を費やしてしまう．
No new mail!
の文字が
(めんどうあれこれを Starred におしこんだだけ，
とゆー説もあるけど)
……
やれやれ．