BISoN: A Bayesian framework for inference of social networks

Animal social networks are often constructed from point estimates of edge weights. In many contexts, edge weights are inferred from observational data, and the uncertainty around estimates can be affected by various factors. Though this has been acknowledged in previous work, methods that explicitly quantify uncertainty in edge weights have not yet been widely adopted and remain undeveloped for many common types of data. Furthermore, existing methods are unable to cope with some of the complexities often found in observational data, and do not propagate uncertainty in edge weights to subsequent statistical analyses. We introduce a unified Bayesian framework for modelling social networks based on observational data. This framework, which we call BISoN, can accommodate many common types of observational social data, can capture confounds and model effects at the level of observations and is fully compatible with popular methods used in social network analysis. We show how the framework can be applied to common types of data and how various types of downstream statistical analyses can be performed, including non‐random association tests and regressions on network properties. Our framework opens up the opportunity to test new types of hypotheses, make full use of observational datasets, and increase the reliability of scientific inferences. We have made both an R package and example R scripts available to enable adoption of the framework.

which the data are collected affects both the interpretation of social network analyses and their accuracy (Whitehead, 2008a).
Networks are usually constructed by taking a normalised measure of sociality, such as the proportion of sampling periods each pair spends engaged in a social behaviour (e.g. the Simple Ratio Index [SRI]; Cairns & Schwager, 1987). These normalised measures ensure that pairs that are observed for longer are not erroneously determined to be more social (Whitehead, 2008a). This is important because observing social events can be challenging and uniform sampling over all pairs is not always possible (Croft et al., 2008).
Though normalised measures of social structure will not be biased by sampling time, they will be accompanied by varying levels of certainty. For example, edge weights will be treated as a certain 0.5 for both a case where individuals have been seen together once and apart once, and equally where individuals have been together 100 times and apart 100 times. Existing solutions to this problem aim to either a priori quantify accuracy and power of social networks (Hart et al., 2021;Whitehead, 2008b), filter out highly uncertain data points (James et al., 2009) or include sampling effort as terms in statistical models . These existing solutions can be useful, but also necessarily lose information and may introduce biases into analyses. Though uncertainty is likely present in most measurements (e.g. age, size or perhaps even sex), the abovementioned studies have highlighted that uncertainty due to sampling effort can be orders of magnitude greater than other sources of measurement uncertainty. Any major sources of uncertainty can drastically impact the performance of statistical models, so it is important that uncertainty due to sampling effort is accounted for properly (Wasserman, 2004). Despite this, methods to fully estimate and propagate uncertainty through social network analyses remain largely underdeveloped.
In this paper, we introduce BISoN: a general unified Bayesian framework for modelling social network data. BISoN captures uncertainty in edge weights, making full use of available data; propagates uncertainty through downstream analyses; and can control for social and non-social effects at any level (individual, dyad, group, observation, etc). Any type of social network analysis can be conducted within our framework, including dyadic and nodal regressions, non-random edge weight tests and estimation of structural network properties. The BISoN framework comprises (a) an edge weight model, dependent on data type, that builds edge weights and networks with uncertainty from empirical data, and (b) downstream analyses that use estimated edge weights to take into account uncertainty in the network. This manuscript is intended as a technical overview of BISoN and assumes some familiarity with the main concepts of Bayesian statistics and social network analysis. Readers unfamiliar with the Bayesian statistical workflow and social network analysis are directed to introductory texts such as Kruschke (2015), McElreath (2020), Carrington et al. (2005), and Croft et al. (2008). We will focus on how our framework can be applied to animal systems but the underlying principle can be applied to any type of network analysis where network edges are inferred.
The BISoN framework presents an opportunity to generate reliable, flexible and rich scientific inference in the study of social systems.

