Generation and applications of simulated datasets to integrate social network and demographic analyses

Abstract Social networks are tied to population dynamics; interactions are driven by population density and demographic structure, while social relationships can be key determinants of survival and reproductive success. However, difficulties integrating models used in demography and network analysis have limited research at this interface. We introduce the R package genNetDem for simulating integrated network–demographic datasets. It can be used to create longitudinal social network and/or capture–recapture datasets with known properties. It incorporates the ability to generate populations and their social networks, generate grouping events using these networks, simulate social network effects on individual survival, and flexibly sample these longitudinal datasets of social associations. By generating co‐capture data with known statistical relationships, it provides functionality for methodological research. We demonstrate its use with case studies testing how imputation and sampling design influence the success of adding network traits to conventional Cormack–Jolly–Seber (CJS) models. We show that incorporating social network effects into CJS models generates qualitatively accurate results, but with downward‐biased parameter estimates when network position influences survival. Biases are greater when fewer interactions are sampled or fewer individuals observed in each interaction. While our results indicate the potential of incorporating social effects within demographic models, they show that imputing missing network measures alone is insufficient to accurately estimate social effects on survival, pointing to the importance of incorporating network imputation approaches. genNetDem provides a flexible tool to aid these methodological advancements and help researchers testing other sampling considerations in social network studies.


