Extending ergmito

In ergmito: Exponential Random Graph Models for Small Networks

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.

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
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")))

ergmito documentation built on July 9, 2023, 7:09 p.m.