knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
The ergmito
's workhorse are two other functions: (1) ergm
's ergm.allstats
which is used to compute the support of model's sufficient statistics, and (2)
ergmito_formulae
which is a wrapper of that same function, and that returns a
list including the following two functions: loglike
and grad
, the functions
to calculate the joint log-likelihood of the model and its gradient.
library(ergmito) data(fivenets) model_object <- ergmito_formulae(fivenets ~ edges + ttriad) # Printing the model object model_object # Printing the log-likelihood function model_object$loglik
Besides of the log-likelihood function and the gradient function,
ergmito_formulae
also returns We can take a look at each component from our previous object:
# The vectors of weights str(model_object$stats_weights) # The matrices of the sufficient statistics str(model_object$stats_statmat) # The target statistic model_object$target_stats
All this is closely related to the output object from the function ergm.allstats
.
The next section shows how all this works together to estimate the model parameters
using Metropolis-Hastings MCMC.
Suppose that we have a prior regarding the distribution of the fivenets
model:
The edges
parameter is normally distributed with mean -1 and variance 2, while
the nodematch("female")
term has the same distribution but with mean 1. We can
implement this model using a Metropolis-Hastings ratio. First, we need to define
the log of the posterior distribution:
# Analyzing the model model_object <- ergmito_formulae(fivenets ~ edges + nodematch("female")) # Defining the logposterior logposterior <- function(p) { model_object$loglik(params = p) + sum(dnorm(p, mean = c(-1,1), sd = sqrt(2), log = TRUE)) }
For this example, we are using the fmcmc R package. Here is how we put everything together:
# Loading the required R packages library(fmcmc) library(coda) # Is it working? logposterior(c(-1, 1)) # Now, calling the MCMC function from the fmcmc R package ans <- MCMC( fun = logposterior, initial = c(0, 0), # 5,000 steps sampling one of every ten iterations nsteps = 5000, thin = 10, # We are using a normal transition kernel with .5 scale and updates are done # one variable at a time in a random order kernel = kernel_normal(scale = .5, scheme = "random") )
We can now visualize our results. Since the resulting object is of class mcmc.list
,
which is implemented in the coda
R package for MCMC diagnostics, we can use
all the methods included in the package:
plot(ans) summary(ans)
Finally, we can compare this result with what we obtain from the ergmito
function
summary(ergmito(fivenets ~ edges + nodematch("female")))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.