| INTRODUC TI ON
Network analysis has revolutionized animal social behavior research by quantifying how dyadic social interactions and relationships are nested in wider group-and population-level social structures (Krause et al., 2014;Pinter-Wollman et al., 2013). Network studies in behavioral ecology have often focused on how the position of an individual within its social network influences its fitness, either via reproductive success (Formica et al., 2012;Oh & Badyaev, 2010) or survival (Blumstein et al., 2018;Ellis et al., 2017;Stanton & Mann, 2012).
Quantifying direct links between social network position and fitness can help us understand how selection acts on social behavioral traits. Furthermore, determining how social behavior is linked to survival can identify demographic consequences of interactions and associations (Clements et al., 2022), which can help develop more realistic models for how social species respond to population declines or environmental change (Snijders et al., 2017). However, while there is growing interest in linking animal social networks with demography (Shizuka & Johnson, 2020), there remain many methodological challenges.
Currently, most studies that link network position and fitness use known fate approaches such as generalized linear models (e.g., Blumstein et al., 2018) or Cox proportional-hazards models (e.g., Ellis et al., 2017). However, application of these approaches is limited in many wild populations where individuals that are alive are not necessarily detected. In these cases, survival is most commonly estimated using hidden Markov models (HMMs; McClintock et al., 2020) that can simultaneously estimate survival and probabilities of capture (Gimenez et al., 2012;Pradel, 2005). These models also have potential as tools in animal social network analysis (Clements et al., 2022;Fisher et al., 2017), especially when not all associations are detected.
However, it is challenging to provide universal guidance on the applicability of these approaches given the diversity of animal social systems and sampling designs used to study them.
Here, we introduce the R package genNetDem to simulate cocapture datasets. We define a co-capture dataset as one in which capture-recapture data also provide information on social structure, such as when individuals are caught or observed in groups (see also . The package generates integrated longitudinal social network and capture-recapture datasets with known statistical relationships. This provides functionality for methodological research, power analyses, and sampling design. Here, we present an overview the package, outline effective workflows and describe key functions. We then provide two case studies to demonstrate its use.
Finally, we identify key next steps in merging social network and demographic analyses, and discuss the role of genNetDem in these.

| genNetD em OVERVIE W
genNetDem is a set of R (R Core Team, 2021) functions that generate longitudinal social network and/or capture-recapture datasets with known underlying properties. Functionality can be split into four broad groups: (a) population features; (b) survival features; (c) social network features; and (d) observation features. The package is modular meaning-specific components can be used in isolation or usergenerated code can be integrated to extend functionality to different ecological or social contexts. Here, we provide an idea of potential workflows when using genNetDem including a detailed example ( Figure 1) and an overview of key functions (Table 1; with more detail provided in the Appendices S1 and S2). genNetDem is available on GitHub (https://github.com/NETDE M-proje ct/genNe tDem).

| genNetD em WORKFLOW
While genNetDem is designed to be modular so that individual components can be adjusted to perform a range of tasks, many of the functions fit well within specific workflows. We illustrate one such common workflow ( Figure 1), but various other applications are demonstrated in package vignettes. The workflow illustrated here generates a population with a known, underlying social structure and then simulates grouping events (or associations) using this underlying social structure alongside demographic change, sampling from the grouping events to simulate an observation process.

Population generation:
genNetDem provides functionality to simulate a population of a given size that can be subdivided into a prespecified number of (underlying) social groups distributed in 2D space.

Generation of trait data:
genNetDem can also be used to simulate trait data for individuals in the population with considerable flexibility in the types of traits that could be included. It is also possible to use existing biological data or external methods of simulating trait data if preferred as long as the datasets are then formatted in an equivalent manner.

Generate social network:
A key feature of genNetDem is a generative model of underlying social network structures that uses provided information on the presence of social groups, the spatial structure of the population, and traits of individuals within it using an adapted stochastic block model (Lee & Wilkinson, 2019). We use social group to refer to the assignment of individuals to prespecified groups when populations are generated, and spatial structure as any additional effects attributed to the distribution of these groups in 2D space. While using this inbuilt functionality is appealing due to the required inputs and outputs being adapted for other parts of the package, it is equally possible to use other tools to simulate the underlying social network structure. For example, users may want to employ standard generative models (e.g., erdos-renyi random graphs and small-world networks) or to take advantage of the growing availability of more advanced and highly flexible generative models for networks. One example is the STRAND R package (Ross et al., 2022) that combines features from the social relations model alongside the stochastic block model.

Simulate interactions:
It is then possible to use genNetDem to simulate social interactions using this underlying social structure. These interactions/ events can incorporate dyadic or nondyadic interactions; hence, our use of grouping events to describe these (higher-order) interactions generated from an underlying dyadic network of social relationships.

Simulate population processes:
genNetDem additionally provides functionality to simulate survival and recruitment to incorporate population dynamics. Survival can be simulated as a function of an individual's social interactions and nonsocial traits enabling genNetDem to provide a powerful tool to better understand links between social behavior and population processes. Currently, recruitment is strongly density-dependent as a tool to maintain an (approximately) constant population size.

Simulate an observation process:
Finally, genNetDem also provides tools to simulate a capture and observation process based on the simulated grouping events (interactions) such that it is inherently influenced by the underlying social structure. These samples can be used to generate co-capture datasets to test the power to detect social network effects on demographic rates (as illustrated in the case studies) or other research questions of interest.
In a typical workflow, these stages can be linked together to generate longitudinal datasets. For example, in Figure 1, we repeat steps 4, 5, and 6 to generate a co-capture dataset that provides a window into how social network structure and demographic rates are linked in our simulated population. Adapted versions of this workflow are used for the case studies below. F I G U R E 1 An example workflow for using genNetDem to simulate integrated network-demographic datasets. This is a simplified version of the approach used in the case studies with the gray box capturing a demographic timestep and step 4 involving one or more behavioral timesteps. Note that while in our case studies, we fix a relationship between social network position and survival prior to the repeated steps, this relationship in the covariates_survival() function could vary between timesteps if desired, hence its inclusion within the repeated steps box. Note also that multiple set of interactions can be generated prior to the simulation of population processes. The modular nature of the functions in the package means that different parts of this workflow can also be used independently. Further usage examples are provided in the package vignettes.

| genNetD em FUN C TIONS
We provide a description of key functions here and summary in Table 1, with more technical details on functions provided in the Appendix S2.

| Population features
The population features provide capability to simulate a population and generate data about individuals in it. There is then functionality that simulates population dynamics based on individual survival probabilities (see Section 4.2) and stochastic recruitment that maintains an approximately stable population size when employed.
The population_generation_basic() function generates data for a group-structured population distributed uniformly in 2D space. The function takes two arguments: n defines the population size and ng the number of groups in the population. When n = ng, individuals are distributed uniformly at random across the defined coordinates.
When n > ng, groups are distributed uniformly at random across the same coordinates with individuals in the same group sharing the same spatial location. Currently, simulated population size is independent of the extent of the area it occupies. Therefore, population density will increase with population size and impact spatial effects on social network structure. This does not represent a problem except when users want to compare the social structures of populations of different sizes. Group membership is currently fixed once an individual is recruited into the population (although future versions are likely to allow more flexibility in group membership). The indiv_info_gen() and TA B L E 1 An overview of the main functions provided by genNetDem alongside information on their main inputs and outputs. Note that where multiple similar functions exist, we include one example in the table, but alternatives are detailed in the main text and Appendix S2.

| Survival features
The covariates_survival() function allows survival probabilities to be calculated for each individual based on individual traits and the position of an individual within a population social network (this could be any network provided to the function; the underlying social networks, simulated interaction network or a separate user-specified network).
Individual traits specified in the dataframe generated by indiv_info_ gen() can be used as covariates. There is also considerable flexibility in which measures of network position can be included as covariates; both the function and R package used can be specified within the function, with functionality for most common packages (e.g., sna: Butts, 2014;igraph: Csardi & Nepusz, 2006;tnet: Opsahl, 2009) incorporated. It is also possible to simulate network covariance in survival whereby closely connected individuals have either more or less similar survival probabilities than expected by chance (this uses an approximation of the underlying network that is positive definite as a covariance matrix). Note that it may also be possible for some network covariance in survival probabilities to arise without this being encoded directly, for example if survival is positively associated with centrality and more central individuals tend to be more connected with each other. Currently, covariates_survival() simulates independent (additive) effects of traits, meaning that, while the effects of multiple traits can be incorporated together, there is no functionality to capture interactions among variables (e.g., network position having different effects in males than females). The simpler basic_survival() function generates population-level survival probabilities in the absence of covariates.

| Network features
There are two core functionalities of the network features: to generate underlying social networks for the population and to generate grouping events (interactions/associations) based on these networks. There are also two network_checker() functions that quantify and visualize how well social networks derived from grouping events match the underlying network used to generate them.
The network_generation_covariates() function generates an underlying network structure based on social group membership (as defined when generating the population), spatial locations and individual traits. Figure 2 shows examples of networks generated. Current functionality is focused on how these traits impact the probability of forming social connections within and between groups separately, thus employing an adapted stochastic block model (Lee & Wilkinson, 2019). Block membership is defined based on the assignment of individuals to prespecified social groups. It is additionally possible for between block edge probabilities to be further modified by the spatial distance between groups ( Figure S13; the spatial structure -implemented by multiplying baseline values for between-block edge probabilities and weights by 1 d_eff distance , where d _ eff is a user specified effect and distance is the distance between groups). Therefore, genNetDem is not directly designed to incorporate some known social processes such as triadic closure or assortativity, for example females being more closely connected to other females, although it could be possible to use group membership and no spatial structure to approximate these effects (and it is also important to note that assortativity or triadic closure can also arise [indirectly] as an emergent property of the selected generative model). It is also currently not possible for interaction effects to be coded directly (e.g., if size effects on connectivity were different for males and females). Edge probabilities and edge weights are modeled independently to allow variables to explain variation in one or both of them. Edge weights are parameterized by fitting a beta distribution to a provided mean and variance, generating edge weights between 0 and 1 in the underlying network. There is also a simpler network_generation_basic() function that uses the same generative model without covariates.
When modeling longitudinal network data, the individual social network positions could vary from relatively stable to highly dynamic  (Butts, 2014) to conduct a matrix regression between the two networks to test the association between edge weights in each (see Appendix S2 for more detail).

| Observation features
The

| C A S E S TUD IE S
We use two complementary case studies to illustrate the use of genNetDem. In the first, we test how our ability to estimate the relationship between network position and survival depends on sampling effort; whether local or global centrality affects survival; and network dynamics. We compare the performance of cross-sectional versus longitudinal imputation of the network position of nondetected individuals and explore the importance of network covariance in survival probabilities. In the second, we demonstrate how a researcher could use genNetDem to compare sampling designs. We test how the power to estimate relationship between network position and survival depends on how sampling effort is distributed through time. Our simulation asks the question as to whether it is better to concentrate resources into intensively monitoring more groups in fewer sampling windows or fewer groups in more sampling windows. We examine whether any differences are impacted by the proportion of individuals detected in each sampled group and the structure of the underlying social network.

| Methods common to both case studies
In both case studies, we use genNetDem to simulate survival and social interactions and then sample from them to generate capture histories. Illustrations of the workflows used a provided in Figures S1 and S2. We fit HMMs to estimate survival and capture probabilities using nimble (de Valpine et al., 2017(de Valpine et al., , 2022

| Modeling approach
We fitted CJS models estimating both capture and survival probabilities (Lebreton et al., 1992)  We assumed that all individuals in the population were marked or individually identifiable prior to the start of the study. Captures and/ or observations (which were functionally equivalent as all individuals were marked) took place in all behavioral timesteps (50 in total). Each grouping event had either a 25%, 50%, or 75% percent chance of being detected (parameter varied between simulation runs) with the detection probability of an individual in a detected group fixed at 0.9.

| Simulation structure
In total, we generated 3240 simulated datasets, varying five parameters that influenced network dynamics (one parameter), network effects on survival (three parameters), and sampling (one parameter).
• • Sampling: we varied the probability of sampling (either capturing or observing) a group at each behavioral timestep to be 0.25, 0.5, or 0.75.
An illustration of the workflow used is in Figure S1. For each combination of parameters (162), we ran 20 replicates.

| Data recorded
In addition to the four types of data described in the combined methods, we also recorded the full social network generated from all interactions within each demographic timestep (including those not observed). We estimated the network measure of interest from these full networks and scaled them within each demographic timestep as for measures from partial networks.

| Analysis of simulation results
Prior to the general analysis outlined above, we assessed whether the model had converged using the posterior median and standard deviation of its estimate for the social effect on survival. We used k means clustering to identify groups of simulation runs where the model was unlikely to have converged. We used k = 6 clusters and retained three of six of these clusters based on the elbow method and visual inspection of the output ( Figure S3). This method identified ~2.5% of models had likely not converged.
To compare the success of models that used network measures calculated from the partial versus full network, we calculated the earth mover's distance (

| Results and discussion
Overall, we show it is possible to estimate social effects survival from partial networks, albeit with substantial limitations in power ( Figure 3, Table 2; Table S1). Estimates of social effects on survival were downward biased meaning that statistical power was limited and only stronger social effects on survival are likely to be detected.
Sampling effort was particularly important and interacted with how imputation was conducted in determining how well models converged and biases in parameter estimates when they did. Estimates of other parameters were unaffected.
Previous research has demonstrated that network measures from sampled, partial networks are correlated with those in the full, unobserved network but that these correlations vary depending on the proportion sampled and network measure calculated (Silk et al., 2015;Smith & Moody, 2013). Furthermore, the regression slope is rarely 1:1 indicating values for measures estimated are not perfectly accurate (Silk et al., 2015). This likely explains many of our results showing the difficulty of detecting social effects on survival in the absence of network imputation or the use of measures from independently (and better) sampled social networks.

| Network variable and covariance structure
When we compared models that used network measures from the full and partial networks, we found downward-biased parameter estimates and reduced statistical clarity of results when partial network measures were used (Table 2, Figure 3). These patterns were more striking when survival was related to a global measure of centrality (betweenness) than a local measure of centrality (strength). We found that including positive or negative covariance in survival probabilities related to social network structure had little effect on estimation or power in the contexts simulated ( Figure S5, Tables S5 and S6). These results are promising in suggesting that this may present a more limited issue in this context than often considered (e.g., Croft et al., 2011;Farine & Carter, 2020;Silk et al., 2017). However, the importance of covariance likely depends substantially on network structure and density, so it would be unwise to generalize these patterns without further work focused specifically on this question.

| Sampling effort, imputation approach, and network dynamics
Lower sampling effort was typically associated with both (a) reduced likelihood of model convergence (Table 3; Table S2), and (b) downward-biased parameter estimates (Figure 3). However, the nature of these relationships depended on the imputation approach  Figure S4).
Models were much less likely to converge when sampling effort was low (25% group capture probability), betweenness centrality from partial networks was used as an explanatory variable and cross-sectional imputation was used to infer missing values (Table 3;   Table S2). Even when 50% of groups were sampled in these situations, there was still a reduction in convergence rate. Note that this was apparent regardless of whether betweenness centrality had a positive or no effect on survival probability. Any other changes in the likelihood of model convergence were of much smaller magnitude, but generally occurred when sampling effort was low (and measures from partial networks were used).  (Table 2; Table S1), indicating that posterior distributions had higher variance when cross-sectional imputation was used.
When measures from the partial network were used instead, there was a much more substantial reduction in both parameter estimates and statistical power apparent even for higher sampling efforts ( Figure 3, Table 2). Reductions in parameter estimates were more substantial and remained linear when longitudinal imputation was used, instead flattening out for cross-sectional imputation so that the difference between 25% and 50% of groups being sampled was less than the difference between 50% and 75% ( Figure 3).
However, similarly to the pattern for full network measures, this was not reflected in changes to statistical power which were broadly equivalent for both, indicating a less precise posterior distribution for cross-sectional imputation. These differences between crosssectional and longitudinal imputation changed how EMDs calculated for the differences between posteriors from the full network and partial network model fits depended on sampling effort (Figure 4).
For cross-sectional imputation, EMDs were highest for low sampling effort (p = .25) while for longitudinal imputation, they peaked at intermediate sampling effort (p = .5). However, in general, EMDs were higher for cross-sectional than longitudinal imputation.
Our results show that when social networks are constructed based on the same co-capture data used to estimate survival, even relatively small drops in sampling effort can lead to downward biases in parameter estimates and statistical power. While TA B L E 2 Proportion of simulation runs where 0 falls outside the 89% HDI for different parameter combinations.

| Estimates of other parameters
Estimates for other parameter values were unaffected by social effects on survival, use of measures from full or partial networks or imputation strategy ( Figures S6-S8). survival probability of 0.8 in females. In this case study, there was recruitment into the population over time (i.e., the population stayed roughly constant over each simulation). There was a 10% chance that a surviving individual rewired its underlying social connections after each demographic timestep, and if it did each connection had a 50% chance of changing. Edges were rewired using the same generative model used to create the initial network.

|
The population was initially unmarked. Captures only occurred in the first behavioral timestep of each demographic timestep with 90% of groups sampled and a 0.9 probability of individuals in a sampled group being detected. Sampling design and effort for subsequent observations depended on parameter choice (see below).

| Simulation structure
In total, we generated 2880 simulated datasets, varying five parameters that influenced network structure (one parameter), network effects on survival (two parameters), and sampling effort/design (two parameters).
• Network structure: we varied underlying network structure so that either (a) there was no group structure and moderate spatial structure driving the probability and weight of edges or (b) the population was divided into 20 groups with the probability of a within-group connection of 0.5 and within-group connection weights having a mean of 0.5 (vs. a baseline of 0.2 and 0.25, respectively, for between-group connections prior to adjusting for distance effects).
• Network effects on survival: (a) we varied the network measure that influenced survival to be either strength or betweenness and (b) we varied the effect size to be 0 (no effect), 0.4 (moderate effect), or 0.8 (strong effect).
• An illustration of the workflow used is in Figure S2. For each combination of parameters (144) we ran 20 replicates.

| Model-fitting
Unlike Case Study 1, each CJS model was conditioned on first capture (as individuals were not assumed to have been captured previously).

| Results and discussion
Overall, survival models performed adequately in detecting social effects on survival (Table 4, Figure 5; Tables S7-S9, Figure S9). When we simulated positive effects of network centrality on survival prob-  (Figure 5), and results were more frequently statistically unclear (Table 4). The results here support those from Case Study 1 and the existing literature (Silk et al., 2015;Smith & Moody, 2013) in highlighting that global measures of network position are more susceptible to sampling effects than local measures.

| Sampling design
There was no clear effect of how grouping events were sampled within each demographic timestep on estimates of social effects on survival ( Figure 5). Unsurprisingly, probability of observing individuals within groups did have some effect, with less downward-biased parameter estimates and more statistical power when sampling within groups was more complete ( Figure S9, Tables S7 and S8), as would be expected.
Lower observation success within sampled grouping events leading to reduced model performance is unsurprising as it leads to missing edges in the sampled network, reducing its correlation with the true (unobserved) network. This finding supports related work focused on calculating network measures (e.g., Franks et al., 2010). Franks et al. (2010) also tentatively recommended that more censuses (behavioral timesteps sampled in our case) were preferable than ensuring a high proportion of grouping events sampled in each census for calculating weighted measures of centrality. However, we found no clear evidence that this extended to our survival analysis, where there were only small differences in model performance and no clear overall trend. It should be noted, however, that the simulation architecture differed between the two papers.

| Social structure
Social structure had a small effect on the ability to detect social effects on survival, with some differences in statistical power between the two structures investigated. While there were minimal differences in posterior medians ( Figure 5), results tended to be statistically clearer when there was no underlying group structure than when the population was divided into 20 groups (Table 4).
Previous studies of sampling in social networks have rarely considered the types of modular social structures common for group-living animal populations (Silk, 2018). The slight negative impact of this group-structure on our ability to detect social effects on survival perhaps suggests that the correlation between network measures F I G U R E 5 Impacts of sampling design (within-plot: sets of boxes of the same color), within-group detection probability (columns) and social structure (rows) on Cormack-Jolly-Seber estimates of social effects on survival probability for a range of simulated effect sizes (colors of boxes). Boxplots show the distribution of posterior medians from multiple simulation runs with the solid line the median, boxes the interquartile range and whiskers the full range of values. We illustrate contexts in which a local measure of centrality (strength) and global measure of centrality (betweenness) are used as explanatory variables. The blue-dotted line indicates the accurate parameter estimate when the true effect size is 0.4 (the equivalent line for 0.8 is not illustrated).
calculated in the sampled and full networks is weaker in these types of networks.

| Estimates of other parameters
Estimates of other parameters were largely unaffected by social effects or sampling design. Strong social effects on survival were associated with slightly lower estimates of mean survival probability, but these differences were caused by differences in simulated survival probabilities rather than model performance ( Figures S10-S12).
While limited in scope, these results provide evidence that including social effects on survival in demographic models is unlikely to impact other parameter estimates substantially (see also Clements et al., 2022).

| FUTURE S TEPS
With the two case studies presented, we only scratch the surface of the potential of genNetDem as a methodological tool for animal social network analyses. Below, we highlight some logical next steps for methodological studies on this topic, focusing on the integration of social networks and demography.
First, while we demonstrated the capacity for genNetDem to generate diverse social structures (Figure 2), this was only a partial focus of our results. Animal social systems vary widely, and while optimal sampling strategies are likely to vary with social structure (Clements et al., 2022;Silk, 2018;Sunga et al., 2021), this has remained understudied. Similarly, while we varied network dynamics in our simulations, individual variation in edge probabilities was limited.
Incorporating greater trait-based or individual variation into network position would likely influence conclusions drawn about imputation approaches, for example. The modular design of genNetDem allows it to be integrated with other tools to simulate social network structure (e.g., Ross et al., 2022), which will help tackle these types of challenges more comprehensively in future. Similarly, Clements et al. (2022) included estimation of network structure within a CJS model to improve estimation of social effects on survival. However, the latter approach used a rather basic generative model for the latent network structure that could be improved on or adjusted for researchers working in different contexts.
Consequently, extending these approaches to incorporate more sophisticated social network models as well as to open populations is a key priority.
Third, to keep our case studies accessible, we examined social effects only in CJS models to estimate survival probability. Clements et al., (2022) highlighted the potential value of incorporating social networks within integrated population models (IPMs), where different data sources could also be used to inform network structure itself. However, especially with improvements to imputation of latent network structures, there is also great potential to incorporate network effects within multistate models more generally. Given the central role of social behavior in mediating interactions between infectious disease dynamics and demographic processes (Silk et al., 2019;Silk & Fefferman, 2021), extending multistate models to incorporate social network structure in this way could provide important new insights into wildlife disease ecology, to provide just one example. genNetDem can provide an ideal sandbox to refine these models for application to wild systems.
Finally, we focus here on dyadic social networks; however, many of the social interactions studied are nondyadic and may include higherorder interactions (Battiston et al., 2021;Greening Jr et al., 2015).
While there has been limited focus on higher-order interactions in animal societies (Musciotto et al., 2022), theory suggests they will impact infectious disease transmission and social contagions (Battiston et al., 2021;Iacopini et al., 2022;Noonan & Lambiotte, 2021) among other ecological and evolutionary processes. Therefore, expanding some of the developments here beyond dyadic networks to consider higher-order effects on survival and imputation of hyperedges (social connections between more than two individuals) will likely represent valuable developments. Because it generates GBIs that incorporate interactions/associations between more than two individuals, gen-NetDem is an ideal starting point for methodological research testing higher-order methods in animal societies.

| CON CLUS IONS
We introduce the R package genNetDem as a flexible tool for simulating combined social and demographic datasets. While we focus on the integration of social network and demographic models, the modular design of the package allows it to be an equally powerful tool for generating social network or capture-recapture datasets in their own right. It therefore provides a general tool for researchers interested in testing key methodological considerations in animal social network studies, especially as the field moves towards longitudinal analysis. It also helps researchers wishing to test the power of specific analyses or sampling designs in their own study systems.

ACK N OWLED G M ENTS
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 101023948. We thank Sarah Clements, Qing Zhao, Mitch Weegman and Dave Hodgson for insightful discussions related to this work. We would also like to thank two reviewers for constructive comments that helped improve the paper.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors have no conflicts of interest to declare.