Participants

The total sample used in the analyses included 1,449 women who had participated in two waves of a Finnish population-based study: the Genetics of Sexuality and Aggression study conducted in 2006 (T1; mean age 25.5 years, SD = 4.9) and 2013 (T2). The data collection procedure is described at length elsewhere7,31. To correct for familial dependency, due to the participants being twins and sisters of twins and thus genetically related, one person per family was randomly selected from the original 2,173 women who had submitted data at both time points, resulting in a sample of 1,729 women. Since network estimation needs full information, missing values were imputed for quantitative variables (described further in the Statistical Analyses section). Thirty-seven participants that had missing data on the categorical variable relationship status at T1 were excluded from the analyses. Furthermore, when comparing networks between groups, it is important that the sample size for each group is the same32. As the smallest subgroup consisted of 483 participants (the three subgroups were decreased, increased, and stable desire; the grouping procedure is described further in the Statistical Analysis section), we randomly removed participants from the two bigger subgroups in order to equalize the subgroup sizes to 483. This resulted in the final total sample of 1,449 women.

Both data collections were approved by the Ethics Committee of Åbo Akademi University, and the study was carried out in accordance with the Helsinki Declaration. Written informed consent was obtained from all participants at both time points. Data can be made available upon request to the corresponding author.

Measures

A short form of the revised Sexual Desire Inventory33 (SDI-2) was used for creating the grouping variable assessing change in sexual desire over time. The following measures were included in the network: The Female Sexual Function Index34 (FSFI), assessing six domains of sexual functioning: sexual desire, sexual arousal, lubrication, orgasm function, sexual satisfaction, and sex-related pain; a short form of the Female Sexual Distress Scale35, assessing sexual distress; two subscales of the Brief Symptom Inventory-1836, assessing depression and anxiety; the body image subscale of the Derogatis Sexual Function Inventory37, assessing body dissatisfaction; the Alcohol Use Disorders Identification Test38, assessing alcohol use; the Desired and Actual Sexual Activity Scale39, assessing discrepancy in desired and actual frequency of sexual behaviours (too much sexual activity and too little sexual activity); and the Sociosexual Orientation Inventory40, assessing sociosexual orientation. Detailed information about the measures and psychometric properties is provided as supplementary material. Other measures included in the networks were age, height, weight, use of hormonal contraception (yes/no), number of (biological) children, and a question inquiring about whether the woman was in a steady sexual relationship (relationship status) at the time of study (yes/no).

Statistical Analyses

Handling of missing data

At T1, 4.1% of the data were missing in the included variables, and at T2, 0.2% of the data were missing. We imputed missing data for quantitative variables with the SPSS Missing Value Analysis regression method, using residual adjustment. We used all quantitative study variables as predictors for T1 and T2, respectively.

Sexual desire groups

We based the grouping on the SDI-2 and not the FSFI desire subscale, as grouping based on a measure that is included in the network can lead to induced problematic artefacts within the covariance structure within the subgroups41. Individual change scores for the SDI-2 were first computed by subtracting the T1 score from the T2 score. Then, the women were divided into three groups based on the change score: those with (a) decreased levels of sexual desire (women with scores more than −0.5 standard deviations from the mean of the change score), (b) increased levels of sexual desire (women with scores more than +0.5 standard deviations from the mean), and (c) a relatively stable level of sexual desire (women with scores within 0.5 standard deviations from the mean). Note that variables included in the estimated networks were, in contrast to the grouping variable, based on only one time point (T1 and T2, respectively). Before the network estimations, we assured that our grouping variable did not correlate too strongly with any variables included in the estimated networks (see supplementary material), making it unlikely that the grouping procedure affected the covariance structure and induced relationships in the networks by conditioning on a common effect42.

Network analyses

We analysed the data both at T1 and T2, using network packages for R (version 3.3.3). In the analyses for T2, we aimed to replicate the analyses of T1. Each analysis was run for each desire group separately.

