Extending ergmito

  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.

model_object <- ergmito_formulae(fivenets ~ edges + ttriad)

# Printing the model object

# Printing the log-likelihood function

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

# The matrices of the sufficient statistics

# The target statistic

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.

Example: Bayesian inference with fivenets

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

# 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:


Finally, we can compare this result with what we obtain from the ergmito function

summary(ergmito(fivenets ~ edges + nodematch("female")))

Try the ergmito package in your browser

Any scripts or data that you put into this service are public.

ergmito documentation built on Aug. 11, 2020, 1:08 a.m.