Participants

The source population for this study consisted of all individuals born in Sweden between January 1, 1980 and December 31, 1999 who had not died or emigrated before the end of the follow-up on December 31, 2013. We extracted data from the Swedish Medical Birth register, the Multi-Generation Register, the National Patient Register, and the National Crime Register. All registers were linked via the unique personal identification number assigned to each Swedish resident at birth.

We identified two samples. The first sample included the oldest full-sibling pair within each family (N = 580,891 pairs), with a mean age of 24.1 years (SD, 5.1; range, 14.1–34.0) at the end of the follow-up. The second sample consisted of 22,682 twin pairs, and after excluding 5512 pairs without zygosity information, the final sample comprised 17,170 pairs, including 5133 monozygotic (MZ) and 12,037 dizygotic (DZ) twin pairs. Zygosity was determined by being of opposite sex, DNA information, or a validated algorithm based on five questions concerning twin similarity (with a probability of correct classification ≥95%) [32]. The mean age of this twin sample at the end of the follow-up was 22.5 years (SD, 5.5; range, 14.1–34.0).

This study was approved by the Regional Ethical Review Board in Stockholm, Sweden. Informed consent was obtained from the twin sample but was not required for de-identified register data by law.

Measures

We derived the p factor from the following 10 diagnoses assigned by psychiatrists after contact with the in- or outpatient psychiatric services: anxiety spectrum disorder (anxiety, obsessive-compulsive disorder, and/or post-traumatic stress disorder), depression, bipolar disorder, eating disorder, drug misuse, alcohol abuse, attention deficit hyperactivity disorder (ADHD), autism, tics, and schizophrenia (containing schizoaffective disorder). Supplementary Table 1 presents related International Classification of Diseases (ICD) codes.

Exposure

The exposure was older siblings’ observed total diagnostic sum score, which served as a proxy for the latent p factor. We turned the sum score into binary dummy codes, whereby each p sum score value was compared to a reference group with p sum score equal to 0 (i.e., 0 vs. 1; 0 vs. 2; etc.). The dummy-coding allowed for examining if the associations between the siblings increased in a linear fashion, even at very high scores (i.e., it allowed for investigating potential non-linearity).

Outcome

The outcome was the younger siblings’ observed total diagnostic sum score. To examine how associated the observed diagnostic sum score was with the corresponding latent p factor, we estimated its reliability, that is, how much variance in the sum score was accounted for by the latent p factor.

To derive the latent p factor, we applied exploratory structural equation modeling (ESEM) to the 10 psychiatric diagnoses [33]. We decided on the number of factors to extract based on scree plot [34], and then rotated the factors toward one general and several uncorrelated specific factors using the Direct Schmid–Leiman transformation [35]. This way, the general factor (p factor) captured the shared variance among all psychiatric diagnoses, whereas the specific factors captured the variance unique to subsets of psychiatric disorders over and above the p factor.

Because the factor indicators were binary diagnoses, we used Item Response Theory (IRT) to estimate how much variance in the total sum score was accounted for by the latent p factor (i.e., its reliability). IRT reliability estimates differ in two ways from those based on Classical Test Theory (which is suitable for continuously distributed factor indicators). First, IRT reliability is conditional on the latent score (e.g., reliability could be high for individuals who are above the latent mean, but low for individuals who are below the latent mean). Second, IRT reliability estimates are usually expressed in a scale-dependent fashion (unlike classical reliability estimates that are commonly expressed as a scale-free R2). To facilitate interpretability, we translated the scale-dependent IRT reliability estimate into a conditional Classical Test Theory estimate, such that the conditional reliability was expressed as a scale-free R2 [36]. An R2 above 0.70 (i.e., that the latent factor accounted for at least 70% of the variation in the corresponding sum score) is generally considered acceptable [37].

To ensure that the sum score of the younger and older siblings captured the same underlying construct, we tested whether the factor loadings were invariant in the younger and older siblings in two ways. First, we fit the aforementioned latent factor model within a two-group model framework (with one group for the younger siblings, and one for the older siblings), in which we allowed the latent factor loadings to vary between groups versus being constrained to equality. We then compared the difference in model fit (using the Comparative Fit Index, CFI, and Root Mean Square Error of Approximation, RMSEA) between the more constrained (i.e., where the loadings were constrained to equality) versus less constrained model (i.e., the model where the loadings were allowed to vary between groups). Based on simulations, Cheung and Rensvold recommended that a ΔCFI < 0.01 was inadequate to conclude that two nested models differed [38]. Second, using the less constrained model in which the loadings were allowed to vary, we examined the similarity in the factor loadings by computing the factor congruence coefficient, with values above 0.95 implying that two factors can be considered equal [39].

Statistical analysesEstimating the association between the exposure and outcome

We regressed the younger siblings’ p sum score onto the older siblings’ dummy-coded p sum score. As the exposure was binary (e.g., 0 vs. 1; 0 vs. 2; etc.), the ensuing betas correspond to mean differences in the younger siblings’ p sum score for each additional diagnoses in the older sibling. In addition to visually examining whether the associations appeared linear into the extreme, we also conducted a linear-by-linear trend test. This test is more suitable than adding a quadratic term in the regression when the exposure is a categorical variable. A significant p-trend value rejects the null hypothesis that the trend is non-linear [40]. All regressions included the younger siblings’ age as a covariate.

