The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
Extending ergmito
George G. Vega Yon
2023-06-13
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
#> ergmito log-likelihood function
#> Number of networks: 5
#> Model: fivenets ~ edges + ttriad
#> Available elements by using the $ operator:
#> loglik: function (params, ..., as_prob = FALSE, total = TRUE) grad : function (params, ...)
# Printing the log-likelihood function
model_object$loglik
#> function (params, ..., as_prob = FALSE, total = TRUE)
#> {
#> if (total) {
#> ans <- sum(exact_loglik(ergmito_ptr, params = params,
#> ..., as_prob = as_prob))
#> if (!as_prob && !is.finite(ans))
#> return(-.Machine$double.xmax * 1e-100)
#> else return(ans)
#> }
#> else {
#> exact_loglik(ergmito_ptr, params = params, ..., as_prob = as_prob)
#> }
#> }
#> <bytecode: 0x55ad78f358a8>
#> <environment: 0x55ad82b1fd60>
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)
#> List of 5
#> $ : num [1:40] 66 12 168 72 48 196 120 1 88 1 ...
#> $ : num [1:40] 66 12 168 72 48 196 120 1 88 1 ...
#> $ : num [1:40] 66 12 168 72 48 196 120 1 88 1 ...
#> $ : num [1:40] 66 12 168 72 48 196 120 1 88 1 ...
#> $ : num [1:40] 66 12 168 72 48 196 120 1 88 1 ...
# The matrices of the sufficient statistics
str(model_object$stats_statmat)
#> List of 5
#> $ : num [1:40, 1:2] 2 1 4 8 7 3 6 12 6 0 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:2] "edges" "ttriple"
#> $ : num [1:40, 1:2] 2 1 4 8 7 3 6 12 6 0 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:2] "edges" "ttriple"
#> $ : num [1:40, 1:2] 2 1 4 8 7 3 6 12 6 0 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:2] "edges" "ttriple"
#> $ : num [1:40, 1:2] 2 1 4 8 7 3 6 12 6 0 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:2] "edges" "ttriple"
#> $ : num [1:40, 1:2] 2 1 4 8 7 3 6 12 6 0 ...
#> ..- attr(*, "dimnames")=List of 2
#> .. ..$ : NULL
#> .. ..$ : chr [1:2] "edges" "ttriple"
# The target statistic
model_object$target_stats
#> edges ttriple
#> [1,] 2 0
#> [2,] 7 4
#> [3,] 4 0
#> [4,] 5 2
#> [5,] 2 0
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))
#> [1] -38.24283
# 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:
summary(ans)
#>
#> Iterations = 10:5000
#> Thinning interval = 10
#> Number of chains = 1
#> Sample size per chain = 500
#>
#> 1. Empirical mean and standard deviation for each variable,
#> plus standard error of the mean:
#>
#> Mean SD Naive SE Time-series SE
#> par1 -1.614 0.4896 0.02190 0.06164
#> par2 1.458 0.5972 0.02671 0.06958
#>
#> 2. Quantiles for each variable:
#>
#> 2.5% 25% 50% 75% 97.5%
#> par1 -2.5556 -1.989 -1.586 -1.246 -0.7719
#> par2 0.2576 1.034 1.439 1.893 2.6050
Finally, we can compare this result with what we obtain from the
ergmito function
summary(ergmito(fivenets ~ edges + nodematch("female")))
#> Warning: The observed statistics (target.statistics) are near or at the
#> boundary of its support, i.e. the Maximum Likelihood Estimates maynot exist or
#> be hard to be estimated. In particular, the statistic(s) "edges",
#> "nodematch.female".
#>
#> ERGMito estimates (MLE)
#> This model includes 5 networks.
#>
#> formula:
#> fivenets ~ edges + nodematch("female")
#>
#> Estimate Std. Error z value Pr(>|z|)
#> edges -1.70475 0.54356 -3.1363 0.001711 **
#> nodematch.female 1.58697 0.64305 2.4679 0.013592 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> AIC: 73.34109 BIC: 77.52978 (Smaller is better.)
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.