Package {Bergm}

Bayesian exponential random graph models

Description

Bergm provides a range of tools to analyse Bayesian exponential random graph models using advanced computational methods: parameter estimation via the approximate exchange algorithm (bergm), calibrated pseudo-posteriors (bergmC), adjusted pseudolikelihoods (ergmAPL), model evidence estimation (evidence), missing data imputation (bergmM) and Bayesian goodness-of-fit (bgof).

Author(s)

Maintainer: Alberto Caimo alberto.caimo1@ucd.ie

Authors:

• Lampros Bouranis

• Robert Krause

Other contributors:

• Nial Friel [contributor]

See Also

Useful links:

• https://acaimo.github.io/Bergm/

Parameter estimation for Bayesian ERGMs

Description

Fits Bayesian exponential random graph models with the approximate exchange algorithm and a population of interacting chains updated by adaptive direction sampling (ADS) (Caimo and Friel, 2011).

Usage

bergm(
  formula,
  prior.mean = NULL,
  prior.sigma = NULL,
  burn.in = 100,
  main.iters = 1000,
  aux.iters = 1000,
  nchains = NULL,
  gamma = 0.5,
  V.proposal = 0.0025,
  startVals = NULL,
  offset.coef = NULL,
  constraints = NULL,
  thin = 1,
  cut.reject = FALSE,
  saveEveryX = NULL,
  saveEveryXName = "partialBergmEstimate.rda",
  seed = NULL,
  adapt = TRUE,
  target.accept = 0.25,
  delayed.acceptance = TRUE,
  da.inflate = 1.5,
  cores = 1,
  log.prior = NULL,
  prior.decay = c(meanlog = log(0.5), sdlog = 0.75),
  verbose = TRUE,
  ...
)

Arguments

Details

Compared with Bergm 5, the sampler draws every auxiliary network through a single call to ergm's compiled sampler on a pre-built network state, supports curved ERGM terms, models with a single parameter, custom priors, burn-in adaptation of the ADS move factor, delayed acceptance and parallel chains. The posterior sample is returned both pooled (Theta, as in previous versions) and per chain (chains), which enables standard convergence diagnostics such as the potential scale reduction factor.

Value

An object of class bergm with elements

plus the model specification and run settings.

References

Caimo, A. and Friel, N. (2011), "Bayesian Inference for Exponential Random Graph Models," Social Networks, 33(1), 41-55. https://arxiv.org/abs/1007.5192

Caimo, A. and Friel, N. (2014), "Bergm: Bayesian Exponential Random Graphs in R," Journal of Statistical Software, 61(2), 1-25. doi:10.18637/jss.v061.i02

Christen, J. A. and Fox, C. (2005), "Markov chain Monte Carlo using an approximation," Journal of Computational and Graphical Statistics, 14(4), 795-810.

Examples

## Not run: 
# Load the florentine marriage network
data(florentine)

# Posterior parameter estimation:
p.flo <- bergm(flomarriage ~ edges + kstar(2),
               burn.in    = 50,
               aux.iters  = 500,
               main.iters = 3000,
               gamma      = 1.2,
               seed       = 1)

# Posterior summaries:
summary(p.flo)

# Curved model, delayed acceptance and 4 parallel workers:
p.laz <- bergm(lazega ~ edges + nodematch("Office") + gwesp(0.5, fixed = FALSE),
               burn.in = 200, main.iters = 2000, aux.iters = 3000,
               delayed.acceptance = TRUE, cores = 4, seed = 1)

## End(Not run)

Posterior point estimates and covariance for bergm objects

Description

Posterior point estimates and covariance for bergm objects

Usage

## S3 method for class 'bergm'
coef(object, ...)

## S3 method for class 'bergm'
vcov(object, ...)

## S3 method for class 'bergm'
as.mcmc(x, ...)

## S3 method for class 'bergm'
as.mcmc.list(x, ...)

Arguments

Value

coef: posterior means; vcov: posterior covariance matrix.

Calibrating misspecified Bayesian ERGMs

Description

Function to transform a sample from the pseudo-posterior to one that is approximately sampled from the intractable posterior distribution (Bouranis, Friel and Maire, 2017).

Usage