Negative control analysis

Past research has shown that whereas mild intellectual disability runs in families, severe intellectual disability seldom does (presumably because it is primarily caused by rare mutations or environmental factors such as traumatic brain injury unique to only one sibling) [27, 41]. Given that past research has shown that the p and g factors are inversely associated [42,43,44,45], if the p factor were mainly attributed to common genetic variants, then one might expect that it should be associated with mild but not with severe and profound intellectual disability. Therefore, as a negative control condition, we examined the familial coaggregation between the p factor and diagnoses of intellectual disability of varying degrees of severity. Specifically, we regressed younger siblings’ p factor onto the older siblings’ intellectual disability diagnosis, where mild (2–3.33 standard deviations [SD] below the g factor mean), moderate (3.33–4.33 SD below the g factor mean), and severe-profound (>4.33 SD below the g factor mean) intellectual disability were compared to a reference group without intellectual disability [27, 28, 46].

DF extremes analysis and twin heritability

We first computed the observed p sum scores for both twins. We then used a DF extremes analysis and a classical twin model to estimate the group heritability and individual differences heritability of these p sum scores, respectively. The DF analysis tests whether extreme and normal variations in the p factor are genetically linked [29,30,31]. We detail this approach in the Supplementary Method. Briefly, estimating group heritability involves identifying twins who score above a cut-off (i.e., probands), and then estimating the degree to which the means of their co-twins regress toward the population mean. If the mean of DZ co-twins regress further to the population mean than that of MZ co-twins, it would imply that p sum scores both above and below the specified cut-off are genetically linked. To facilitate comparison with the g factor, we used the same cut-offs as those for intellectual disability to define the proband groups, namely mild (2–3.33 SD above the p sum score mean), mild-profound (>2 SD above the p sum score mean), severe-profound (>4.3 SD above the p sum score mean), and profound (> 5.33 SD above the p sum score mean). In addition, if the group heritability estimates (hg2) are similar to those of the individual differences heritability (h2), this further suggests that p sum scores above and below the specified threshold likely have the same etiology [29, 30]. Therefore, we also applied the classical twin model to decompose the variance of the p sum score into additive genetic effects (A), shared environment effects (C), and nonshared environment effects (E) [47], and compared the individual differences heritability to the group heritability estimates.

Sensitivity analyses

We conducted five sensitivity analyses to examine the robustness of the findings. First, we regressed the younger siblings’ latent p onto the older siblings’ dummy-coded observed p sum score (see Supplementary Fig. 1 for model diagram). The advantage of this approach was two-folded. Measurement error in the outcome can generate larger standard errors. As the latent factor model is estimated to have perfect reliability, this could lead to smaller standard errors. In addition, given the multidimensional nature of the psychiatric conditions, the sum score is likely not only associated with the latent p factor, but also with the specific factors to a smaller degree. In contrast, the latent p is fixed to be uncorrelated with the specific factors, such that the association between the observed p sum score and a latent p cannot be confounded by variance attributed to specific psychopathology factors.

Second, we performed a modified familial coaggregation analysis by expanding the number of conditions used to derive the p factor from 10 to 15. The motivation for this sensitivity analysis was to allow for a more fine-grained measurement model. However, the downside was that the number of indicators for each specific psychiatric factor was uneven. In particular, there were more indicators for the internalizing factor, such that some might end up with a high p sum score by having several anxiety-related diagnoses. Specifically, we decomposed the anxiety spectrum disorder into three separate diagnoses (anxiety, obsessive-compulsive disorder, and post-traumatic stress disorder), separated schizoaffective disorder from the schizophrenia diagnosis, and included oppositional defiant disorder and court convictions of violent and/or property crimes (e.g., homicide and theft) [48] to capture a broader range of externalizing behaviors. The ICD codes for the additional psychiatric diagnoses can be found in Supplementary Table 2.

Third, to examine whether the associations between the p sum scores across family members might be impacted by rare deleterious mutations or severe environmental factors such as traumatic brain injury, we excluded sibling and twin pairs in which at least one member of each pair had diagnoses of severe or profound intellectual disability, and then we re-ran the familial coaggregation analyses and DF extremes analysis.

Fourth, aside from the p sum score, we also used p factor scores as exposures and outcomes. Whereas sum scores create a scale by applying unit weights to each indicator (e.g., indicator 1, 2, 3, etc., are simply summed into a scale), factor scores allow the weights to vary (e.g., indicator 1 might contribute 0.5 units, indicator 2 might contribute 0.75 units, etc., to the scale score). Both approaches have their respective advantages [49,50,51,52].

Fifth, as some disorders might have a later age of onset, we re-ran the models in a subsample in which the participants were 28–34 years old.

Data were analyzed from February 2022 to December 2022 using software SAS 9.4 [53], Mplus 8.3 [54], and R 4.0.5 [55] with GPArotation [56] package.