The EGA+GNM framework: An integrative approach to modelling behavioural syndromes

Behavioural syndromes refer to correlated suites of behavioural traits exhibiting consistent among‐individual variation, i.e. personality. Factor analysis (FA) is currently the dominant method for modelling behavioural syndromes in humans and animals. Although FA is useful for inferring the latent causes underlying trait correlations, it does not account for the pairwise behavioural interactions that also contribute to syndrome structure. Given that latent factors and pairwise interactions are likely ubiquitous causes of trait covariation, both should be modelled simultaneously. Currently, however, behavioural ecologists lack an integrative framework for describing and inferring such behavioural syndromes. Generalized network modelling (GNM), representing an integration of FA and Gaussian graphical modelling (GGM), meets this challenge. We provide a theoretical introduction to GNM as well as a method for detecting latent factors in GGMs called exploratory graph analysis (EGA). We then propose the novel EGA+GNM framework for modelling multiple sources of trait correlations and ensuring more robust causal inferences. To empirically demonstrate the utility of this framework, we compare models derived from EGA+GNM and FA using observational measures of social and arousal behaviour in common marmosets Callithrix jacchus. Using information‐theoretic model comparison, we find support for EGA+GNM models compared to models generated by FA. Two EGA+GNM models suggest that while latent factors contribute to the emergence of clustered sociability and arousal behaviours, correlations among these traits may also be partially explained by pairwise interactions. Additionally, these behavioural clusters are hypothesized to be causally linked by a positive pairwise interaction between allogrooming and activity level. These results support our claim that EGA+GNM provides a superior and integrative framework for describing behavioural syndromes. Consequently, by simultaneously modelling both latent factors and pairwise interactions, behavioural ecologists can better understand the evolutionary causes and consequences of animal personality. A formal overview of the EGA+GNM framework and a R tutorial demonstrating its application are provided in the electronic Supporting Information .