bergmC(
  formula,
  prior.mean = NULL,
  prior.sigma = NULL,
  burn.in = 10000,
  main.iters = 40000,
  aux.iters = 3000,
  V.proposal = 1.5,
  thin = 1,
  rm.iters = 500,
  rm.a = 0.001,
  rm.alpha = 0,
  n.aux.draws = 400,
  aux.thin = 50,
  estimate = c("MLE", "CD"),
  seed = 1,
  map.method = c("newton", "rm"),
  map.iters = 20,
  init = NULL,
  verbose = TRUE,
  ...
)

Arguments

Details

The pseudo-posterior is sampled by a compiled random-walk Metropolis kernel. The mode of the true posterior, which anchors the calibration, is found by default with a Newton-type stochastic approximation that uses the simulated covariance of the sufficient statistics and a persistent auxiliary network chain (map.method = "newton"); this typically converges in about ten iterations, whereas the Robbins-Monro recursion of previous versions (map.method = "rm") needs hundreds. The pseudo-posterior mode is found exactly by Newton-Raphson.

Value

An object of class bergm, with additional elements Theta_star (estimated posterior mode), Theta_PL (pseudo-posterior mode), W (curvature adjustment matrix) and map.trajectory.

References

Bouranis, L., Friel, N., & Maire, F. (2017). Efficient Bayesian inference for exponential random graph models by correcting the pseudo-posterior distribution. Social Networks, 50, 98-108. https://arxiv.org/abs/1510.00934

Examples

## Not run: 
# Load the florentine marriage network
data(florentine)

# Calibrated pseudo-posterior:
cpp.flo <- bergmC(flomarriage ~ edges + kstar(2),
                  aux.iters  = 500,
                  burn.in    = 500,
                  main.iters = 10000,
                  V.proposal = 2.5)

# Posterior summaries:
summary(cpp.flo)

## End(Not run)

Parameter estimation for Bayesian ERGMs under missing data

Description

Function to fit Bayesian exponential random graphs models under missing data using the approximate exchange algorithm with Bayesian data augmentation (Koskinen, Robins and Pattison, 2010).

Usage

bergmM(
  formula,
  burn.in = 100,
  main.iters = 1000,
  aux.iters = 1000,
  prior.mean = NULL,
  prior.sigma = NULL,
  nchains = NULL,
  gamma = 0.5,
  V.proposal = 0.0025,
  seed = NULL,
  startVals = NULL,
  offset.coef = NULL,
  constraints = NULL,
  thin = 1,
  saveEveryX = NULL,
  cut.reject = FALSE,
  saveEveryXName = "partialBergmEstimate.rda",
  imputeAllItr = TRUE,
  imputeLast = TRUE,
  nImp = NULL,
  missingUpdate = NULL,
  imputeData = NULL,
  attributeNames = NULL,
  miceIt = 5,
  onlyKeepImputation = FALSE,
  adapt = TRUE,
  target.accept = 0.25,
  log.prior = NULL,
  prior.decay = c(meanlog = log(0.5), sdlog = 0.75),
  verbose = TRUE,
  ...
)

Arguments

Details

Every chain of the population carries its own imputed network, and missing tie variables are updated by a Gibbs-type step from the ERGM restricted to the missing dyads given the chain's current parameter. This is a valid data-augmentation sampler for each chain (in Bergm 5 a single imputed network was shared by all chains). Imputation steps act directly on ergm's compiled network states, which makes imputing at every iteration cheap; network objects are only built for the returned imputations.

References

Caimo, A. and Friel, N. (2011), "Bayesian Inference for Exponential Random Graph Models," Social Networks, 33(1), 41-55. https://arxiv.org/abs/1007.5192

Koskinen, J.H., Robins, G.L., and Pattison, P.E. (2010), "Analysing exponential random graph (p-star) models with missing data using Bayesian data augmentation," Statistical Methodology 7(3), 366-384.

Krause, R.W., Huisman, M., Steglich, C., and Snijders, T.A. (2020), "Missing data in cross-sectional networks-An extensive comparison of missing data treatment methods", Social Networks 62: 99-112.

Examples

## Not run: 
# Load the florentine marriage network
data(florentine)

