Introduction
Internalizing symptoms are a particular concern in the early teenage years, with increased prevalence of anxiety and depression, and frequent first lifetime onset in this age range (NHS Digital, 2018; Rapee et al., Reference Rapee, Oar, Johnco, Forbes, Fardouly, Magson and Richardson2019; Solmi et al., Reference Solmi, Radua, Olivola, Croce, Soardo, Salazar de Pablo and Fusar-Poli2022). Understanding the role of risk factors to aid prevention and early identification of adolescent internalizing symptoms is therefore a major research priority. To date, such work has focused on simple formulations of mental health difficulty, drawing on diagnoses, latent constructs, and total scores (for reviews see e.g. Brumariu & Kerns, Reference Brumariu and Kerns2010; Evans, Li, & Whipple, Reference Evans, Li and Whipple2013; Rapee et al., Reference Rapee, Oar, Johnco, Forbes, Fardouly, Magson and Richardson2019). While these approaches have provided insight, they are not in line with current theory in which the interacting complex nature of symptoms is conceptualized as a network (Borsboom, Reference Borsboom2017). We therefore set out to analyze the effect of previously identified risk factors on properties of internalizing symptom networks in mid adolescence to gain insight into whether global properties vary according to risk exposure.
The network theory of mental disorders holds that rather than a single process causing a given set of symptoms, a state of disorder arises from pathological patterns in how symptoms interact (Borsboom, Reference Borsboom2017; Robinaugh, Hoekstra, Toner, & Borsboom, Reference Robinaugh, Hoekstra, Toner and Borsboom2020). Two key hypotheses have been offered to explain differences between normal and disordered networks, namely altered connectivity and structure. Global connectivity refers to the sum total of all the relationships in the network (i.e. between all included symptoms), with disordered networks hypothesized to show stronger overall relationships between symptoms (Borsboom et al., Reference Borsboom, Rhemtulla, Cramer, van der Maas, Scheffer and Dolan2016; Robinaugh et al., Reference Robinaugh, Hoekstra, Toner and Borsboom2020). The altered structure hypothesis predicts that pathological and healthy networks may vary in terms of how they are connected, with certain symptoms interacting with one another in one group but not, or to a lesser extent, in the other (Borsboom et al., Reference Borsboom, Rhemtulla, Cramer, van der Maas, Scheffer and Dolan2016). These hypotheses have been applied in similar work considering stressors and substance abuse (Lin, Fried, & Eaton, Reference Lin, Fried and Eaton2020b), the effects of screen time (Lin, Eaton, & Schleider, Reference Lin, Eaton and Schleider2020a), and genetic/environmental risk (van Loo et al., Reference van Loo, Van Borkulo, Peterson, Fried, Aggen, Borsboom and Kendler2018).
While networks might often ideally be considered using intraindividual (i.e. longitudinal) data since this provides direct insight into mental health processes, where individual differences such as risk factors are of interest, cross-sectional networks are considered appropriate (Borsboom et al., Reference Borsboom, Deserno, Rhemtulla, Epskamp, Fried, McNally and Waldorp2021). We aimed to consider altered connectivity and structure after conditioning on a given risk factor using cross-sectional data as others have done (Lin et al., Reference Lin, Eaton and Schleider2020a, Reference Lin, Fried and Eaton2020b; McElroy, Fearon, Belsky, Fonagy, & Patalay, Reference McElroy, Fearon, Belsky, Fonagy and Patalay2018). Interindividual networks used in this way could provide insight into how symptoms are associated with one another at the between-person level given a known a risk factor. For instance, those who are bullied may be more likely to feel sad if they are also worried than those not in the bullied risk group. This approach is also consistent with risk-factor literature which typically focuses on measurable variables that are likely to be associated with increased risk on average.
We therefore sought to capitalize on network theory and methods to robustly assess the relationship between global network properties and exposure to a range of established risk factors. Risk factors are measurable attributes, characteristics, and exposures that predict negative outcomes in a given domain of functioning (Furber, Leach, Guy, & Segal, Reference Furber, Leach, Guy and Segal2017). Very little work has considered whether adolescent internalizing networks are globally sensitive to known risk factors, with only a few studies considering age and sex, with most considering only structural differences (Abend et al., Reference Abend, Bajaj, Coppersmith, Kircanski, Haller, Cardinale and Pine2021; Black, Panayiotou, & Humphrey, Reference Black, Panayiotou and Humphrey2022; McElroy et al., Reference McElroy, Fearon, Belsky, Fonagy and Patalay2018; Monk, McLeod, Mulder, Spittlehouse, & Boden, Reference Monk, McLeod, Mulder, Spittlehouse and Boden2023). In addition, studies of other populations and known risks have yielded equivocal findings, with inconsistent and sometimes null results for both connectivity and structure (Lin et al., Reference Lin, Fried and Eaton2020b; van Loo et al., Reference van Loo, Van Borkulo, Peterson, Fried, Aggen, Borsboom and Kendler2018). There is therefore a clear need to clarify the relationship between risk factor exposure and global network properties.
Building on this, we also consider cumulative risk theory, the core tenet of which is that negative outcomes are better predicted by compound risk exposure than single risk factors in isolation. In other words, it is the number rather than the nature of risk factors that best predicts vulnerability to psychopathology. The underpinning assumption is that exposure to increasing numbers of risk factors produces chronic stress (i.e. increased allostatic load), ultimately leading to dysfunction (Evans et al., Reference Evans, Li and Whipple2013). The cumulative risk model has garnered considerable empirical support (see Evans et al., Reference Evans, Li and Whipple2013 for a review of studies focusing on children and adolescents). Furthermore, a recent analysis of number of life stressors and substance use networks also showed some cumulative effects (Lin et al., Reference Lin, Fried and Eaton2020b). Notably, this and other work has provided evidence of curvilinear trends (quadratic or cubic models) for cumulative risk exposure. For example, Ashworth and Humphrey (Reference Ashworth and Humphrey2020) found evidence of threshold effects in which marked increases in behavior and reading problems were seen after exposure to specific numbers of risk factors. However, as found in the work on substance abuse networks (Lin et al., Reference Lin, Fried and Eaton2020b), models can also suggest plateau effects at certain numbers of risk factors.
Given a lack of studies considering network properties in relation to risk exposure, we sought well-evidenced risk factors (from other analytical approaches), drawing on reviews/meta-analyses where possible, to aid interpretation of results (though we also note a key pragmatic consideration – availability in the #BeeWell dataset used in the current study). Based on this, the following risk factors, which have been found to be significantly associated with internalizing symptoms among adolescents, were selected: age (Deighton et al., Reference Deighton, Lereya, Casey, Patalay, Humphrey and Wolpert2019; Deighton, Yoon, & Garland, Reference Deighton, Yoon and Garland2020; NHS Digital, 2018), gender identity and sexuality (Connolly, Zervos, Barone, Johnson, & Joseph, Reference Connolly, Zervos, Barone, Johnson and Joseph2016; Black, Humphrey, & Marquez, Reference Black, Humphrey and Marquez2023; Plöderl & Tremblay, Reference Plöderl and Tremblay2015), low parent/carer (Gariépy, Honkaniemi, & Quesnel-Vallée, Reference Gariépy, Honkaniemi and Quesnel-Vallée2016) and peer (Rueger, Malecki, Pyun, Aycock, & Coyle, Reference Rueger, Malecki, Pyun, Aycock and Coyle2016) support, income deprivation (Deighton et al., Reference Deighton, Lereya, Casey, Patalay, Humphrey and Wolpert2019; Reiss, Reference Reiss2013), having special educational needs (SEN; Deighton et al., Reference Deighton, Lereya, Casey, Patalay, Humphrey and Wolpert2019; Patalay & Fitzsimons, Reference Patalay and Fitzsimons2016), bullying victimization (Moore et al., Reference Moore, Norman, Suetani, Thomas, Sly and Scott2017), and poor physical health (Patalay & Fitzsimons, Reference Patalay and Fitzsimons2016).
The current study
Based on the above, we set out to consider the effects of risk factors individually and cumulatively on global properties of internalizing symptom networks at the between-person level. We hypothesized that the presence of a given individual risk factor would be associated with greater network connectivity and altered structure. Similarly, we hypothesized that higher cumulative risk exposure would result in higher connectivity and altered network structure. We also explored whether the relationship between cumulative risk exposure and network connectivity was linear or curvilinear. To test these hypotheses, we undertook secondary analysis of the first annual wave (2021) of the #BeeWell study (#BeeWell Research Team, 2021), which comprised N = 37 978 adolescents aged 12–15 across 10 Local Authorities in the Greater Manchester city-region in England.
Method
Participants
After removing participants with missing on all variables in the current analysis, since we drew from a wider survey, and excluding 269 from other year groups, the final sample size was N = 34 564 adolescents aged 12–15, drawn from 159 schools. The sample were comparable to national figures for sex (50.12% male compared to the national figure of 50.26%); distribution across year groups (53.27% in Year 8 compared to 50.89% nationally); English as an additional language (20.90% compared to 19.50% nationally); and special educational needs (13.18% compared to 14.10% nationally; HM Goverment, 2022). In terms of deprivation, 25.57% were from the top income deprivation decile considering neighborhoods (most deprived), and only 3.77% were in the bottom decile compared to 10% per decile nationally (Ministry of Housing Communities and Local Government, 2019). Ethnicity was also somewhat different in the current sample compared to national figures. For instance, 65.02% were white compared to the national figure of 70.80% (HM Goverment, 2022). Further details of demographic characteristics can be seen in Table 1 with national figures drawn from official summaries given for illustrative purposes.
Table 1. Sample characteristics
Note: National data derived from HM Goverment (2022).
a National data are not available for FSM eligibility in the last six years. 20.9% of pupils aged 11–16 are currently eligible nationally.
b Deciles are calculated by ranking the 32 844 neighborhoods in England from most deprived to least deprived and dividing them into 10 equal groups, ranging from the most deprived 10% of neighborhoods to the least deprived 10% of neighborhoods nationally (Ministry of Housing Communities and Local Government, 2019).
Measures
Internalizing symptoms
Though the overall internalizing subscale of Me and My Feelings has 10 items (Deighton et al., Reference Deighton, Tymms, Vostanis, Belsky, Fonagy, Brown and Wolpert2013), we opted to use only the eight that captured conceptually distinct symptoms to avoid noise in connectivity due to similar content. Inclusion of conceptually similar items is at odds with network theory and can bias parameters (Fried & Cramer, Reference Fried and Cramer2017; Hallquist, Wright, & Molenaar, Reference Hallquist, Wright and Molenaar2021). The measure contains two items each for sleep and worry which have been shown elsewhere to have local dependence (Black, Panayiotou, & Humphrey, Reference Black, Panayiotou and Humphrey2019), further suggesting the inappropriateness of including both of each pair in a network analysis (Christensen, Garrido, & Golino, Reference Christensen, Garrido and Golino2023). In both cases we opted for the more general items: ‘I worry a lot’ rather than the school worry item, and ‘I have problems sleeping’ rather than the night waking item. While some work suggests averaging similar items (de Ron et al., Reference de Ron, Robinaugh, Fried, Pedrelli, Jain, Mischoulon and Epskamp2022), we had only two items per symptom, had validity reasons to prefer one item over the other, and were also mindful that combining items favors sensitivity (while we sought specificity since we focused on global network properties).
Risk factors
Risk factors were coded as present/absent (i.e. binary). For some variables this was consistent with the way there were collected, while for others we had to construct a cut off (see online Supplementary Table S1). This approach is consistent with cumulative risk theory and modeling (Evans et al., Reference Evans, Li and Whipple2013), and allowed us to reduce the number of structural comparisons since methods for handling many groups are limited for networks. Since we aimed to capture those most at risk, we aimed to create cut-offs likely to represent high levels of risk exposure while preserving sample size. For instance, for income deprivation we opted to combine a marker of low family income (FSM eligibility) with neighborhood deprivation, since the child-level FSM metric held as part of #BeeWell is considered somewhat inclusive (Boliver, Gorard, & Siddiqui, Reference Boliver, Gorard and Siddiqui2022).
Cumulative risk
We operationalized cumulative risk by summing the risk variables above (0–9). Inspection of this revealed low cell size for groups at five risk factors and higher. We therefore collapsed these into a 5+ group (N = 604), consistent with established practice in cumulative risk research (e.g. Ashworth & Humphrey, Reference Ashworth and Humphrey2020).
Analysis
R code with fake data and bootstrapped matrices are available to run the analysis and reproduce results where possible (see supplementary material).
Network estimation
Networks consist of nodes (individual symptoms or items), which are represented in diagrams by circles, and edges, the relationship between a given pair of nodes, represented by lines, with thickness reflecting the relative magnitude of an effect. Consistent estimation procedures were used across networks, though sample size varied depending on the risk factor. Since dense networks, as can be expected for adolescent internalizing symptoms (Black et al., Reference Black, Panayiotou and Humphrey2022; McElroy et al., Reference McElroy, Fearon, Belsky, Fonagy and Patalay2018), are best estimated using non-regularized methods (Burger et al., Reference Burger, Isvoranu, Lunansky, Haslbeck, Epskamp, Hoekstra and Blanken2023), and there was some missing data, we opted to use full information maximum likelihood estimation via the psychonetrics package (Epskamp, Reference Epskamp2021). Combined with pruning (removing non-significant edges) at α = 0.01, this has been shown to have specificity of 0.95 and sensitivity of 0.80 for skewed ordinal data with a relatively small number of nodes and N = 5000 (Isvoranu & Epskamp, Reference Isvoranu and Epskamp2023).
Comparison of covariance structure
Several methods are available to consider network invariance across groups (Burger et al., Reference Burger, Isvoranu, Lunansky, Haslbeck, Epskamp, Hoekstra and Blanken2023). All are relatively new, underpinned by very few simulation studies (Kan, van der Maas, & Levine, Reference Kan, van der Maas and Levine2019; van Borkulo et al., Reference van Borkulo, van Bork, Boschloo, Kossakowski, Tio, Schoevers and Waldorp2022), reflecting the novelty of the psychological network field in general. Given this, and to address our hypotheses of focusing on global properties, we drew on several approaches. While such novel analyses can represent increased degrees of freedom, we provide reproducible code and synthetic data, in line with best practice (Epskamp, Reference Epskamp2019). First, we considered whether there was global evidence of covariance structure differences via recursive partitioning, using the model-based algorithm (Jones, Mair, Simon, & Zeileis, Reference Jones, Mair, Simon and Zeileis2020). If a significant difference was found here, we estimated networks in each group separately to begin structural comparisons (see below).
Bootstrapping procedure
To allow subsequent structural and connectivity comparisons, data were bootstrapped with 1000 repetitions at the smallest (cumulative) risk-group sample size, following the procedure set out by Lin et al. (Reference Lin, Fried and Eaton2020b). This enabled us to balance sample sizes to improve power (Yoon & Lai, Reference Yoon and Lai2018), avoid bias in standard errors (Delacre, Lakens, & Leys, Reference Delacre, Lakens and Leys2017), provide some insight into the distribution (i.e. reliability) of results (Epskamp, Borsboom, & Fried, Reference Epskamp, Borsboom and Fried2018), and model the relationship between connectivity and group membership (via the bootstrapped distribution).
Structural comparison
The adjacency matrices (made up of 0s and 1s to represent the pattern of estimated edges) for each group were compared, based on edges that were included in > 50% of bootstrap samples. This enabled us to check if the same pattern of edges could be compared between groups. This ≥ 50% threshold followed other similar work (Lin et al., Reference Lin, Fried and Eaton2020b), and avoids known problems with thresholding on bootstrapped confidence intervals where model selection has already taken place (Burger et al., Reference Burger, Isvoranu, Lunansky, Haslbeck, Epskamp, Hoekstra and Blanken2023). If adjacency matrices were not identical, networks were considered to be structurally different, and we visualized qualitative differences. If adjacency matrices were equal, confirmatory multi-group network testing could be performed to see if imposed equality constraints across groups could be maintained based on comparison of the Akaike information criterion and Bayesian information criterion (Black et al., Reference Black, Panayiotou and Humphrey2022; Kan et al., Reference Kan, van der Maas and Levine2019).
Connectivity comparison
Based on the bootstrapped networks, we calculated global strength as a measure of connectivity, the absolute sum of all edge weights, and modeled the effect of group membership (risk present v. risk absent) on this outcome, and treated cumulative risk as continuous consistent with theory and previous work (Evans et al., Reference Evans, Li and Whipple2013; Lin et al., Reference Lin, Fried and Eaton2020b). We report unstandardized effects for individual risk factors given consistent metrics were used and challenges inherent in standardizing binary variables, but additionally include standardized effects for the continuous cumulative risk model. For cumulative risk we then also fitted quadratic and cubic models to assess curvilinear trends and considered such higher-order polynomial models to be preferred only when the change in R 2 was > 0.01 given that such models are likely to show better fit (Lin et al., Reference Lin, Fried and Eaton2020b). In addition, we did not interpret absolute R 2 for a given model since this is affected by the bootstrapping procedure. We report all values for transparency and to consider change as described.
Results
Structural differences for individual risk factors
Recursive partitioning to determine differences in covariance structure for each of the nine risk factors individually showed significant differences, with Bonferroni-corrected p < 0.05 for every variable except income deprivation (see also online Supplementary Table S2).
Structural differences were therefore considered for every risk factor except income deprivation. As expected for data-driven network analysis, fit was near perfect when networks were estimated separately in each group for the eight risk factors considered, and networks were typically densely connected with low degrees of freedom (see online Supplementary Table S3).
Comparison of adjacency matrices, based on edges present ≥ 50% of bootstrapped networks, suggested the pattern of edges was not the same across risk groups for any risk factor, precluding further statistical comparison of the magnitude of edges. Figure 1 shows average structures across risk factors and groups based on bootstrapping at the sample size of the smaller risk group. The discrepancy in adjacency matrices can be seen in Fig. 1 with for instance the shy-worry and sleep-scared edges clearly missing in the risk group for parent/carer relationships compared to the no-risk group. As indicated in the annotations in Fig. 1 (based on number of unique estimated edges in the averaged adjacency matrices divided by the number of possible edges, (8 × 7)/2 = 28), density was slightly lower in the risk group for six of the eight risk factors and for the remaining two it was identical across risk present/absent groups.
Figure 1. Network structures for risk factors by group.
Note: The layout is fixed as the average across all networks plotted to aid visual comparison. The maximum edge width (the relative size of the partial correlation) is standardized within risk factors. Density refers to bootstrapped density (edges present in ≥ 50%). Min and Max refer to the minimum and maximum edge weights (partial correlations) for a given network. SEN, special educational needs.
Connectivity for individual risk factors
Connectivity was estimated for every individual risk factor, with the bootstrap procedure applied to all variables, including income deprivation. Results of the individual regression models with group membership predicting connectivity can be seen in Table 2. These indicate that significantly higher connectivity was observed for the risk group for all risk factors except income deprivation and SEN, for which negative coefficients were observed. Effects appeared small except for gender/sexual minorities and parent/carer support.
Table 2. The effect of risk factors on connectivity: results for individual connectivity regression models
Note: CIs, confidence intervals; SEN, special educational needs; **p < 0.01.
Structural differences across levels of cumulative risk
Recursive partitioning across the six levels of cumulative risk (0–5+), suggested significant differences in the covariance structure of all six levels (see online Supplementary Fig. S2). Networks were therefore estimated at the smallest sample size (5+ risk factors, N = 630) with 1000 bootstrap replications for each. Fit of the whole sample (not bootstrapped) model in each of the five levels of risk selected based on recursive portioning is in online Supplementary Table S3.
Comparison of adjacency matrices based on edges present in ≥ 50% bootstrapped networks revealed no pair of cumulative risk levels to share the same adjacency matrix, precluding further comparison of edge magnitude. Networks across different levels of cumulative risk are visualized in Fig. 2.
Figure 2. Network structures at different levels of cumulative risk.
Note: The maximum edge width (the relative size of the partial correlation) and layout are standardized across levels of cumulative risk. Density refers to bootstrapped density (edges present in ≥ 50%). Min and max refer to the minimum and maximum edge weights (partial correlations) for a given network.
Connectivity across levels of cumulative risk
Results of the regression model comparing levels of connectivity for all levels of cumulative risk suggested connectivity was strongly related to cumulative risk level (see Table 3). However, fitting higher-order models did not result in R 2 change > 0.01, with quadratic R 2 = 0.794**, 95% CI [0.79–0.80] and cubic R 2 = 0.800**, 95% CI [0.79–0.81], meaning we did not find support for a curvilinear effect of cumulative risk on network connectivity.
Table 3. The effect of cumulative risk on connectivity: results for cumulative risk regression model
Note: **p < 0.01; CIs, confidence intervals.
Discussion
The current study drew on a large sample to evaluate whether two related network hypotheses (altered structure and increased connectivity for risk present groups; Borsboom et al., Reference Borsboom, Rhemtulla, Cramer, van der Maas, Scheffer and Dolan2016) were supported considering adolescent internalizing symptoms. Structural differences were observed between risk-present/risk-absent groups and across levels of cumulative risk in all cases except for income deprivation. However, network density and visualization suggested differences in structure were likely small. Similarly, significant effects were seen for all individual risk factors on connectivity though these mostly appeared small, and SEN and income risk groups had slightly lower connectivity than their risk absent counterparts, counter to the connectivity hypothesis. Noteworthy effects on connectivity were seen for gender and sexual minority risk groups as well as those with perceived low parent/carer support. Finally, and consistent with cumulative risk theory, a strong linear effect of cumulative risk on connectivity was observed.
Altered structure for individual risk factors
Drawing on bootstrapping procedures, we found small differences in the average pattern of interactions between symptoms across risk-present/risk-absent groups (except for income deprivation). However, we were unable to quantify these differences in detail since the very finding of a different pattern precluded multi-group comparison as would be conducted in a less data-driven framework such as factor analysis (Brown, Reference Brown2015). Consistent with our analytical decisions described above, we therefore highlight a need to further methods for structural network comparisons, given the pertinence to questions such as the current study, and a clear interest in this area (e.g. Abend et al., Reference Abend, Bajaj, Coppersmith, Kircanski, Haller, Cardinale and Pine2021; Monk et al., Reference Monk, McLeod, Mulder, Spittlehouse and Boden2023). While alternative methods to those used here to determine significant differences are theoretically available (van Borkulo et al., Reference van Borkulo, van Bork, Boschloo, Kossakowski, Tio, Schoevers and Waldorp2022), these are practically limited to regularized estimation which is not appropriate in many cases, including the current study. We highlight this issue since though multi-group testing can be used with non-regularized estimation (Kan et al., Reference Kan, van der Maas and Levine2019), this relies on a common pattern. However, it is likely, given the conditional nature of networks, that small differences could shift this across groups. Nevertheless, in the current study, our additional deployment of recursive partitioning, which has been shown to perform well in large samples estimating small networks as we did (Jones et al., Reference Jones, Mair, Simon and Zeileis2020), supported differences between groups and levels of cumulative risk.
A further consideration that arises from the current study is how altered structure should be interpreted at an interindividual level. The structure hypothesis has typically been introduced contrasting individuals, suggesting that different structures might explain (possibly known) different states (e.g. Borsboom et al., Reference Borsboom, Rhemtulla, Cramer, van der Maas, Scheffer and Dolan2016). However, at the group level (i.e. to usefully consider risk factors) stricter hypotheses may need to be developed. For instance, are there particular symptoms or combinations of symptoms which are considered to be more risky (theoretically) and could these become the focus of analyses? This would limit a general vulnerability of covariance analyses in which many results can flexibly be interpreted to support loose theories (Fried, Reference Fried2020), but requires much more work considering networks and individual symptoms, possibly also drawing on intraindividual studies. We did not make predictions about this, consistent with similar work which has also been exploratory (Abend et al., Reference Abend, Bajaj, Coppersmith, Kircanski, Haller, Cardinale and Pine2021; Monk et al., Reference Monk, McLeod, Mulder, Spittlehouse and Boden2023), given the novelty of the study and limitations in structural comparison methods noted above. However, future work might consider, for instance, whether stronger effects are found between key symptoms such as worry and unhappiness (Black et al., Reference Black, Panayiotou and Humphrey2022), though this would also require guidance on how to compare networks with different edges, as were found here. Such work could build on the proof of concept provided by the current study.
Individual risk factors and connectivity
Since we report unstandardized effects, we focus on the relative size of coefficients across risk variables (given consistent metrics were used). Age, bullying, physical health, peer support, SEN, and income deprivation, had absolute effects b < 0.10, compared to a typical reference group mean ~3.20. One interpretation is that these effects are trivial, and the negative effects for SEN and income deprivation are null, just as for these other risks. This would be consistent with literature suggesting neighborhood deprivation is more consistently associated with externalizing than internalizing symptoms (Visser et al., Reference Visser, Bolt, Finkenauer, Jonker, Weinberg and Stevens2021), and that SEN may show surprisingly limited predictive power for internalizing symptoms, varying by gender (Deighton et al., Reference Deighton, Humphrey, Belsky, Boehnke, Vostanis and Patalay2018). These findings suggest that the way levels of internalizing symptoms predict one another, whether the risk is present or not for these variables, is fairly consistent. However, much more work to understand these effects is needed. Much stronger effects were seen for gender (b = 0.22) and sexual (b = 0.12) minorities, as well low parent/carer support (b = 0.26). These effects are noteworthy since while the latter speaks to the consensus that home relationships are key to functioning (Brumariu & Kerns, Reference Brumariu and Kerns2010; Gariépy et al., Reference Gariépy, Honkaniemi and Quesnel-Vallée2016), the former two risk factors are relatively understudied. Nevertheless, there is emerging evidence these may be particularly at-risk groups due to issues such as stigma, bullying and the stress these create (Connolly et al., Reference Connolly, Zervos, Barone, Johnson and Joseph2016; Plöderl & Tremblay, Reference Plöderl and Tremblay2015), consistent with findings here.
Effects of cumulative risk
The cumulative risk framework enabled us to extend our analysis to include compound risk exposure, which given the interindividual network/bootstrapping procedure could not otherwise be accommodated through a multivariable approach. In terms of structure, much sparser networks were estimated for all levels of cumulative risk compared to across risk-present/risk-absent groups. This likely relates to the smaller sample size used for bootstrapping in the cumulative risk analysis, given that sensitivity is reduced for our selected estimation procedure with lower sample size (Isvoranu & Epskamp, Reference Isvoranu and Epskamp2023). Density should therefore not be compared across individual risk group networks and those for levels of cumulative risk. Nevertheless, recursive partitioning and density suggested differences between every level of cumulative risk, and no pair of levels contained the same pattern of edges. Again, work is needed to formulate hypotheses about how structure might be expected to vary with levels of cumulative risk. Nevertheless, the current study demonstrates a clear justification for further work in this area.
Consistent with this, a strong linear trend for risk exposure predicting connectivity was found. The fact we found support for a linear, as opposed to a curvilinear, trend, suggests cumulative risk exposure can be interpreted as acting additively on global network properties (Evans et al., Reference Evans, Li and Whipple2013). Our study also suggests support for cumulative risk theory, in that it appeared robust to network as opposed to sum/latent score approaches. While previous network analysis had considered multiple stressors (Lin et al., Reference Lin, Fried and Eaton2020b), this study extended the existing method beyond a single life stress inventory, consistent with cumulative risk research. This paves the way for further work analyzing interindividual risks in network analysis, to consider insight into issues such as number v. nature of risks, and draw on additional variables without making large and perhaps poorly defined networks by including risks in the network (Neal & Neal, Reference Neal and Neal2023). The current method also has the advantage of providing insight without conditioning on symptom severity which can create spurious results: Conditioning on a risk factor, which is correlated with elevated symptoms but not caused by these does not induce the bias associated with conditioning on (e.g.) disorder status (de Ron, Fried, & Epskamp, Reference de Ron, Fried and Epskamp2021).
The individual and cumulative risk connectivity models do not directly answer the question of whether the number or nature of risks is most important. This is because the individual models cannot control for anything else (e.g. compete with other risks), making the comparison unfair. Nevertheless, the strong linear effect of cumulative risk on connectivity suggests the number of risks, regardless of their nature, plays a substantial role in how symptoms interact on average. This suggests more risk work should consider effects via network methods, rather than remain rooted in latent/sum-score frameworks. In sum, the current study, as the first (to our knowledge) to apply cumulative risk in this way leveraging secondary data, provides an important framework for future work.
Limitations
This study drew on a large sample size and robust methods to provide novel insights as discussed above. Nevertheless, we acknowledge several limitations. First, we dichomotized multi-categorical risk variables, without validated cut-offs, resulting in information loss. However, we did this to facilitate structural comparisons to address our hypotheses, carefully considered approaches to capture those most at risk, and consistent with cumulative risk theory. Similarly, we were limited by variables in the dataset (e.g. it only covered two age groups), and much more research is therefore needed considering different measures for internalizing symptoms and risks to assess the generalizability of our findings. Second, we did not consider balanced sample sizes in recursive partitioning, which likely affects power (Jones et al., Reference Jones, Mair, Simon and Zeileis2020). Nevertheless, groups were based on a balance of sample size and theory, also supported by our very large total sample size. We also followed recursive partitioning analysis with further structural comparison based on balanced sample sizes, mitigating this risk to some extent. Third, the summing of risks to create the cumulative metric did not take into account missingness, meaning our cumulative risk variable could underestimate cumulative risk. However, missingness was generally low. Fourth, we applied the same estimation procedure across analyses which likely reduced sensitivity for the cumulative risk networks. However, this facilitated a consistent approach and maintained specificity. Finally, though our sample was large and often close to national averages, peculiarities of the Greater Manchester context (e.g. the distribution of ethnicity and income) mean more work considering the generalizability of results is needed.
Conclusion
The current study demonstrates the utility of, and a robust approach to considering the effect of risks on global network properties. It contributes to network and risk theories. A strong effect was seen for cumulative risk on network connectivity, consistent with the ideas that pathological states are associated with higher connectivity, and that the number of risks, regardless of their nature, is important. In terms of individual risks factors, gender/sexual minority status and low parent/carer support had the biggest impact on connectivity, suggesting these are particularly important for identification and prevention. More work is needed to support the interpretation of interindividual structural differences. Nevertheless, altered structures were seen in most cases suggesting this could also play a role in explaining internalizing symptom states.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0033291723002362.
Acknowledgements
We thank the young people who contributed to designing and providing data for #BeeWell, as well as the team at the Greater Manchester Combined Authority and our partner organizations.
Funding statement
#BeeWell is funded by several organizations, see https://gmbeewell.org/partners/funders/.
Ethical standards
The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008. The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional guides on the care and use of laboratory animals.
Open access