| B ISoN-BAYE S IAN INFEREN CE OF SO CIAL NE T WORK S
In the following section, we outline the three edge weight models we have developed for binary, count, and duration data (see Table 1 for definitions of these types of data). For brevity, the models are presented without detailing any specific priors but when fitting these models, priors should always be specified, and the choice of prior should depend on the context and the question in hand (van de Schoot et al., 2021). These edge weight models can be used to generate estimates for network edge weights ij . We describe how to obtain edge weight estimates for binary, count and duration data below. Once estimates for edge weights ij are obtained, they can be used in downstream analyses. We will briefly outline three types of analyses that propagate uncertainty through the analysis. These analyses are: (1) testing for non-random edge weights, (2) dyadic regression and (3) nodal regression. When presenting the downstream analyses, for brevity, we will refer to them in the context of the count data model but these downstream analyses are freely applicable to all types of data. An overview of the BISoN modelling workflow is shown in Figure 1. The models we present here are simple examples that can be extended and refined in many ways. These methods will be useful and reliable, providing they are used appropriately in the context of the data, the scientific question and with appropriate diagnostic tools (Kruschke, 2015;McElreath, 2020 where X ij is the number of presence/absence social events that occurred between i and j, and D ij is the number of sampling periods where it is considered possible for a social event to have been observed between i and j. The exact definition will depend on context, but this may often be the number of periods in which at least one of the individuals of the dyad was seen. The edge weight is therefore equivalent to the probability of i and j engaging in a social event in any given sampling period. This can also be interpreted as an estimate of the proportion of time i and j spend engaging in a social event. Let us refer to this probability (or proportion) as p ij . Without making any additional assumptions beyond that of the standard SRI calculation, an equivalent formulation is Note that the maximum likelihood estimate for p ij in this process The binomial process is the natural process for multiple Bernoulli trials, which makes this formulation a natural extension of point estimates for edge weights. In a Bayesian context, the parameter p ij is treated as a random variable, and inherently captures any uncertainty around its value. A similar formulation was used by (Farine & Strandburg-Peshkin, 2015) to estimate uncertainty over edge weights using a beta conjugate prior over the p ij parameters (as detailed in Fink, 1997).
Social data can sometimes be influenced by non-social effects at the observation level. For example, varying levels of habitat visibility depending on location, variance in observer reliability, or time of day can all affect the recorded observations of social events. These effects cannot be modelled by aggregating social events X ij at the dyad level. To solve this, we propose to model these effects by decomposing the binomial process over the aggregated observations to distinct Bernoulli processes over unaggregated observations, as follows: F I G U R E 1 Schematic overview of the BISoN framework. Different types of observation data in various levels of aggregation are fed into the appropriate edge weight model for the type of data. The edge weight model generates posterior distributions over edge weights, quantifying uncertainty over them. Posterior edge weights can be visualised in a sociogram with uncertainty, or extracted as an edge list with credible intervals. Posterior edge weights can also be extracted for downstream analyses such as regression analyses, non-random edge weight tests and more.
where p (n) ij is the probability of a social event occurring between i and in the nth observation of i and j, and ij is the edge weight between i and j. The ellipsis represents additional terms that can be included (if desired) to model various effects. The term logit p (n) ij can be seen as analogous to the predictor terms in regression models. In the locationdependent visibility example, including an observation-level location effect will account for this visibility effect. When additional effects are included, the edge weight will become an estimate of relative sociality in the presence of non-social effects. We use a linear combination of effects in these examples for simplicity, but the effects could equally be modelled as non-linear functions if deemed necessary in context.
Alternatively dyad-or individual-level effects can also be included to separate out different types of social effect. The binary data model can also be used for group-based (gambit of the group) data (Whitehead & Dufault, 1999). If being used for group data, an additional term in the model can be used to account for non-independence within group observations. See the Supporting Information (S1) for further discussions on these points.

| Edge weight model: Count data
The same concept as above can also be applied to the case of count data, where counts of social events between individuals per sampling period are recorded. First, note that for count data, point estimates of edge weights are often also defined as W ij = X ij ∕ D ij but where X ij is the total number of social events that occurred between i and j, and D ij is the total amount of time a social event could have been observed between i and j.
The edge weight is therefore equivalent to the rate at which i and j engaged in a social event per unit time. Let us refer to this rate as ij . When events occur with a fixed rate, the number of events expected per unit time is described by the Poisson process (Wasserman, 2004). Since we assume that the rate of events characterises edge weights well, then without making any additional assumptions, we can model social events as where ij is the underlying rate parameter. When modelled in a Bayesian framework, this provides a natural way to model the uncertainty around edge weights because the maximum likelihood estimate (Held & Bové, 2014). Again, this model does not allow observation-level effects to be modelled, so we must decompose the aggregated data into the original sampling periods. Fortunately, the sum of Poisson-distributed random variables is also described by a Poisson distribution. This means that the full model can take a similar form to the previous model: where (n) ij is the rate at which social events occur between i and j in the nth observation of i and j, D (n) ij is the length of time of the nth sampling period for the dyad, and ij is again the edge weight between i and j. As above, the ellipsis represents potential additional effects that could be included in the model if necessary. This is discussed in a later section and in the Supporting Information (S2).

| Edge weight model: Duration data
The edge weight model for duration data is slightly more complex than those of the binary and count data models, because it models two different types of data simultaneously: the duration of social events, and the frequency with which they occur. With duration data, the edge weight is again calculated using W ij = X ij ∕ D ij , but where X ij is the amount of time i and j were seen engaged in a social event and D ij is the maximum amount of time i and j could have been seen engaged in a social event. The total amount of time X ij is the sum of times from K ij social events. The edge weight is then modelled as the proportion of time t ij each dyad could engage in a social event that is actually spent engaging in the social event: where ij is the mean event duration for the dyad and ij is the mean rate of events per unit time for the dyad. Inter-arrival times in a Poisson process are distributed according to the exponential distribution, so we use this to model the duration of social events once they have started. The count of events can again be considered to be a Poisson process. Bringing these two assumptions together gives the following model (see Supporting Information S2 for the full derivation): The mean event time ij is modelled inherently in t (n) ij and can be recovered from the fitted model afterwards using ij = ij ∕ ij .

| Accounting for social and non-social effects
By modelling social events at the observation-level rather than the dyad-level, it is possible to account for observation-level effects that may influence the estimated edge weights. An example of this is where the study population might move through various locations with different visibilities, leading to location-dependent missed social events. Assuming that variation in sociality between location is not of direct interest to the researcher, we can include a location effect in the model. An example of how this can be achieved in the count model is where L (n) ij is a location effect corresponding to the location of the social events observed for the dyad ij in the nth observation period. The location effect has an adaptive normal prior centered at zero (often known as a random effect). The location effect is designed to capture the differences in visibility at different locations.
Note that if more direct information such as vegetation density is available that can directly model, for example, visibility, this will always be preferable to weaker proxy variables. In theory, the effect of visibility on edge weights will always be negative but we choose to model it without bounds. This is because the true visibilities in different locations are unknown, so the edge weights cannot be estimated exactly, only proportionally. It is therefore unnecessary to place an upper bound on the visibility effect, making the model easier to fit.
Depending on the assumed causal process that generates observations, it will not always be possible to control for non-social effects in this way, as the social effect of interest may partially be borne out through apparent non-social effects, such as space use. For example individuals may genuinely be more social in some locations than others, and controlling for location would remove the social effect of interest. A discussion on the considerations around modelling non-social effects is included in the Supporting Information (S7).

| DOWNS TRE AM ANALYS IS WITH UN CERTAINT Y
The previous section has shown how uncertainty can be estimated over edge weights. Once these estimates have been obtained, it is then possible to conduct various types of statistical analysis on the network while propagating uncertainty through the entire analysis. This approach will lead to lower 'false positive' error rates because network quantities will have higher variances. Despite this, because the uncertainty captured by BISoN reflects true uncertainty in the data, it will also have a 'false negative' error rate no higher than comparable analyses that do not quantify uncertainty. In this section we will discuss how several common types of social network analysis can be applied to social networks constructed with BISoN edge weight models. In particular, we outline how three specific types of network analysis can be conducted: (1) testing for non-random edge weights, (2) dyadic regression and (3) nodal regression.

| Testing for non-random edge weights
BISoN enables a Bayesian test analogous to the classic non-random association test developed by (Bejder et al., 1998) from a species co-occurrence test developed by (Manly, 1995), and originally proposed as the 'chance sociogram' by (Moreno & Jennings, 1938). At its core, this test asks whether the structure of an observed network is an artefact of sampling alone or if there is genuine variation in edge weights between individuals. If the structure of an observed network is simply due to sampling, then there is no variation in edge weight, and each dyad is equally likely to engage in social events as all other dyads. The original test from (Bejder et al., 1998) is extendable to test for non-random associations within categories (such as sex combination or age difference). However, it requires continuous variables such as age difference to be discretised into categories. Discretising variables is generally not ideal, and removes information, weakens statistical power, and can lose useful insights (Altman & Royston, 2006).
Our proposed test can account for variables such as age without discretisation. This can be achieved by adding variables to the predictor in the following way: where Age is a parameter describing the effects of age on the log rate of social events, and A ij is the age difference for the dyad ij. In addition to modelling continuous variables, this approach to the non-random edge weight testing also makes it possible to compare multiple competing hypotheses of edge structure, model the effects of different dyad-level covariates, and compare hypotheses relating to global network metrics such as network density or clustering coefficients.

| Dyadic regression
A common type of analysis in studies of social networks is to con- where 0 is the intercept parameter, 1 is the slope parameter, u i are parameters accounting for the effect of some edges sharing the same node, is the standard deviation of the residuals, and u is a hyper-

| Nodal regression
The final common type of network analysis we will cover here is nodal regression, where a regression is performed to analyse the relationship between a nodal network metric (such as centrality) and nodal traits (such as age and sex). These analyses are usually used to assess how network position depends on various biological factors, but can also be used where network position is a predictor. Since node metrics are derivative measures of the network, uncertainty in edge weights should ideally propagate through to node metrics, and on to coefficients in regression analyses, giving an accurate estimate of the total uncertainty in inferred parameters. The core of the nodal regression is similar to dyadic regression, taking the form where 0 is the intercept parameter, 1 is the slope parameter, and is the standard deviation of the residuals. M i ( ) denotes the node metric estimates when applied to the edge weights estimated in the top part of the model.
Calculating node metrics within the model may present a practical challenge when using standard model fitting software, as many node metrics cannot be described purely in terms of closed-form equations (Freeman, 1978). In this case, splitting up the model along the dashed line becomes an important step, as it allows the edge weights to be estimated using a standard piece of software, and leaves the regression part of the model to be fit using a sampler that supports numeric estimates of network centralities. As part of the code we have made available, we have written a custom Metropolis-Hastings sampler that maintains the joint distribution of edge weights, rather than sampling edge weights independently.
This ensures that the effect of any structure in the observation data (such as location effects) is maintained and propagated through to the node metric estimates, and subsequently to regression coefficient estimates.

| E X AMPLE: SYNTHE TI C DATA S E T
To demonstrate our framework, we present an example based on a synthetic dataset of social events simulating count data. The where X (n) ij is the count of social events between i and j in their nth observation period, (n) ij is the social event rate from the corresponding period, D (n) ij is the length of the corresponding sampling period (fixed to 1 here for simplicity), ij is the dyad-level edge weight, and L (n) ij is an observation-level effect corresponding to one of the six locations. The The model fit was checked visually by examining the chains and computing the Gelman-Rubin R convergence statistic (Gelman & Rubin, 1992). Once the model fit was verified, the parameter estimates of edge weights ij between dyads were extracted. Figure 2a  The 95% credible interval for the control parameter Control was − 1.27, − 0.36 , with a median of − 0.83, and for the treatment parameter Treatment was 0.35, 1.29 , with a median of 0.85. One of the many benefits of Bayesian inference is that posterior distributions of different parameters can themselves be used gain statistical insights. In this case, the question of real interest is the difference between the two treatment conditions: Treatment − Control . This quantity can simply be calculated by subtracting the posteriors corresponding to the treatment parameter from those corresponding to , F I G U R E 2 (a) Sociogram visualisation of a network generated from uncertain edge weights. The width of edges denotes their weights, and the bands around the edges are proportional to the width of their 95% credible intervals. Low edge weights are arbitrarily removed from this plot (though not the analysis) for visualisation purposes. (b) Posterior distributions of node-level eigenvector centralities for each node in the network. Nodes 1-4 are those assigned to the treatment and nodes 5-8 to the control. the control parameter. We found that the difference in standardised centrality between treatment conditions, Treatment − Control , had a 95% credible interval of 0.99, 2.33 and a median of 1.69. The posterior distribution of the difference is also shown in Figure 2b. This indicates that the treatment condition appears to correlate strongly with a higher centrality.
This brief example outlines just some of the types of analysis that can be carried out in our framework. We have included several additional examples in the code included with this paper.

| E X AMPLE: SOUTHERN RE S IDENT K ILLER WHALE S
To demonstrate BISoN on empirical data, we used a publiclyavailable dataset of near-surface physical contact interactions between 22 southern resident killer whales (Orcinus orca) . We fitted a BISoN count edge model to the data and ran a regression analysis to study the associations between node strength and both age and sex. See Supporting Information (S4) for replicable code of this analysis. We opted for a standard linear model with a Gaussian family distribution, implemented using the bison_brm function in bisonR. From previous work, we expected a negative association between node centrality and age, and we also expected males to be less central than females .
Our findings were in line with this, and we found a negative association of age with centrality, with an estimated age coefficient of −0.57 (95% CI: −0.80, −0.35). We also found that males had a lower centrality than females, with an estimated sex coefficient of −9.99 (95% CI: −15.32, −4.67). We hope this findings serve as a useful empirical example of how BISoN can be used to run regression analyses.

| DISCUSS ION
In this paper, we presented BISoN, a general, highly flexible framework for modelling and analysing social networks. BISoN uses models of data collection processes to quantify uncertainty in edge weights. We see BISoN as a generalised synthesis and extension of a wide range of existing methods. In particular, for binary data, BISoN implements a variant of the beta binomial measurement error model, which has been previously used in animal social network analysis (Farine & Strandburg-Peshkin, 2015;Fink, 1997;Fuller, 1987).
BISoN's mechanism for modelling multiple effects was inspired by generalised linear models, and can be seen as an extension of social relations models (Kenny & La Voie, 1984;McCulloch & Searle, 2004). This concept was introduced to the animal social network literature by (Whitehead & James, 2015) as generalised affilitation indices. The BISoN framework shares some conceptual similarities with the recently developed STRAND framework, but implements measurement models specifically designed for common types of empirical data, and allows propagation of uncertainty to downstream analyses.
One of the main benefits of our framework is its inherent flexibility but this also introduces a number of considerations when using it. Setting priors is an important part of Bayesian modelling, and it is important that they reflect realistic prior expectations of parameter values (van de Schoot et al., 2021). The priors included in the example code are not intended as good defaults, and instead the priors should be determined by careful consideration of the meaning of the parameter, prior expectations of likely values for the parameter, and using prior predictive checks (Conn et al., 2018). It is also important to note that using additional terms and altering the predictor part of the model can change the interpretation of model parameters, which may require the priors to be altered as well.
Another consideration is that network centrality measures based on binarised edges, such as degree and transitivity, have a slightly different interpretation in this framework. This is because, even though a particular dyad may never have been seen interacting in, say 100 samples, it is still possible that they would interact in the the 101st A final consideration of our modelling framework is that, though non-social effects can be modelled in the predictor, this alone cannot remove some non-social effects, and could instead absorb the social effect of interest, depending on the underlying causal system.
We recommend that these effects be treated with caution, and used only with explicit assumptions about the underlying causal system.
These considerations are also discussed in further depth in the Supporting Information (S7) using causal models.
The versatility of our proposed framework makes the approach powerful, but comes with a number of limitations. Though highly flexible, one important constraint on BISoN is that knowledge of sampling effort will be required to properly model the processes that generate observations. On a practical level, BISoN models will usually need to be fit using a Markov chain Monte Carlo (MCMC) sampler. MCMC can be slow and computationally expensive, especially with large datasets. If additional effects are not included, this issue can be partially overcome by aggregating data at the dyadlevel and fitting collapsed versions of the models (see Supporting Information S8). Furthermore, MCMC methods do not guarantee convergence, so additional checks will be required to make sure the sampler has performed correctly (Cowles & Carlin, 1996;Draper, 2008). Together these limitations make the process of model building, fitting and checking a more complex undertaking than simply building networks from point estimates. However, in return for these investments, researchers will have access to a powerful suite of tools that can make maximum use of hard-won data.
The BISoN framework is open-ended and can be extended in many ways. One particularly useful aspect of the Bayesian approach is the ability to treat data with uncertainty. This allows missing data to be modelled within the statistical model without the need for repeated model fitting routines. Furthermore, missing data can be estimated from other variables, which could be especially useful in populations where traits such as age are often estimated or unknown (Little & Rubin, 2002). Another potential direction for the framework is to move beyond static sociality measures and instead explicitly model the dynamics of edge weights, as discussed by (Pinter-Wollman et al., 2014). Dynamic edge weights could be modelled by incorporating autoregressive time series models into the framework (Wei, 2013). This would make it possible to model the general trends in edge weights, how changes in edge weights propagate through the network, and even model flow processes such as disease or information transmission.

| CON CLUS ION
Uncertainty affects all social network analyses, yet until now has been difficult to quantify. The BISoN framework we have introduced is designed to help quantify uncertainty and ensure it is preserved through subsequent analyses to give robust and reliable statistical inference. Our framework is general and widely applicable to almost any type of social relational data. Estimated edge weights from our models can be used with many different types of downstream analysis while propagating uncertainty. We believe the ideas presented in this paper are only the first step in developing powerful, flexible and robust statistical models for social networks. BISoN has the potential to enable us to ask a wide range of nuanced questions, and obtain a deeper insight into the nature of sociality.

AUTH O R CO NTR I B UTI O N S
Jordan Hart conceived of the initial idea behind BISoN, which was further developed with input from Michael Nash Weiss, Daniel Franks and Lauren Brent. The example code was written by Jordan Hart. The manuscript was written by Jordan Hart with input from Michael Nash Weiss, Daniel Franks and Lauren Brent.