# Create missing data
set.seed(14021994)
n <- dim(flomarriage[, ])[1]
missNode <- sample(1:n, 1)
flomarriage[missNode, ] <- NA
flomarriage[, missNode] <- NA

# Posterior parameter estimation:
m.flo <- bergmM(flomarriage ~ edges + kstar(2),
                burn.in    = 50,
                aux.iters  = 500,
                main.iters = 1000,
                gamma      = 1.2,
                nImp       = 5)

# Posterior summaries:
summary(m.flo)

## End(Not run)

Bayesian goodness-of-fit diagnostics for ERGMs

Description

Function to calculate summaries for degree, minimum geodesic distances, and edge-wise shared partner distributions to diagnose the Bayesian goodness-of-fit of exponential random graph models.

Usage

bgof(
  x,
  sample.size = 100,
  aux.iters = 10000,
  n.deg = NULL,
  n.dist = NULL,
  n.esp = NULL,
  n.ideg = NULL,
  n.odeg = NULL,
  plot = TRUE,
  cores = 1,
  seed = NULL,
  ...
)

Arguments

Details

For every sampled posterior parameter one network is simulated from the model, starting from the observed network. The goodness-of-fit statistics are monitored inside ergm's compiled sampler (as additional zero-coefficient model terms), so no network objects are created and the posterior predictive distribution of the model statistics is obtained at no extra cost. Simulations can run in parallel with cores > 1.

Value

An object of class bgof (returned invisibly when plotted): simulated (sim.*) and observed (obs.*) proportions for each goodness-of-fit statistic, plus sim.model/obs.model and posterior predictive p-values (pval.model) for the model statistics.

References

Caimo, A. and Friel, N. (2011), "Bayesian Inference for Exponential Random Graph Models," Social Networks, 33(1), 41-55. https://arxiv.org/abs/1007.5192

Caimo, A. and Friel, N. (2014), "Bergm: Bayesian Exponential Random Graphs in R," Journal of Statistical Software, 61(2), 1-25. doi:10.18637/jss.v061.i02

Examples

## Not run: 
# Load the florentine marriage network
data(florentine)

# Posterior parameter estimation:
p.flo <- bergm(flomarriage ~ edges + kstar(2),
               burn.in    = 50,
               aux.iters  = 500,
               main.iters = 1000,
               gamma      = 1.2)

# Bayesian goodness-of-fit test:
bgof(p.flo,
     aux.iters   = 500,
     sample.size = 30,
     n.deg       = 10,
     n.dist      = 9,
     n.esp       = 6)

## End(Not run)

Adjustment of ERGM pseudolikelihood

Description

Function to estimate the transformation parameters for adjusting the pseudolikelihood function (Bouranis, Friel and Maire, 2018).

Usage

ergmAPL(
  formula,
  aux.iters = NULL,
  n.aux.draws = NULL,
  aux.thin = NULL,
  ladder = NULL,
  estimate = c("MLE", "CD"),
  seed = 1,
  mle = NULL,
  cores = 1,
  llk.method = c("bridge", "ladder"),
  bridge.control = NULL,
  ...
)

Arguments

Details

The log-likelihood at the MLE, and hence the adjustment constant, is estimated by default (llk.method = "bridge") by bridge sampling from the dyad-independent submodel of the formula, using ergm.bridge.dindstart.llk. The reference model has an exact likelihood and is close to the target, which keeps the estimate accurate on large networks. llk.method = "ladder" uses the stepping-stone ladder \theta_t = t\,\hat\theta_{MLE}, t \in [0,1], of previous versions, traversed in both directions with warm-started chains and averaging forward and reverse estimates (which removes most of the bias of a single sweep); it needs many rungs when the network is large, because the path starts at the uniform random graph. With cores > 1 the two sweeps run in parallel. The returned object can be passed to evidenceCJ and evidencePP through their apl argument to avoid recomputing it.

Value

An object of class ergmAPL: a list with formula, Theta_MLE, Theta_PL, W, logC, ll_true, logztheta and ll_adjpseudo.

References

Bouranis, L., Friel, N., & Maire, F. (2018). Bayesian model selection for exponential random graph models via adjusted pseudolikelihoods. Journal of Computational and Graphical Statistics, 27(3), 516-528. https://arxiv.org/abs/1706.06344