Network estimation and interpretation: First, we estimated the network models of the three groups and compared the network structure visually and by comparing differences in edge significance and explained variance of variables. We further checked whether edge estimates were stable across bootstrapped estimations based on subsamples of the data.

We estimated the network models with the mgm package (i.e., our networks are mixed graphical models, MGMs23, which model relationships according to the distributional assumptions of the respective variables; continuous, ordinal, or categorical). In our MGMs, all relationships represent pairwise interactions (i.e., k = 2, interactions). Furthermore, in all estimations of relationships between two variables, their relationships with all other variables included in the network are controlled for. Thus, the absence of a relationship between two variables indicates that those two variables are conditionally independent given all other variables. In order to minimize the number of estimated parameters and limit the likelihood of estimating false positives, we used regularization and EBIC model selection (i.e., the software estimates several models with differing levels of sparsity and chooses the best fitting model according to the Extended Bayesian Information Criterion, which is similar to the Bayesian Information Criterion but with an additional term that takes into account the size of all possible models), and set the hyperparameter to γ = 0.5 as suggested in the literature43. Further technical details about MGMs can be found elsewhere23.

We visualized the estimated network models using the qgraph package44. The network layout was determined with the Fruchterman-Reingold algorithm for each network separately, so that strongly connected nodes attract each other whereas disconnected nodes repulse each other44. However, to make the networks visually comparable we used the average layout across the three individual layouts when plotting the network models.

In the interpretation of relationships and respective differences between groups, we also used edge stability and significance plots (bootnet package22). Edge stability plots indicate how accurately we can estimate edges as well as how stable the order of relative edge magnitudes is based on the data at hand. In the edge stability plot, we visualized 95% bootstrapped confidence intervals based on re-estimations of edges on resampled observations in the data (with replacement). Note that edge significance plots do indicate which edges are significantly different from each other in magnitude (α = 0.05), but not whether edges are significantly different from 0.

We also calculated how much of a variable’s variance can be explained by variables connected to it in the network (i.e., nodewise predictability45). High predictability of a variable indicates that most of the variance of that variable can be predicted by the variables it is directly connected to.

Network connectivity comparison: Second, we compared the connectivity of the three group networks formally, using the three tests included in the NetworkComparisonTest46 (NCT): the test for differences in global strength (i.e., do the networks differ in overall strength of all relationships?), structure invariance (i.e., are there any significant differences in the nature and strength of individual edges?), and, in case the networks were structurally variant, the test for differences of each edge (i.e., which individual edges show significant differences across groups?). We should note that the NCT does not estimate relationships among variables exactly like the MGM. However, the NCT is the only statistical network comparison method currently available.

Node centrality: Third, we explored whether there were differences in the centrality of variables (i.e., nodes) between the three groups. We planned to retrieve the three most popular centrality statistics: betweenness, closeness, and strength centrality. However, since bootstrapped stability analyses indicated that the closeness and betweenness estimates were too instable across subsamples of the data (see supplementary material), we limited the interpretation of centrality statistics to strength centrality (i.e., the number and strength of a node’s direct relationship with other nodes). In the interpretation of centrality, we also used bootstrapped strength centrality statistics and strength centrality significance tests (i.e., we checked whether nodes were central across subsamples of the data and whether they were significantly more central than other nodes, which is the case when a node is more central than other nodes in 95% of bootstrapped subsamples).

Network clusters: Lastly, we ran a cluster detection algorithm to explore differences in the structure of connectivity across the three groups from an additional perspective. Clusters of nodes represent more strongly connected subnetworks in the larger network. A cluster can pinpoint a group of nodes that were to be affected most quickly when a node included in the respective subnetwork changes states. We used the walktrap algorithm, which identifies clusters of nodes through random walks across the network connections (igraph package47). We ran several estimations with increasing numbers of steps and chose the number of steps which resulted in the first stable number of clusters (see supplementary material for further detail).