MoNAn is the software implementation of the statistical model outlined in:
Block, P., Stadtfeld, C., & Robins, G. (2022). A statistical model for the analysis of mobility tables as weighted networks with an application to faculty hiring networks. Social Networks, 68, 264-278.
AND
Hoffman, M., & Block, P. (2026). Models for emergent structures in mobility: specification and individual-level interpretation. Journal of the Royal Statistical Society Series A: Statistics in Society, qnag057.
Link to the study in Social Networks.
Pre-print with minor differences to the journal version.
Pre-print with minor differences to the journal version.
The model allows for the analysis of emergent structures (known from network analysis) in mobility tables, alongside exogenous predictors. It combines features from classical log-linear models for mobility tables and Exponential Random Graph Models (ERGMs). In the absence of emergent structures, the models reduces to the former, when ignoring characteristics of individuals, it is an ERGM for weighted networks.
Announcements about workshops etc. can be found here.
The MoNAn Manual is available here on github in the manual folder or on SocArXiv. Note that this is the manual for version 1 of MoNAn. This version still represents the functionality and use of the user-side of MoNAn, but shall be replaced soon by a manual for version 2.
This is version 2.1.1 of MoNAn. Compared to version 1, the simulation runs in C++ (estimating the example below takes less than a minute instead of several minutes), several effects have been corrected, and estimation can also use Geyer-Thompson maximum likelihood; see NEWS.md for all changes. Code written for MoNAn 1.0.0 or later works with version 2, except for a few removed effects (see NEWS.md), but results of models with corrected effects have minimal change.
We cannot guarantee that the package and the documentation is free of bugs or errors; further, you might not be able to do what you want. In that case, or if you are unsure please write the package maintainer under his institutional email address.
Finally: Programming Version 2 of MoNAn extensively used Claude Code, which produced the code for the majority of the C++ implementation of the MCMC. The C++ implementation is checked in all steps and effects that it exactly reproduces the human-coded R-based implementation. Additionally, Claude Code was used for optimization; similarly, output from the machine was checked against previous human-coded implementation.
From mid Oct, you can install MoNAn v2.1.1 from CRAN (earlier version 2.0.0):
install.packages("MoNAn")The most recent version (2.1.1) can be installed from GitHub:
# install.packages("remotes")
remotes::install_github("stocnet/MoNAn")MoNAn contains C++ code, so installing it from GitHub requires a C++ compiler (on Windows Rtools, on macOS the Xcode Command Line Tools).
In this section, we outline a simple example with synthetic data stored in the MoNAn package.
library(MoNAn)The example case uses synthetic data that represents mobility of 742 individuals between 17 organisations. The artificial data includes an edgelist, i.e., a list of origins (column 1) and destinations (column 2) of all individuals sorted by origin.
mobilityEdgelist[1:10,]
#> [,1] [,2]
#> [1,] 1 1
#> [2,] 1 1
#> [3,] 1 2
#> [4,] 1 2
#> [5,] 1 2
#> [6,] 1 2
#> [7,] 1 3
#> [8,] 1 3
#> [9,] 1 3
#> [10,] 1 3The data further includes (artificial) characteristics of individuals, which we will call sex. However, using continuous covariates is also possible, although for some effects this would not be meaningful.
indSex[1:10]
#> [1] 1 1 0 0 0 0 1 0 1 0Characteristics of the organisations are the region in which they are located (binary, northern/southern island) and the organisations’ size that represents a composite measure of the number of workers, assets, and revenue.
orgRegion
#> [1] 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0
orgSize
#> [1] 3.202903 5.178353 7.592492 7.368718 9.098950 7.977469 7.751859
#> [8] 10.454544 7.811348 4.958048 7.740106 6.112728 6.640621 2.850886
#> [15] 8.955254 7.597493 9.313139First, we create MoNAn data objects from the introduced data files. A necessary part is that we define and name the nodesets, i.e., we name who (“people”) is mobile between what (“organisations”).
# extract number of individuals and organisations from the mobility data
N_ind <- nrow(mobilityEdgelist)
N_org <- length(unique(as.numeric(mobilityEdgelist)))
# create monan objects
people <- monanEdges(N_ind)
organisations <- monanNodes(N_org)
transfers <- monanDependent(mobilityEdgelist,
nodes = "organisations",
edges = "people")
sameRegion <- outer(orgRegion, orgRegion, "==") * 1
sameRegion <- dyadicCovar(sameRegion, nodeSet = c("organisations", "organisations"))
region <- monadicCovar(orgRegion, nodeSet = "organisations")
size <- monadicCovar(orgSize, nodeSet = "organisations")
sex <- monadicCovar(indSex, nodeSet = "people", addSame = F, addSim = F)We combine the data objects into the process state, i.e., a MoNAn object that stores all information about the data that will be used in the estimation later. We must include one outcome variable (transfers) and two nodeSets (people, organisations).
myState <- monanDataCreate(transfers,
people,
organisations,
sameRegion,
region,
size,
sex)
# inspect the created object
myState
#> dependent variable: transfers with 742 people mobile between 17 organisations
#>
#> covariates of organisations
#> cov. name range mean
#> sameRegion 0-1 0.52
#> region 0-1 0.41
#> size 2.85-10.45 7.09
#>
#> covariates of people
#> cov. name range mean
#> sex 0-1 0.53The predictors in the model are called “Effects”. To specify a model, first an empty effects object is created that only uses the state, and then effects are added one-by-one. Adding an effect requires the effects name and any additional parameters it needs.
# create an effects object
myEffects <- createEffects(myState)
myEffects <- addEffect(myEffects, loops)
myEffects <- addEffect(myEffects, concentration_AC, alpha = 4)
myEffects <- addEffect(myEffects, reciprocity_AC, alpha = 4)
myEffects <- addEffect(myEffects, dyadic_covariate, attribute.index = "sameRegion")
myEffects <- addEffect(myEffects, alter_covariate, attribute.index = "size")
myEffects <- addEffect(myEffects, resource_covar_to_node_covar,
attribute.index = "region",
resource.attribute.index = "sex")
myEffects <- addEffect(myEffects, loops_resource_covar, resource.attribute.index = "sex")
# There is also a simpler way using pipes (|>) and using the more intuitive
# node.attribute & edge.attribute instead of the older
# attribute.index & resource.attribute.index:
myEffects <- createEffects(myState) |>
addEffect(loops) |>
addEffect(concentration_AC, alpha = 4) |>
addEffect(reciprocity_AC, alpha = 4) |>
addEffect(dyadic_covariate, node.attribute = "sameRegion") |>
addEffect(alter_covariate, node.attribute = "size") |>
addEffect(resource_covar_to_node_covar,
node.attribute = "region",
edge.attribute = "sex") |>
addEffect(loops_resource_covar, edge.attribute = "sex")
# inspect the created object
myEffects
#> Effects
#> effect name cov. organisations cov. people parameter
#> loops - - -
#> concentration_AC - - 4
#> reciprocity_AC - - 4
#> dyadic_covariate sameRegion - -
#> alter_covariate size - -
#> resource_covar_to_node_covar region sex -
#> loops_resource_covar - sex -The estimation algorithm requires multiple parameters that guide the simulation of the chain. Most values are generated by default from the state and the effects object.
myAlg <- monanAlgorithmCreate(myState, myEffects, nsubN2 = 3,
multinomialProposal = FALSE)Now we can estimate the model. The first two lines indicate the data (state), effects, and algorithm. The third line specifies the initial estimates, where the previously obtained pseudo-likelihood estimates can be used. The remaining lines define whether parallel computing should be used, how much console output is desired during the estimation, and whether the simulations from phase 3 should be stored.
Running the model takes less than a minute for this data.
myResDN <- monanEstimate(
myState, myEffects,
myAlg,
initialParameters = NULL,
parallel = TRUE, cpus = 4,
verbose = TRUE,
returnDeps = TRUE,
fish = FALSE
)In case pseudo-likelihood estimates have been obtained previously, this can be specified by
myResDN <- monanEstimate(
myState, myEffects, myAlg,
initialParameters = initEst,
parallel = TRUE, cpus = 4,
verbose = TRUE,
returnDeps = TRUE,
fish = FALSE
)Check convergence to see whether the results are reliable. In case the maximum convergence ratio is above 0.1 (or 0.2 for less precise estimates), another run is necessary.
max(abs(myResDN$convergenceStatistics))
#> [1] 0.076954If convergence is too high, update algorithm, re-run estimation with previous results as starting values and check convergence:
# estimate mobility network model again based on previous results to improve convergence
# with an adjusted algorithm
myAlg <- monanAlgorithmCreate(myState, myEffects, multinomialProposal = TRUE,
initialIterationsN2 = 200, nsubN2 = 1, initGain = 0.02, iterationsN3 = 1000)
# for users of other stocnet packages, we can also use monan07 to run an estimation
# (it is an alias for estimateMobilityNetwork)
myResDN <- monan07(
myState, myEffects, myAlg,
prevAns = myResDN,
parallel = TRUE, cpus = 4,
verbose = TRUE,
returnDeps = TRUE,
fish = FALSE
)# check convergence
max(abs(myResDN$convergenceStatistics))
#> [1] 0.076954In case convergence is still poor, updating the algorithm might be necessary. Otherwise, we can view the results. The first column is the estimate, the second column the standard error, and the third column the convergence ratio. All values in the final column should be below 0.1 (see above).
myResDN
#> Results
#> Effects Estimates StandardErrors
#> 1 loops 2.94634644 0.16919809
#> 2 concentration_AC 4 1.40425495 0.20763966
#> 3 reciprocity_AC 4 0.98362353 0.17304975
#> 4 dyadic_covariate sameRegion 0.54538799 0.10733080
#> 5 alter_covariate size 0.01433735 0.01807817
#> 6 resource_covar_to_node_covar region sex -0.64685593 0.15982754
#> 7 loops_resource_covar sex -0.37458067 0.21207675
#> Convergence
#> 1 0.01310609
#> 2 0.01702124
#> 3 0.03168851
#> 4 0.06918573
#> 5 -0.07695400
#> 6 0.01815871
#> 7 -0.03263542The following two diagnostics indicate the extent to which the chain mixes (i.e., whether the thinning was chosen appropriately). The autoCorrelationTest indicates the degree to which the values of the dependent variable of consecutive draws from the chain in phase 3 are correlated. Here lower values are better. Values above 0.5 are very problematic and indicate that a higher thinning is needed.
autoCorrelationTest(myResDN)
#> [1] 0.09969808The output of extractTraces indicates the correlation of statistics between subsequent draws from the chain in phase 3. The plot should show data points randomly scattered around the target line, as shown below. If patterns in the traces are discernible, a higher thinning is needed.
traces <- extractTraces(myResDN, myEffects)
par(mfrow = c(1,2))
plot(traces)
Based on an estimated model, a score-type test is available that
shows whether statistics representing non-included effects are well
represented. If this is not the case, it is likely that including them
will result in significant estimates. The test is conditional on the
statistics of the effects in the model; it requires an estimation with
returnDeps = TRUE.
myEffects2 <- createEffects(myState) |>
addEffect(transitivity_AC, alpha = 8)
test_ME.2 <- scoreTest(myResDN, myEffects2)
test_ME.2
#> Results
#> Effects tStatistic pValue R2
#> 1 transitivity_AC 8 1.49 0.1362 0.979
#>
#> Score test conditional on the model's effects; small p-values: including the effect
#> is likely to result in a significant estimate.
#> R2: share of the variance of the effect's statistic explained by the model's effects.The p-value is above 0.05, so the test gives no clear evidence that transitive clustering is missing from the model, although the positive t-statistic points in that direction. R2 shows that the transitivity statistic is largely explained by the effects already in the model (in particular concentration and reciprocity), which makes its parameter harder to estimate.
Akin to ERGMs, goodness of fit testing is available to see whether auxiliary statistics are well captured by the model.
myGofIndegree <- monanGOF(ans = myResDN, gofFunction = getIndegree, lvls = 1:70)
plot(myGofIndegree, lvls = 20:70)
myGofTieWeight <- monanGOF(ans = myResDN, gofFunction = getTieWeights, lvls = 1:20)
plot(myGofTieWeight, lvls = 1:20)
The package also provides the possibility to exemplarily simulate mobility networks based on the data and the specified effects and parameters.
mySimDN <- monanSimulate(
myState,
myEffects,
parameters = c(2, 1, 1.5, 0.5, 0.1, -1, -0.5),
allowLoops = TRUE,
burnin = 45000,
thinning = 15000,
nSimulations = 10
)
mySimDN[[1]]