| INTRODUC TI ON
Phenotypic integration, which refers to the ultimate and proximate bases of organismal trait covariation, is currently a central topic in evolutionary ecology (Armbruster, Pélabon, Bolstad, & Hansen, 2014;Murren, 2012;Pigluicci & Preston, 2004). Within behavioural ecology, growing attention to phenotypic integration has been apparent in the study of animal personality, which describes consistent amongindividual variation in behaviour Réale, Reader, Sol, McDougall, & Dingemanse, 2007). Research in a variety of taxa has shown that behavioural traits exhibiting personality often correlate across contexts, forming larger clusters termed behavioural syndromes (Sih, Bell, Johnson, & Ziemba, 2004).
Understanding the structure of behavioural syndromes is crucial for quantifying their influence on evolutionary and ecological processes Wolf & Weissing, 2012) and establishing the proximate mechanisms of personality upon which natural selection acts (Araya-Ajoy & Holtmann, Lagisz, & Nakagawa, 2017;Van Oers & Mueller, 2010).
Latent factors are expected to represent the effects of causal processes common to a set of observed traits, such as additive genetic and permanent environmental effects Dochtermann, 2011;Dochtermann, Schwab, & Sih, 2015;Reddon, 2012). The latent factor 'Openness' in bonobos Pan paniscus, for example, which encompasses a syndrome of play behaviour, activity, and neophilia (Martin & Suarez, 2017;Staes et al., 2016), is partially accounted for by variation in the vasopressin receptor gene Avpr1a (Staes et al., 2016). Latent state-behaviour feedback processes can also cause animal personality and behavioural syndromes to emerge . For instance, individual differences in life-history productivity may feedback with and induce correlations among behaviours influencing resource acquisition (Biro & Stamps, 2008).
Although useful for generating models of unobserved common causes for multiple traits, FA is limited in its capacity to capture processes of pairwise interaction between behaviours (Cramer et al., 2012;Goold, Vas, Olsen, & Newberry, 2016;Schmittmann et al., 2013). Pairwise interactions here refer to direct associations between behaviours or states closely proxied by particular behavioural measures (e.g., social dominance, Favati, Leimar, Radesäter, & Løvlie, 2014) that reflect directional or reciprocal causes (e.g., energetic trade-offs or positive feedback processes). For instance, sympatric predation pressure can lead to selection for a behavioural syndrome of aggressiveness, boldness, and exploratory behaviour which can be well described by a latent factor model (Bell & Sih, 2007;Dingemanse et al., 2010). Pairwise correlations between these personality traits can also emerge, however, from feedback processes such as state-dependent safety due to differential body size (Luttbeg & Sih, 2010), immunological capacity (Kortet, Hedrick, & Vainikka, 2010), and contest behaviours and outcomes such as winner-loser effects (Briffa, Sneddon, & Wilson, 2015). The cooccurrence of common causal factors and pairwise interactions within a population may result in multiple pathways to personality and trait correlations, warranting greater consideration of the direct interactions among behaviours and the states they proxy.
Partial correlation network models, also known as Gaussian graphical models (GGMs), have been developed to directly infer pairwise interactions between correlated personality traits (Costantini et al., 2015;Epskamp & Fried, 2016). While GGMs can provide nuanced information about the causes of phenotypic integration (Goold et al., 2016), they do not represent the latent common causes captured by FA. Given that both latent factors and pairwise interactions are likely ubiquitous causes of integrated phenotypes such as behavioural syndromes (Dochtermann, 2011;Murren, 2012;Sih et al., 2015), these processes should be effectively distinguished and modelled simultaneously. Currently, however, behavioural ecologists lack an integrative analytic approach capable of capturing these patterns in their data.
In this paper, we present a novel statistical approach-the EGA+GNM framework-that integrates and overcomes the limitations of current latent factor and network approaches to modelling behavioural syndromes (see Table 1). In particular, EGA+GNM combines generalized network modelling (GNM), a technique for synthesizing FA and GGMs (Epskamp, Rhemtulla, & Borsboom, 2017), with exploratory graph analysis (EGA), which provides a method for detecting latent factors in GGMs (Golino & Epskamp, 2017). To demonstrate the empirical utility of the proposed framework, we compare behavioural syndrome models derived from EGA+GNM and traditional FA techniques using observational measures of social and arousal behaviour in common marmosets Callithrix jacchus.
A formal treatment of our framework (S1) and an extensive R tutorial (S2) are provided in the electronic Supporting Information.
Although we focus on animal personality and behavioural syndromes, the EGA+GNM framework can be applied more broadly to help better understand any integrated phenotype.

| G R APHIC AL MODEL S
Developing causal accounts of trait correlations is crucial for moving beyond superficial characterizations of integrated behavioural phenotypes and uncovering their ecological and evolutionary bases K E Y W O R D S behavioural syndrome, factor analysis, integrated phenotype, marmoset, multivariate, network, personality, social behaviour (Armbruster et al., 2014). In behavioural ecology, it is often impossible to directly demonstrate causal effects, due to the general necessity of observational measurements in naturalistic settings (Niemelä & Dingemanse, 2014;Owens, 2006). Fortunately, however, causal inference can be cautiously pursued in observational contexts through the synthesis of targeted empirical investigation, prior scientific knowledge, and appropriate statistical tools (Grace, 2006;Pearl, 2009;Shipley, 2016;Spirtes & Zhang, 2016). In particular, because causal hypotheses imply specific patterns of conditional independence between traits, observed trait correlations can support or suggest against a causal hypothesis (Grace, 2006;Pearl, 2009).
Conversely, any pattern of observed trait correlations implies some unresolved causal structure (Shipley, 2016), so that causal models can be inferred from observational datasets (Spirtes & Zhang, 2016).
These causal structures can be represented using probabilistic graphical models (Koller & Friedman, 2009;Lauritzen, 1996). A variety of techniques have been developed to estimate and compare the statistical relationships implied by graphical models, as well as to generate plausible graphical models from observed patterns of correlation and conditional independence among traits. Figure 1 provides an overview of the graphical models considered here.