Examples

## Not run: 
data(florentine)
apl <- ergmAPL(flomarriage ~ edges + kstar(2))
apl

## End(Not run)

Wrapper function for evidence estimation

Description

Function to estimate the evidence (marginal likelihood) with Chib and Jeliazkov's method or Power posteriors, based on the adjusted pseudolikelihood function.

Usage

evidence(evidence.method = c("CJ", "PP"), ...)

Arguments

References

Bouranis, L., Friel, N., & Maire, F. (2018). Bayesian model selection for exponential random graph models via adjusted pseudolikelihoods. Journal of Computational and Graphical Statistics, 27(3), 516-528. https://arxiv.org/abs/1706.06344

Examples

## Not run: 
# Load the florentine marriage network:
data(florentine)

# MCMC sampling and evidence estimation:
CJE <- evidence(evidence.method = "CJ",
                formula     = flomarriage ~ edges + kstar(2),
                main.iters  = 30000,
                burn.in     = 2000,
                aux.iters   = 1000,
                num.samples = 25000,
                V.proposal  = 2.5,
                ladder      = 100,
                seed        = 1)

# Posterior summaries:
summary(CJE)

# MCMC diagnostics plots:
plot(CJE)

# Log-evidence (marginal likelihood) estimate:
CJE$log.evidence

# Compare two models, reusing one pseudolikelihood adjustment per model:
a1 <- ergmAPL(flomarriage ~ edges + kstar(2))
a2 <- ergmAPL(flomarriage ~ edges + triangle)
e1 <- evidence("PP", formula = flomarriage ~ edges + kstar(2), apl = a1)
e2 <- evidence("PP", formula = flomarriage ~ edges + triangle, apl = a2)
e1$log.evidence - e2$log.evidence  # log Bayes factor

## End(Not run)

Evidence estimation via Chib and Jeliazkov's method

Description

Function to estimate the evidence (marginal likelihood) with Chib and Jeliazkov's method, based on the adjusted pseudolikelihood function.

Usage

evidenceCJ(
  formula,
  prior.mean = NULL,
  prior.sigma = NULL,
  aux.iters = 1000,
  n.aux.draws = 50,
  aux.thin = 50,
  ladder = 30,
  main.iters = 30000,
  burn.in = 5000,
  thin = 1,
  V.proposal = 1.5,
  num.samples = 25000,
  seed = 1,
  estimate = c("MLE", "CD"),
  apl = NULL,
  cores = 1,
  verbose = TRUE,
  ...
)

Arguments

Details

The adjusted pseudo-posterior is sampled with a compiled random-walk Metropolis kernel and all pseudolikelihood evaluations needed by the estimator are vectorised in compiled code. The Metropolis-Hastings acceptance probabilities in the Chib-Jeliazkov identity are computed with the full (pseudo-)posterior, i.e. including the prior.

References

Caimo, A., & Friel, N. (2013). Bayesian model selection for exponential random graph models. Social Networks, 35(1), 11-24. https://arxiv.org/abs/1201.2337

Bouranis, L., Friel, N., & Maire, F. (2018). Bayesian model selection for exponential random graph models via adjusted pseudolikelihoods. Journal of Computational and Graphical Statistics, 27(3), 516-528. https://arxiv.org/abs/1706.06344

Examples

## Not run: 
# Load the florentine marriage network:
data(florentine)

# MCMC sampling and evidence estimation:
CJE <- evidenceCJ(flomarriage ~ edges + kstar(2),
                  main.iters  = 2000,
                  burn.in     = 200,
                  aux.iters   = 500,
                  num.samples = 25000,
                  V.proposal  = 2.5)

# Posterior summaries:
summary(CJE)

# MCMC diagnostics plots:
plot(CJE)

# Log-evidence (marginal likelihood) estimate:
CJE$log.evidence

## End(Not run)

Evidence estimation via power posteriors

Description

Function to estimate the evidence (marginal likelihood) with Power posteriors, based on the adjusted pseudolikelihood function.

Usage

