knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
This is a short tutorial to show how to employ the
ghypernet R package for the estimation of gHypEG.
In the following, we use Zachary's Karate Club dataset as baseline.
The data is provided within this package and can be simply loaded.
It corresponds to a weighted adjacency matrix whose entries report the number of interactions between individuals.
The graph specified by the adjacency matrix is undirected and does not have self-loops.
For this reason, we specify the two parameters
library(ghypernet) data("adj_karate") directed <- FALSE selfloops <- FALSE print(adj_karate[1:10,1:10])
Fitting a gHypEG Model
GHypEG Models can be fitted using the general function
However, we provide some default functions to fit specific models.
These functions take an adjacency matrix (not sparse) as main argument, or an
igraph graph object.
They further allow to specify whether the model should be directed or not, and whether it should have selfloops or not.
regularm allows to fit the simplest possible model that can be formulated with gHypEGs: gnp-like graph where only the average degree is preserved.
(regularmodel <- regularm(graph = adj_karate, directed = directed, selfloops = selfloops))
scm allows to fit the soft-configuration model, where all degrees are preserved.
(confmodel <- scm(graph = adj_karate, directed = directed, selfloops = selfloops))
From the models, we can generate random realisations that can be used, e.g., as null-models to test hypothesis about the data.
This can be achieved by means of the function
rghype, specifying the number of realisations to take.
random_graphs_rm <- rghype(nsamples = 10, model = regularmodel, seed = 123) random_graphs_scm <- rghype(nsamples = 10, model = confmodel, seed = 123)
Finally, we can perform model selection and hypothesis testing by comparing statistics of two models.
We can follow two different approaches, comparing information scores such as AIC, or performing a likelihood ratio test.
One simple question we can ask, for example, is whether the degree sequence of the Karate Club can be simply explained by a regular model, or there is the need to use a configuration model.
AIC scores can be obtained using the standard R function
Likelihood-ratio tests are performed using the function
However, we provide some wrappers for the function
lr.test for some commonly done tests like the one above.
Comparing AIC scores for the regular model and the configuration model yields the following result:
AIC(regularmodel) AIC(confmodel) # difference in AICs, high value means more complex model is better AIC(regularmodel) - AIC(confmodel)
This shows that the there is strong evidence for the employing the configuration model, because the degree sequence of the graph strongly deviates from the regular one.
To increase confidence in this result, we can compare the result above with that obtained from random data generated from the regular model.
# Generate regular models and configuration models for random graphs regularmodels <- lapply(X = random_graphs_rm, FUN = regularm, directed=directed, selfloops=selfloops) confmodels <- lapply(X = random_graphs_rm, FUN = scm, directed=directed, selfloops=selfloops) # Compute AICs AIC_regularmodels <- sapply(X = regularmodels,FUN = AIC) AIC_confmodels <- sapply(X = confmodels,FUN = AIC) # differences in AIC, high value means more complex model is better AIC_regularmodels - AIC_confmodels
As can be seen by the experiment above, in the case of random graphs generated from the regular model, the AIC of the configuration model is higher than that of the regular model, confirming that the added complexity is not justified.
The (empirical) likelihood-ratio test gives a similar result, with the benefit of providing a p-value.
conf.test provides a simple way to perform the test without the need of computing the models.
conf.test(graph = adj_karate, directed = directed, selfloops = selfloops, seed=123)
Similarly to what done above, we can perform the test using the random graphs. Now we expect large p-values.
tests <- lapply(X = random_graphs_rm, FUN = conf.test, directed = directed, selfloops = selfloops, seed=123) sapply(X = tests, FUN = function(x) x$p.value)
Block-Constrained Configuration Model
The next model that can be estimated is the block constrained configuration model.
The Karate Club has a well-known partitioning into two communities that can be loaded from the package.
We fit a bccm using the
bccm function, specifying the vertex labels.
data("vertexlabels") (blockmodel <- bccm(adj = adj_karate, labels = vertexlabels, directed = directed, selfloops = selfloops)) print(blockmodel$blockOmega)
By default, the function fits a 'full' block structure, where every parameter for in- out-blocks relations are different.
In this case this corresponds to three parameters: one for block 1, one for block 2, one for the edges between the two blocks.
However, in some cases we see that the parameters of in-blocks relations are similar to each other, as can be seen above looking at the diagonal entries of the block matrix.
In this case, a full block structure may overfit the data, while a simpler model could be better.
This can be fitted setting the parameter
This corresponds to a model where there are only two parameters, one for in-block edges, one for out-block edges, irrespectively of the number of blocks.
The result is shown below.
(blockmodel_2 <- bccm(adj = adj_karate, labels = vertexlabels, directed = directed, selfloops = selfloops, homophily = TRUE))
The first test we need to perform is whether either of the model just fitted significantly improve the fit to the data. We do so by performing a likelihood ratio test between the configuration model and the bccms.
lr.test(nullmodel = confmodel, altmodel = blockmodel, seed = 123) lr.test(nullmodel = confmodel, altmodel = blockmodel_2, seed = 123)
Unsurprisingly, the test gives low p-values, confirming the presence of a block structure. Again, we can perform a similar analysis using AIC:
AIC(confmodel) AIC(blockmodel_2) AIC(blockmodel)
From this, we note that while specifying the full bccm gives an improvement in the score compared to the two parameters model, such improvement is relatively small, and appears to not justify the increased complexity. We can verify again this with a likelihood-ratio test:
lr.test(nullmodel = blockmodel_2, altmodel = blockmodel, seed=123)
The result shows that the added complexity of the second model is not justified by the data. We investigate this further by comparing the results with what obtained from random variations.
# First generate random sample from blockmodel random_graphs_bccm2 <- rghype(nsamples = 100, model = blockmodel_2, seed = 123) # Generate the two models for random graphs blockmodels <- lapply(X = random_graphs_bccm2, FUN = bccm, labels = vertexlabels, directed=directed, selfloops=selfloops) blockmodel_2s <- lapply(X = random_graphs_bccm2, FUN = bccm, labels = vertexlabels, directed=directed, selfloops=selfloops, homophily = TRUE) # Compute AICs AIC_blockmodels <- sapply(X = blockmodels, FUN = AIC) AIC_blockmodel_2s <- sapply(X = blockmodel_2s, FUN = AIC) # mean difference in AIC, high value means more complex model is better summary(AIC_blockmodel_2s - AIC_blockmodels)
This confirms the fact that the asymmetry in the two blocks can ascribed to random variations.
Goodness of Fit
Finally, our framework allows to evaluate the goodness-of-fit of model to data.
Similarly to a multinomial goodness-of-fit, we perform a test of the chosen against the maximally complex model that can be formulated.
This model provides the baseline against which comparing simpler models in terms of their goodness-of-fit.
The 'full model' is specified such that the observed graph is the expected one.
It is fitted using the function
ghype with the flag
The likelihood ratio test can be performed either manually or using the function
fullmodel <- ghype(graph = adj_karate, directed = directed, selfloops = selfloops, unbiased = FALSE) lr.test(nullmodel = blockmodel_2, altmodel = fullmodel, seed = 123) gof.test(model = blockmodel_2, seed = 123)
The reason for this result is the fact that, although the bccm is a good fit for the empirical graph, the latter is characterised by some particular structures that are not encoded by the bccm. For example, the empirical graph is characterised by few 'bridges' between the two communities, managed by relatively low degree vertices. Instead, the bccm assumes that all vertices connect weakly the two communities, with an amount of edges proportional to the degree.