| Factor analysis
Factor analysis (FA) involves the estimation of unobserved variables from their expected effects on measured variables (Bollen, 2002;Skrondal & Rabe-Hesketh, 2004). In the context of animal personality research, FA is used to reduce correlated personality traits to a smaller set of latent factors, which are hypothesized to represent the causal underpinnings of the observed behavioural syndrome.
The causal model underlying FA can be represented using a so-called directed acyclic graph. In Figure 1a, for example, the latent factor U 1 is a common cause of the observed behaviours V 1 , V 2 , V 3 and V 4 .
Importantly, this basic FA model assumes that any observed correlations result solely from the causal effects of the latent factor, as indicated by the absence of direct edges between the observed behaviours. This assumption, referred to as 'local independence' (Bollen, 2002), is the central theoretical motivation linking the statistical estimation and causal interpretation of factor analysis (Haig, 2005).
FA model parameters can be estimated in both an exploratory and confirmatory manner. Exploratory factor analysis (EFA) is particularly useful for hypothesis generation but is limited by its reliance on analytic rotation for estimation. Analytic rotation refers to a process wherein an initially unidentified model is adjusted to minimize a complexity criterion, which subsequently produces unique parameter estimates (Browne, 2001). The various rotation techniques available for EFA are largely guided by ad hoc statistical preferences, such as choices between solutions with more complex trait loadings or factor correlations (Browne, 2001). Applications of EFA therefore often reflect historical practice rather than biologically motivated considerations.
Confirmatory factor analysis (CFA), in contrast, requires the specification of identifiable models based upon past theory, providing a means for more rigorous hypothesis testing and model comparisons both within and between datasets (Dingemanse et al., 2010;Martin & Suarez, 2017). CFA also facilitates the inclusion of latent factors into larger causal hypotheses through the broader framework of structural equation modelling (SEM) (Grace, 2006;Shipley, 2016). Although CFA facilitates behavioural syndrome model comparison and selection, the SEM framework can constrain researchers to ignore unexpected but important phenotypic relationships that TA B L E 1 Analytic techniques for modelling behavioural syndrome structure Statistical method (abbr.)

