ネットワーク分析に明るい方であれば、2019年に書かれたこの野心的な記事を覚えている人も多いのではないだろうか。
statnet チームが提供している R のパッケージに頼らずに、Exponential Random Graph Model (ERGM) をできるだけスクラッチで実装しようという壮大な試みである。
実装した関数による推定値は、残念ながら {ergm} で推定した値と乖離してしまったという結びになっている。
その原因解明は多くの読者にとって待ち望んでいたことに違いない。
そこで時間を持て余した筆者が、ゴールデンウィークを使って ERGM を (MPLE を除いて) 一から、かつ {Rcpp} を用いて高速に推定できるよう実装した。*1
以下の repository にコードをパッケージとして公開しているので、興味のある方は遊んでみていただけると幸いである。
flomariagge ~ edges + triangle という単純なモデルを {myergm}、{ergm} で 100 回ずつ推定し、その推定値の分布をプロットしたものが以下の図である。
両者は概ね一致しているように思われる。
先の記事ではなぜ推定がうまくいかなかったのか？
気になるのが、一番最初にあげた記事でなぜ ERGM の推定値が {ergm} のそれと比較して大きく異なってしまったのかである。
その原因が、以下の関数にあることを特定した。
#対数尤度比を計算する関数 loglikelihood_ratio <- function(theta,fit,n){ sample_length <- length(fit$edges) sample_index <- (sample_length-n):sample_length sampled_edge_count <- fit$edges[sample_index] sampled_triangle_count <- fit$triangle[sample_index] theta0 <- c(fit$theta0_edges,fit$theta0_triangle) g_obs <- c(count_edges(fit$obs_graph),count_triangle(fit$obs_graph)) first = (theta - theta0) %*% g_obs exp_vec <- c() for(i in 1:n){ sam <- exp((theta - theta0) %*% c(sampled_edge_count[i],sampled_triangle_count[i])) exp_vec <- c(exp_vec,sam) } second <- log(sum(exp_vec)/length(exp_vec)) loglikelihood_ratio <- first - second return(loglikelihood_ratio) }
sample_index <- (sample_length-n):sample_length とあるので、log likelihood ratio を計算するときに、MCMC で生成されたネットワークの数を としたとき、
番目から
番目までの連続した全てのサンプルを使っていることになる。
一方、{ergm} には control.ergm() という関数があり、MCMC をチューニングできるが、そのデフォルトは以下のようになっている。
MCMC.interval = 1024, MCMC.burnin = MCMC.interval * 16, MCMC.samplesize = 1024,
まず、1024 * 16 を burnin とする。
その後、MCMC で生成されたグラフ同士の間隔が 1024 となるように (MCMC.interval = 1024)、1024 個のグラフをサンプリングする (MCMC.samplesize = 1024)。
そのサンプリングされた 1024 個のグラフを推定値を得るために用いるのである。
これを踏まえると、最初の記事で推定がうまくいかなかった原因は、MCMC でサンプリングされたグラフ同士の間隔が 0 であり、グラフ同士の serial correlation が強いことに原因がありそうである。
そこで、{myergm} でグラフ同士の相関が強く残るようにサンプリングして推定した場合の推定値の分布を見てみよう。
推定がうまくいかなかったケースを、仮に arabou_ergm と呼ぶことにすると、それは以下のように実装できる。
# Number of simulations R <- 100 df <- foreach(i = 1:R, .combine = "rbind") %do% { # Estimate the parameters. theta_arabou_ergm <- myergm::myergm_MCMLE(model = flomarriage ~ edges + triangle, seed = i, p_one_node_swap = 0, p_large_step = 0, p_invert = 0, MCMC_samplesize = 1024, MCMC_interval = 1, second_round = FALSE) # Return the output. return(data.frame("iter" = i,"param" = c("edges", "triangle"), "arabou_ergm" = theta_arabou_ergm)) }
{ergm} で推定したケースと比較したとき、それぞれの分布は以下の通りになる。
arabou_ergm の方が推定値にかなりばらつきが出ていることがわかる。
私が実装した {myergm} は、厳密には {ergm} の実装とは異なる点がある。
例えば、{ergm} の MCMC のデフォルトの sampler は "tie-no-tie" といって、最初にリンクが存在する dyad か、リンクが存在しない dyad どちらの状態を変化させるかをそれぞれ確率 1/2 で決定し、リンクを加える/消す。
一方 {myergm} では、tie-no-tie sampler は実装しておらず、代わりに Mele (2017) が提案したような large step を許容した sampler、つまりある確率で 1 回のステップで複数の dyad の状態を変化させるような実装になっている。
ただし、これらの MCMC sampling の手法は推定値に大きく影響を与えるものではなく、大事なのは推定に用いるグラフの sampling の方であったようだ。 他にも ERGM には複雑怪奇なことが多いので*2、それらは今後の検討課題としたい。
References
- Hummel, R. M., Hunter, D. R., & Handcock, M. S. (2012). Improving simulation-based algorithms for fitting ERGMs. Journal of Computational and Graphical Statistics, 21(4), 920-939.
- Hunter, D. R., Goodreau, S. M., & Handcock, M. S. (2013). ergm. userterms: A Template Package for Extending statnet. Journal of statistical software, 52(2), i02.
- Mele, A. (2017). A structural model of dense network formation. Econometrica, 85(3), 825-850.