evidencePP(
  formula,
  prior.mean = NULL,
  prior.sigma = NULL,
  aux.iters = 1000,
  n.aux.draws = 50,
  aux.thin = 50,
  ladder = 30,
  main.iters = 20000,
  burn.in = 5000,
  thin = 1,
  V.proposal = 1.5,
  seed = 1,
  temps = NULL,
  estimate = c("MLE", "CD"),
  apl = NULL,
  cores = 1,
  verbose = TRUE,
  ...
)

Arguments

Details

Each tempered adjusted pseudo-posterior is sampled with a compiled random-walk Metropolis kernel, and the zero-variance control variates and the corrected trapezoidal rule (Friel, Hurn and Wyse, 2014) are computed with vectorised compiled code.

References

Bouranis, L., Friel, N., & Maire, F. (2018). Bayesian model selection for exponential random graph models via adjusted pseudolikelihoods. Journal of Computational and Graphical Statistics, 27(3), 516-528. https://arxiv.org/abs/1706.06344

Friel, N., Hurn, M., & Wyse, J. (2014). Improving power posterior estimation of statistical evidence. Statistics and Computing, 24(5), 709-723.

Examples

## Not run: 
# Load the florentine marriage network:
data(florentine)

PPE <- evidencePP(flomarriage ~ edges + kstar(2),
                  aux.iters   = 500,
                  aux.thin    = 50,
                  main.iters  = 2000,
                  burn.in     = 100,
                  V.proposal  = 2.5)

# Posterior summaries:
summary(PPE)

# MCMC diagnostics plots:
plot(PPE)

# Log-evidence (marginal likelihood) estimate:
PPE$log.evidence

## End(Not run)

Lazega lawyers network data

Description

Lazega lawyers network data

Usage

lazega

Format

An object of class network.

Source

This network dataset comes from a network study of corporate law partnership that was carried out in a Northeastern US corporate law firm in New England from 1988 to 1991. It represents collaborative relations among the 36 attorneys (partners and associates) of this firm. Nodal attributes include: Age, Gender, Office, Practice, School, and Years.

References

Lazega, E. (2001), "The Collegial Phenomenon: The Social Mechanisms of Cooperation Among Peers in a Corporate Law Partnership," Oxford University Press.

Examples

## Not run: 
 par(mfrow = c(1, 2), oma = rep(0, 4))
 CC <- hcl.colors(3, "Teal")
 set.seed(22)
 plot(lazega,
      vertex.col = CC[lazega %v% "Office"], 
     vertex.cex = 2)
 legend("topright",
        pch    = 21,
        pt.bg  = CC,
        legend = c("Boston", "Hartford", "Providence"),
        title  = "OFFICE")

## End(Not run)

Plot BERGM posterior output

Description

This function creates MCMC diagnostic plots for bergm objects: posterior density, trace plot (one line per chain when available) and autocorrelation for every free parameter.

Usage

## S3 method for class 'bergm'
plot(x, lag.max = 50, ...)

Arguments

Examples

## Not run: 
data(florentine)
p.flo <- bergm(flomarriage ~ edges + kstar(2), main.iters = 1000, seed = 1)
plot(p.flo)

## End(Not run)

Plot Bayesian goodness-of-fit diagnostics

Description

Plot Bayesian goodness-of-fit diagnostics

Usage

## S3 method for class 'bgof'
plot(
  x,
  n.deg = NULL,
  n.dist = NULL,
  n.esp = NULL,
  n.ideg = NULL,
  n.odeg = NULL,
  ...
)

Arguments

Posterior predictive simulation from a bergm fit

Description

Draws parameters from the posterior sample and simulates one network for each, starting from the observed network.

Usage

## S3 method for class 'bergm'
simulate(
  object,
  nsim = 1,
  seed = NULL,
  aux.iters = 10000,
  output = c("network", "stats"),
  cores = 1,
  ...
)

Arguments

Summary of BERGM posterior output

Description

This function summarises MCMC output for bergm objects.

Usage

## S3 method for class 'bergm'
summary(object, returnTable = FALSE, ...)

Arguments

Value

Invisibly (or, with returnTable = TRUE, visibly) the table of posterior summaries. When per-chain samples are available, the table also reports the potential scale reduction factor (Rhat).