Description References
Factor analysis (FA) Generate and test continuous latent variable models Araya-Ajoy and Dingemanse (2014), Dingemanse et al. (2010) and Martin and Suarez (2017) Structural equation modelling (SEM) Test causal models containing latent variables Grace (2006) and Shipley (2016) Exploratory structural equation modelling (ESEM) Extend exploratory FA procedures into SEM Marsh et al. (2010Marsh et al. ( , 2014 Gaussian graphical modelling (GGM) Generate and test pairwise partial correlation network models Costantini et al. (2015), Goold et al. (2016) and Epskamp and Fried (2016) Exploratory graph analysis (EGA) Detect latent variables in GGMs using community detection algorithms Golino and Demetriou (2017) and Golino and Epskamp (2017) Generalized network modelling(GNM) Generate and test GGMs for latent factor and residual trait correlations Epskamp et al. (2017) Generalized linear mixed-effects modelling (GLMM) Extend linear regression to incorporate random effects and non-Gaussian responses Bolker et al. (2009), Dingemanse and Dochtermann (2013) and Nakagawa et al. (2017) Notes. The method abbreviations listed here are used throughout the paper. The primary references provide general overviews of these models in the context of biological applications and/or personality research.
may have been uncovered through exploratory investigation. Most SEM software provide modification indices for semi-exploratory CFA model revision, but these post hoc adjustments often exhibit low generalizability (Boomsma, 2000;MacCallum, Roznowski, & Necowitz, 1992). Yet failing to address such unexpected associations can also lead to structures well supported by EFA exhibiting appreciably worse fit in CFA (e.g., Vassend & Skrondal, 1997). Moreover, residual zero-order correlations may indicate spurious associations rather than direct causal relationships, but neither FA nor (E)SEM can identify such confounding effects using partial correlations. Therefore, although CFA and SEM provide a more theoretically robust approach to modelling behavioural syndromes than EFA and ESEM, the assumption of local independence limits the capacity of all FA techniques to adequately describe the complex and multiply determined phenotypes uncovered in animal personality research.

| Generalized network modelling
Generalized network modelling (GNM) overcomes these limitations of FA and (E)SEM . GNM formalizes CFA and SEM within the broader framework of mixed graphical models containing both directed and undirected edges. In particular, GNM represents two distinct generalizations of SEM, whereby the relationships between both latent factor and residual trait values are represented as Gaussian graphical models (GGMs). GGMs are undirected graphs, also known as Markov random fields, that encode conditional independence relationships among multivariate normal random variables (Lauritzen, 1996). A GGM is therefore a network of partial correlations, where each correlation controls for the effect of all other correlations in the network. In Figure 1b, In contrast to the causal effects indicated by directed edges in FA and SEM, the undirected edges in a GGM can be conceptualized as pairwise interactions between observed traits, which may represent both directional and reciprocal causal processes (Costantini et al., 2015;Epskamp & Fried, 2016;Goold et al., 2016). Importantly, the undirected edges encoded by a GGM completely specify a unique set of conditional independence relationships, so that there are no equivalent GGMs possible given a particular dataset and estimation procedure. This is in contrast to directed causal graphs such as in CFA and SEM, which typically facilitate a vast set of causally distinct but statistically equivalent representations for the same dataset (Raykov & Marcoulides, 2001). GGMs thus represent a gateway between correlational and causal modelling , as all possible causal models must be consistent with the underlying

GGM.
Within Epskamp et al.'s (2017) GNM framework, GGMs can be utilized to estimate partial correlation networks describing both latent factor and residual trait associations. In Figure 1d, for example, a residual partial correlation between behaviours V 1 and V 4 remains after accounting for the causal effect of latent factor U 1 on V 2 , as well as any other residual correlations between behaviours V 1 , V 2 , V 3 and V 4 . By integrating SEM and GGMs, GNMs therefore provide researchers with appreciable flexibility for simultaneously investigating multiple levels of phenotypic structure.

| Exploratory graph analysis
Both GNM and FA share the a priori assumption that observed trait correlations are at least partially accounted for by latent common causes. As noted above, although a well-fitting latent factor model provides initial plausibility for this hypothesis, a variety of alternative but often unconsidered causal models will make similar or equivalent predictions of trait correlations. For example, human intelligence is often explained with a single latent causal factor 'g', but alternative models positing developmental feedback between distinct cognitive processes can predict statistically equivalent performance outcomes (Van der Maas et al., 2006). An overemphasis upon unobserved causes may therefore obscure direct causal interactions between behavioural traits (Cramer et al., 2012;Goold et al., 2016;Schmittmann et al., 2013). This general problem of model equivalence also underscores the importance of theoretically informed model testing (Keith, Caemmerer, & Reynolds, 2016;Skrondal & Rabe-Hesketh, 2004) and subsequent empirical investigation to uncover the potential causal mechanisms represented by latent factors (Shipley, 2016). Researchers often lack sufficient information about the plausibility of factor models of behavioural syndrome structure, however, such that the application of GNM  in both exploratory and confirmatory contexts will benefit from additional justification for the presence of latent factors.
Exploratory graph analysis addresses this issue by providing a method for detecting latent variables in GGMs using community detection algorithms (Golino & Epskamp, 2017 can subsequently be applied to these sparse GGMs to uncover plausible latent factors. Importantly, by partitioning edges within and outside latent clusters, EGA can also assist in the identification of residual pairwise interactions. Figure 1c, for instance, represents a hypothetical EGA procedure, resulting in a GGM cluster consistent with a latent factor underlying behaviours V 1 , V 2 and V 3 , as well as an independent partial correlation between behaviours V 2 and V 4 . Additionally, the degree of network modularity determined by community detection algorithms can be used as an estimate of phenotypic modularity, which refers to the degree of semi-autonomous trait clustering within a complex character (Murren, 2012). This EGA procedure is more theoretically motivated than traditional EFA techniques and strongly outperforms them in the accurate recovery of factor models (Golino & Demetriou, 2017;Golino & Epskamp, 2017).

| THE EGA+GNM FRAMEWORK
Taken together, GNM  in coordination with EGA (Golino & Epskamp, 2017) provides an integrative statistical framework that can overcome the limitations of FA for exploring and modelling the causal structure of integrated behavioural phenotypes. Through EGA, researchers can more rigorously assess whether latent factors underlie observed trait correlations and identify pairwise interactions independent of this structure; with GNM, the model(s) suggested by EGA-including the special cases of pure CFA or GGM structures-can be appropriately estimated. Figure 1d, for example, represents the hypothetical EGA procedure in Figure 1c as a GNM.
Our proposed EGA+GNM framework consists of a four-step process for generating and comparing plausible graphical models of behavioural syndrome structure from repeatedly measured behavioural traits: 1. Assess trait repeatability.

Estimate among-individual trait correlations.
3. EGA: convert correlations to a regularized GGM and apply a community detection algorithm.

GNM: estimate and compare models suggested by EGA.
See Supporting Information S1 for a formal overview of our framework. EGA (3) can be directly implemented for crosssectional data with the 'EGA' package for the R statistical environment (Golino & Epskamp, 2017), and GNM model estimation and comparison (4) can be conducted using the 'lvnet' package . Longitudinal data are required, however, to estimate the degree of personality exhibited across traits (1) and effectively distinguish (co)variation resulting from among-and within-individual correlations (2). These steps are crucial for accurate causal modelling, as raw phenotypic correlations can produce biased estimates of the size and direction of among-individual trait correlations (Dingemanse, Dochtermann, & Nakagawa, 2012).
Given that the GGM assumes multivariate normal data (Epskamp & Fried, 2016), multi-response GLMMs can also be used to appropri- for low power datasets, which are common in behavioural ecology (Nakagawa, 2004) and often provide an exploratory basis for future confirmatory analyses in larger samples.

| EMPIRICAL DEMONSTRATION
We now provide an empirical demonstration of the EGA+GNM framework using observational measures of social and arousal behaviours in a laboratory population of common marmosets (see Table 2). This empirical application serves as a comparison of EGA+GNM and traditional FA techniques, as well as an example of the utility a Bayesian EGA+GNM implementation provides in the context of exploratory research using a low power sample. We con- individual variation in these measures. We further hypothesized a priori a behavioural syndrome structure consisting of two latent factors causing sociability (contact sitting, allogrooming, and social proximity) and arousal (activity, scent-marking, and gnawing) behaviours to correlate. Following the procedure outlined above, we began by assessing the degree of repeatability in the measured behaviours.
We then extracted the among-individual trait correlations and subsequently applied both EGA and traditional FA approaches in an exploratory manner to generate multiple graphical models of behavioural syndrome structure. Finally, we used GNM to compare these models and assess which structures provide more plausible representations of the causal pathways connecting the measured behaviours.

| Repeatability assessment
We assessed repeatability using Bayesian GLMMs (McElreath, 2016). Beta GLMMs appropriate for continuous proportions were utilized for the duration measures of social behaviour, while Poisson GLMMs were employed for the count measures of arousal behaviour. These models were estimated using the 'brms' package (Bürkner, 2017) (Martin & Suarez, 2017). By partitioning variance in average behaviour across observational periods, both shortterm (R shortterm ) and long-term (R longterm ) repeatability can also be calculated, which represent the total proportion of model variance within and across observational periods due to individual differences (Araya-Ajoy et al., 2015). After adjusting for fixed effects, R longterm corresponds to the common measure of adjusted repeatability (Nakagawa & Schielzeth, 2010). The delta method  was used to derive the observation-level variance for the Beta and Poisson latent scales, which are necessary to calculate these ratios.
All models included month of observation as a fixed effect to control for potential temporal effects, as well as subject and social group identity as random effects. A so-called series random effect was also included for repeatability estimation, which captured the variance in average individual responses across each observational period (Araya-Ajoy et al., 2015). We used weakly regularizing priors, which place lower expected probability on large values, to achieve more conservative estimates. Regularizing priors are particularly appropriate for analyses conducted with small samples, as they reduce the risk of inferring an effect that does not exist or is in the wrong direction relative to flat or highly diffuse priors (Gelman & Tuerlinckx, 2000).

| Model generation
We fit a multi-response model correlating subject intercepts across behaviours to estimate among-individual trait correlations, and we used the resulting posterior median correlation matrix for exploratory model generation (see Figure 2b). In addition to a null model (M1) and our initial hypothesis (M2), we generated six models from the application of FA approaches and the proposed EGA+GNM framework (see Figure 3 for an overview of the model set). FA methods provided support for two-factor solutions, with an unspecified causal association among the factors (M3), an effect of the latent sociability factor on activity level (M4), two factors with a residual zero-order correlation between allogrooming and activity (M5), and some causal effect of both latent factors across all traits (M6). The application of EGA+GNM (see Figure 2b) suggested two candidate models representing a sparse GGM with missing edges in each behavioural cluster (M7), and two factors with a residual partial correlation between allogrooming and F I G U R E 2 EGA+GNM model generation. (a) Reaction norm, short-term, and long-term repeatability estimates determine the proportion of variance explained by the behavioural syndrome structure, and each is represented here by its posterior median value ± SD. (b) The median correlation network portrays zero-order correlations among the behavioural traits and was utilized for model generation with FA.
The EGA network, which shows the regularized partial correlations among traits conditional on the network structure, was used to generate GGM and GNM models. Edge colours correspond to the direction of the correlation, with blue and grey corresponding to positive and negative correlations respectively. The width and translucency of edges encode the relative magnitude of correlations, which are faded with a nonlinear gradient to enhance visibility. The node colours for the EGA network correspond to clusters detected by the community detection algorithms. See Table 2 for behavioural trait abbreviations activity (M8). Note that this GNM, while superficially similar to the residual correlation factor model (M4), only assumes that the residual GGM is sparse. Zero-order correlations among residual trait values are therefore not directly constrained, allowing greater local dependence without additional model complexity.

| Model comparison and selection
We fit our model set in the 'lvnet' package  and utilized information-theoretic model comparison to assess the relative fit of each model across the posterior of personality trait correlations. The posterior percentage of admissible factor solutions was calculated as a measure of model stability for M2-M6 and M8.
Inadmissible factor model solutions contain zero or negative residual trait estimates, which often occur because of missing model parameters or small residual variances near zero (Kolenikov & Bollen, 2012). Inadmissible posterior solutions were discarded and the resultant posterior median EBIC values (ẼBIC) were compared across models. Reported ΔẼBIC therefore reflect the difference in Ẽ BIC between alternative models and the model with Ẽ BIC min . We considered ΔẼBIC > 2 as providing minimal evidence and ΔẼBIC > 10 as providing strong evidence for reduced model quality relative to the best fitting model (Burnham, Anderson, & Huyvaert, 2011). Uncertainty in the parameters of the best supported models was quantified by estimating the posterior probability of observing a positive effect size (p >0 ). no support compared to the EGA+GNM solutions, supporting the claim that partial correlation networks are generally more informative than zero-order correlations (Costantini et al., 2015;Epskamp & Fried, 2016;Goold et al., 2016). The pure FA models were also highly unstable compared to the GNM (M2-M7 admissible sample range:

| Results and discussion
9%-19%; M8: 91%), which suggests that the factor model alone does not provide an appropriate representation of trait covariation.
Our results collectively suggest that common causal factors contribute to the emergence of clustered sociability and arousal F I G U R E 3 Model set for comparison. For each model, circles represent latent variables (η) and boxes represent observed variables. Single-headed arrows encode the expectation of directional causation, double-headed arrows describe zero-order correlations arising from an unspecified causal process, and undirected lines represent partial correlations conditional on the residual correlations among all observed traits. See Table 2 for behavioural trait abbreviations behaviours, but also that correlations among these traits may be partially explained by pairwise interactions. This is reflected in the structure of the sparse GGM model ( see Figure 4). The uncertainty in these posterior estimates reflects the low power of our sample, but the results overall provide consistent evidence for positive effects. The sociability and arousal clusters discovered here are well supported by and consistent with past marmoset personality research. In particular, the sociability cluster is consistent with the behavioural functions of the "Sociability" (Inoue-Murayama et al., 2018) and "Agreeableness" syndromes (Iwanicki & Lehmann, 2015;Koski et al., 2017) previously described with rating methods, while the arousal cluster is consistent with "Inquisitiveness" (Koski et al., 2017), "Openness" (Iwanicki & Lehmann, 2015), and the locomotor component of "Stress-Activity" (Šlipogor et al., 2016).
Gnawing often co-occurs with scent-marking, which likely enhances the adhesion of the deposited chemical cues (Massen, Šlipogor, & Gallup, 2016). These behaviours are expressed more frequently at points of direct (Lazaro-Perea, Snowdon, & de Fátima Arruda, 1999) and indirect  olfactory contact between groups and may function both for establishing territorial boundaries and eliciting among-group mating opportunities (Heymann, 2006;Lazaro-Perea et al., 1999;Lledo-Ferrer, Peláez, & Heymann, 2011). The potential mediational effect of gnawing within the arousal cluster may therefore reflect the causal influence of activity level on gnawing behaviour, which tends to occur while an individual moves along the perimeter of their group territory and subsequently facilitates scent-marking. In those cases where F I G U R E 4 Parameters for the best supported behavioural syndrome models. (a) Posterior factor loading and residual partial correlation network coefficients for the best supported GNM model (M8). (b) Posterior partial correlation coefficients for the sparse GGM model (M7). The posterior probability of observing a positive effect size (p >0 ) is listed for each model parameter. Additionally, each posterior is shaded by its empirically estimated tail probability distribution, also known as the complimentary cumulative density function, which encodes the probability of observing an outcome at least as extreme as the current value. Darker colours therefore represent regions closer to the expected posterior value. Samples with <1% of the maximum density have been suppressed for visual clarity. See Table 2 for behavioural trait abbreviations scent-marking reflects territorial defense, all three behaviours may be caused by latent arousal factors induced through direct or indirect cues of conspecific presence.
The pairwise interaction found between allogrooming and activity level is consistent with previously reported links between facets of activity and sociability in hominids, including both behavioural and rating measures of human (Wilson & Dishman, 2015), chimpanzee Pan troglodytes (Pederson, King, & Landau, 2005), bonobo (Garai, Weiss, Arnaud, & Furuichi, 2016), gorilla Gorilla beringei (Eckardt et al., 2015), and orangutan Pongo pygmaeus and P. abelii (Weiss, King, & Perkins, 2006) personality. As a preliminary hypothesis for further research, this association between allogrooming and activity, independent of contact sitting and social proximity, may reflect personality in social sampling behaviour, as more gregarious and proactive individuals interact with numerous social partners and maintain a larger proportion of weak social network ties (Aplin et al., 2013;Sih & Del Giudice, 2012).
Individuals may therefore exhibit consistent differences in both their aggregate sociability and how they distribute their social behaviour across available partners (Aplin et al., 2015), resulting in differential associations among the observed social behaviours and their causal interaction with activity.

| CON CLUS ION
We presented a theoretical introduction to exploratory graph analysis (Golino & Epskamp, 2017) and generalized network modelling , and we argued that an integration of these approaches-the EGA+GNM framework-will enhance descriptive modelling and causal inference in behavioural syndrome research.
As demonstrated in our exploratory empirical example, GGM and GNM techniques can provide additional nuance and causal insight beyond the factor analytic approaches prominent in personality research (Araya-Ajoy & Dingemanse et al., 2010;Martin & Suarez, 2017;Weiss, 2017). The EGA+GNM framework facilitates a theoretically motivated model generation procedure integrating exploratory and confirmatory approaches to model comparison and selection. By employing EGA+GNM within a Bayesian framework, the uncertainty in this process can be effectively represented and carried across stages of analysis. The models best supported by our data would not have been uncovered through the application of traditional factor analytic techniques, and we therefore encourage other researchers to apply our EGA+GNM framework in subsequent animal personality research, as well as in research on other integrated phenotypes more broadly.

ACK N OWLED G EM ENTS
We are very grateful to Tanja

DATA ACCE SS I B I LIT Y
The data and R code used for our empirical demonstration are available through the Dryad Digital Repository https://doi.org/10.5061/ dryad.5f88r9m (Martin et al., 2018).