The Power to Resolve Cultural Transmission and Sibling Interaction Using Polygenic Scores

doi:10.21203/rs.3.rs-5461256/v1

In the classical twin design, the assumption that the additive genetic (A) and shared environment (C) variance components are uncorrelated may not hold. If there is positive AC covariance, the C component is overestimated. Two processes have been studied that lead to AC covariance: Cultural transmission (e.g., genetic nurture), when the parents’ genotype contributes to the effective environment of the child, and sibling interaction, when the genotype of one sibling contributes to the effective environment of another. Several designs use polygenic scores of parents or siblings to detect AC covariance, but these models cannot unambiguously identify the source. A combined model has been proposed, but its power to identify both processes has not been well-studied yet. This study uses simulated data from genotyped twins and their parents to investigate the power to disentangle these processes. Results demonstrated that we can detect AC covariance using either genotyped-sibling or genotyped-parent data, but we cannot resolve its source and risk making wrong inferences. However, these sources of AC covariance can be resolved using genotyped data of both siblings and parents. This emphasizes the need for whole-family genotyping and modeling.

AC covariance

classical twin design

gene-environment interplay

polygenic scores

statistical power

Integrating polygenic scores (PGSs) into traditional family and twin designs has created new possibilities to address these designs’ assumptions (McAdams et al. 2023). In the classical twin design (CTD), a standard assumption is that the covariance between the additive genetic factors (A) and the shared-environment factors (C) is zero (cov[AC] = 0). This assumption is at odds with theories in social science concerning, for instance, cognitive development and educational attainment, which posit substantive gene-environment correlation (e.g., Plomin & von Stumm 2018; Sauce & Matzel 2018). The consequence of assuming cov(AC) = 0 in the CTD depends on its sign and its magnitude. Positive unmodeled AC covariance results in an overestimation of C and underestimation of A variance (Eaves 1976; Purcell 2002; Verhulst et al. 2017).

Two well-known processes within the family lead to AC covariance: cultural transmission, and sibling interaction. First, in the case of cultural transmission, parental heritable behavior contributes to the offspring’s environment relevant to the phenotype of interest. This is known as genetic nurture (Kong et al. 2018). However, the term genetic nurture suggests positive parental effects, as in a stimulating nurturing environment, but cultural transmission also includes negative parental effects, for instance parental aggression. To assess cultural transmission, the CTD has been extended to include additional family members, such as the parents of the twins in the nuclear twin family design (Fulker 1988; Keller et al. 2009) or the offspring of twins in the children-of-twins design (D’Onofrio et al. 2003; Eaves et al. 2005; McAdams et al. 2018). Second, in the case of sibling interaction, the behavior of one twin contributes to the effective environment of the other twin. Such sibling interaction has also been included in the twin design (Carey 1986; Dolan et al. 2014; Eaves 1976; Eaves et al. 1977). If the siblings’ behavior has a heritable component, this gives rise to AC covariance. Again, if present, but unmodeled, positive AC covariance will result in an overestimation of the C and underestimation of the A variance components.

One approach to detect AC covariance is the use of PGSs (Dolan et al. 2021). This approach uses single-nucleotide polymorphisms (SNPs) aggregated into an individual’s PGS to offer a summary measure for someone’s genetic tendencies for a particular trait. One proposed design, in which the offspring phenotype is regressed on both the parental alleles that were transmitted, and the parental alleles that were not transmitted to the offspring, is the transmitted/non-transmitted (T/NT) design (Bates et al. 2018; Kong et al. 2018). An effect of the non-transmitted alleles suggests the presence of cultural transmission. Later studies demonstrated that the T/NT regression model is equivalent to the regression model in which the offspring phenotype is regressed on the offspring PGS and the parental PGSs (Okbay et al. 2022; Tubbs et al. 2020). Conversely, an effect of the parental PGSs (beyond the effect of the offspring’s PGS) suggests the presence of cultural transmission. This approach has the advantage that it does not require alleles to be identified as transmitted or non-transmitted.

Another approach, which exploits genetic and phenotypic information in dizygotic (DZ) twins (or full siblings), is within-family regression (Selzam et al. 2019). This approach involves regressing the individual’s phenotype on the individual’s PGS and the DZ-pair mean PGS. If the effect of the DZ-pair mean is significant, this means that an individual’s phenotype is predicted by the PGS of their sibling (beyond the effect of their own PGS), demonstrating AC covariance. This effect is usually attributed to sibling interaction, but this attribution may be wrong.

The detection of AC covariance per se using either regression method does not necessarily identify the process that gave rise to the AC covariance. At first glance, the T/NT design and the equivalent parental PGS regression model seem to serve the aim of detecting cultural transmission. However, observed AC covariance may not result from the parents but rather from the siblings. Similarly, the within-family regression method (Selzam et al 2019) uses sibling data to detect AC covariance but identified AC covariance does not necessarily stem from sibling interaction. Another study combined the use of parental and sibling PGSs, but did not study statistical power (Young et al. 2022; see also Demange et al. 2022; Tubbs et al. 2020). This combined model (i.e., combining parental and sibling PGSs) allows to dissect AC covariance into the part due to cultural transmission and the part due to sibling interaction. The present aim is to apply this model to explore the power to detect and resolve AC covariance stemming from cultural transmission and sibling interaction in the CTD (monozygotic [MZ] and DZ twins) and in the DZ (or full siblings) only design. We formulated the models as structural equation models (SEM) and calculated the statistical power of the relevant likelihood ratio tests, given various parameter settings and predictive power of the PGSs. We also considered the formulations of the models as linear regression models, using generalized estimating equations (GEE), as this method is easier to implement than SEM.

The outline of this paper is as follows. We first present the method in terms of the relevant regression equations and the corresponding SEM. Subsequently, we present the simulation design, and the results of applying this design to address the issues of resolution and power. We end this paper with a discussion of results and limitations.

The regression model used here combines the two models of Okbay et al. (2022) and of Selzam et al. (2019). The model by Okbay regresses the child’s phenotype on the child’s PGS and parental PGSs. Disregarding covariates, the model can be expressed as:

(1) ph_ik = b₀ + b₁ × PGS_ik + b₂ × (PGS_im + PGS_if ) + ε_ik,

where ph_ik denotes the phenotypic score in child k in family i, PGS_ik, PGS_im, and PGS_if are the PGSs in the child, the mother and the father, respectively, and ε_ik is the residual. In this model, the effects of the father's and the mother's PGS are equal (i.e., b₂). If desired, this assumption can be relaxed. In within-family regression, the phenotypic variance is decomposed into a between-family effect, represented by the family mean, and a within-family effect, represented by the deviation of the individual PGS from the family mean. We used an equivalent, but slightly altered, model: Instead of the deviation from the family mean, we used the raw PGS of the individual to represent the within-family genetic effect. In the presence of genotyped siblings, we regressed twin phenotype on their own PGS and the average PGS of them and their co-twin, such that:

(2) ph_ik = f₀ + f₁ × PGS_ik + f₂ × PGS_i. + ε_ik ,

where PGS_i. denotes the average PGS of the twin pair.

Given Eq. 1, the rejection of b₂ = 0 suggests that cov(AC) ≠ 0, which is often attributed to cultural transmission (Balbona et al. 2021; Tubbs et al. 2020). Given Eq. 2, rejection of f₂ = 0 suggests that cov(AC) ≠ 0, but does not necessarily identify the source (or process) of AC covariance. A significant f₂, so an effect of the sibling mean, could be due to either cultural transmission or sibling interaction. Conversely, the rejection of b₂ = 0 could be due to sibling interaction, rather than cultural transmission. We investigated this possibility of misattribution using the regression model that combines PGSs of parents and offspring and the phenotypes of the offspring (Demange et al. 2022; Young et al. 2022):

(3) ph_ik = g₀ + g₁ × PGS_ik + g₂ × PGS_i. + g₃ × (PGS_im + PGS_if) + ε_ik ,

Below, we used exact data simulation (see van der Sluis et al. 2008) according to the model, where the residual ε_ik is bivariate normally distributed N(0, S_ε²). The 2x2 covariance matrix S_ε² is a function of residual genetic and (shared and unshared) environmental influences relevant to the phenotype. In this model, the test of AC covariance is based on the parameters g₂ and g₃. In our simulation g₃ is attributable to the process of cultural transmission, while g₂ is attributable to the process of sibling interaction. We do not consider the main effects of covariates, such as genetic principal components (Price et al. 2006) sex, or age, as these are not relevant to the present issue.

The estimation of the parameters in these regression models can be done using analyses that can account for correlated residuals (e.g., S_ε²). We can use GEE regression analysis, linear mixed modeling, or SEM. In all three approaches, the regressors are treated as fixed regressors, such that the distributional assumption concerns the distribution of the residuals (ε_ik). However, in SEM, we can also choose to treat the regressors as random variables. Treated as such, the distributional assumption concerns the joint distribution of the PGS and the phenotypes, as in this approach a covariance structure model is fitted to the PGS and the phenotypes, and the PGS and phenotypes are assumed to be multivariate normally distributed. Using the SEM approach has the considerable practical advantage that it can handle missing PGS and missing phenotypic data (i.e., using full information maximum likelihood estimation). We note that the random/fixed terminology is confusing. Therefore, we emphasize that the treatment of the regressors as random variables does not imply random effects: the effects of the random regressors are fixed (to be estimated) effects. Random effects are relevant in multilevel modeling, in general, and with linear mixed modeling, particularly in genetics (e.g., Yang et al. 2011).

Procedure

Data simulation was based on the model presented in Fig. 1. We limited the analyses to one set of values of the three variance components (ACE) but varied the predictive power of the PGS (i.e., the proportion of A variance that is accounted for by the PGS), and the effect sizes for cultural transmission, and sibling interaction. Explicitly, we first simulated data for a single set of variance components (a² = .4, c² = .3, e² = .3), one sample size (n_MZ = n_DZ = 4,000 pairs), and set the PGS to explain 10% of A variance. We varied the two sources of AC covariance (settings: 0, 0.05, 0.1) in a 3x3 design. Note that the variance components sum up to 1 (i.e., .4 + .3 + .3), but the total phenotypic variance equals a² + c² + e² + 2 · cov(AC), meaning it will exceed 1 in the presence of positive AC covariance. Subsequently, we varied the predictive power of the PGS (settings: 2%, 5%, 10%, 15%) and conducted analyses for a DZ-only sample (n_DZ = 4,000 pairs). The predictive power is defined in terms of a percentage of total additive genetic variance. Below we present the effect size in terms of both PGS predictive power and in explained phenotypic variance. We note that researchers may be interested in a broader range of settings for the ACE variance components. While these are not the focal point of the current study, we consider a broader range of settings in the supplement (Appendix A).

We adopted the SEM approach using OpenMx. We present results for three OpenMx models: The combined model based on Eq. 3, including the parents’ and twins’ PGSs, the cultural transmission-only model based on Eq. 2 (nested under the Eq. 3 model by fixing g₃ to zero) and the sibling interaction-only model based on Eq. 1 (nested under the Eq. 3 model by fixing g₂ to zero). We focused on how the models that only assume one source of AC covariance performed in comparison to the combined model in terms of power and effect sizes. Hereafter, we refer to the coefficient for cultural transmission as g and the coefficient for sibling interaction as b.

We fitted the models to DZ and MZ twin data, i.e., data typically present in the twin and family registries, but we note that MZ twin data are not required to fit these models. The residual covariance matrix S_ε² is expected to differ in the MZ and DZ samples, due to residual genetic influences (i.e., not explained by the PGS). We also fitted the models to the DZ twin data only, as in practice, phenotypic data may be limited to full siblings.

Power calculations were based on the likelihood-ratio test statistic, given the degrees of freedom of the test, the non-centrality parameter of the χ² distribution, and the a of 0.05 (e.g., Hewitt & Heath 1988). We also present estimated path coefficients for cultural transmission and sibling interaction. We express the effect size in terms of the increase in phenotypic variance in MZ and DZ twins that is due to AC covariance (in comparison to the baseline ACE model with no AC covariance). This measure is zygosity-dependent in the presence of sibling interaction and will then be higher in MZ twins, assuming positive AC covariance (Carey 1986; Eaves 1976; Eaves et al 1978; see Appendix C). We also present results for four different setting of the predictive power of the PGS, which is expressed as the percentage of A variance captured by the PGS. From these values, we can calculate the proportion of overall phenotypic variance that is captured by the PGS. In absence of AC covariance (g = b = 0), and given a² = .4 and the values for the predictive power of the PGS (settings: 2%, 5%, 10%, 15%), the PGS captures 0.8%, 2%, 4%, and 6% of the total phenotypic variance, respectively. This value increases in the presence of AC covariance and is, again, zygosity-dependent in the presence of sibling interaction. For example, when g = .1 and b = .1, the PGS captures 4.9% of MZ phenotypic variance and 5.2% of DZ phenotypic variance (see Appendix C).

The SEM model can be extended by fitting an ACE model to the residual (background) covariance matrices S_ε² in the MZ and DZ twin data. Fitting the extended model is not necessary for the purpose of this study, but it is relevant to note that this background model is misspecified in the presence of AC covariance. The estimate of the residual C variance is expected to be biased in the presence of AC covariance, because the PGS only captures part of the A variance, and thus only part of the AC covariance. For reference, we present extended results for the SEM approach in a wider range of settings (Appendix A). We also provide results of the GEE regression analyses (Appendix B).

Materials

All analyses were performed in R (R Development Core Team 2023) used in RStudio (RStudio Team 2023). We used the mvrnorm() function of the MASS package (Ripley et al 2023) for exact data simulation, geepack (Højsgaard et al 2022) for GEE regression (Appendix B), and OpenMx (Neale et al 2016) for SEM. For ease of use, analysis steps were combined in the package gnomesims (Bernardo 2024; see Appendix D).

We present power results for the combination of the 3x3 settings for cultural transmission and sibling interaction (Table 1). The power was calculated based on the results of the three OpenMx analyses (see m1-m3, Table 1) of simulated data of MZ and DZ twins and assuming a PGS with a predictive power of 10%. For detailed results and DZ-only regression see Appendix A. As expected, the power to reject b = 0 and/or g = 0 is equal to the alpha (α = .05), if these parameters are truly zero (i.e., cov[AC] = 0). In presence of a single source of AC covariance, only the combined model results in correct inference concerning the source of AC covariance. Notably, when we apply a model that incorrectly assumes the presence of a single source of AC covariance that is not actually present in the dataset, the power exceeds 0.05 - which could result in the erroneous conclusion that sibling interaction is present when, in reality, only cultural transmission is present, or vice versa. The extent to which power exceeds .05 depends on the confounding source of AC covariance. For the smaller setting (g = 0.05 or b = 0.05) one will only detect a false effect in around 1 in 6 cases, while for the larger setting (g = 0.1 or b = 0.1) one will detect a false effect in half of the cases.

Table 1 Power Estimates in Factorial Design With Two Factors (CT, SI) Each With 3 Levels (9 Cells).^a

CT	SI	CT(m1)	SI(m2)	CT(m3)	SI(m3)	MZ s²_Ph	DZ s²_Ph
0	0	0.05	0.05	0.05	0.05	0.00%	0.00%
0.05	0	0.408	0.161	0.3	0.05	6.82%	6.82%
0.1	0	0.917	0.474	0.787	0.05	14.65%	14.65%
0	0.05	0.16	0.399	0.05	0.292	6.57%	3.41%
0	0.1	0.505	0.929	0.05	0.8	13.65%	7.32%
0.05	0.05	0.756	0.754	0.285	0.281	13.9%	10.74%
0.1	0.05	0.989	0.945	0.761	0.271	22.22%	19.06%
0.05	0.1	0.953	0.992	0.27	0.783	21.47%	15.15%
0.1	0.1	0.999	0.999	0.734	0.763	30.3%	23.97%

^aCT = cultural transmission, SI = sibling interaction, m1 = model only assuming CT, m2 = model only assuming SI, m3 = model assuming both CT and SI, MZ s²_Ph = percentage of variance inflation due to AC covariance for MZ twins, DZ s²_Ph = percentage of variance inflation due to AC covariance for DZ twins. We have marked instances of only one source of AC covariance in bold.

The results obtained using GEE regression were very similar to the results in OpenMx (Appendix B) and may, therefore, be a suitable alternative to SEM. As mentioned above, GEE is arguably easier to implement than the SEM model. Nevertheless, SEM has the considerable advantage of full information maximum likelihood estimation, which can accommodate missing data.

Next, we present parameter estimates for the same combination of AC covariance settings (Table 2). We note that the presence of one source of AC covariance biases the estimate for the other source, if we apply the incorrect model. Only when applying the combined model do the coefficients equal the parameters used for data simulation. The same holds true in the presence of two sources of AC covariance – when using a model that only assumes one source, the influence of that source is overestimated, when the other source is not considered (Fig. 2).

Table 2 Parameter Estimates in Factorial Design With Two Factors (CT, SI) Each With 3 Levels (9 Cells).^b

CT	SI	CT (m1)	SI (m2)	CT (m3)	SI (m3)	MZ s²_Ph	DZ s²_Ph
0	0	0	0	0	0	0.00%	0.00%
0.05	0	0.05	0.029	0.05	0	6.82%	6.82%
0.1	0	0.1	0.058	0.1	0	14.65%	14.65%
0	0.05	0.027	0.05	0	0.05	6.57%	3.41%
0	0.1	0.058	0.1	0	0.1	13.65%	7.32%
0.05	0.05	0.078	0.079	0.05	0.05	13.90%	10.74%
0.1	0.05	0.129	0.108	0.1	0.05	22.22%	19.06%
0.05	0.1	0.11	0.129	0.05	0.1	21.47%	15.15%
0.1	0.1	0.162	0.157	0.1	0.1	30.30%	23.97%

^b CT = cultural transmission, SI = sibling interaction, m1 = model only assuming CT, m2 = model only assuming SI, m3 = model assuming both CT and SI, MZ s²_Ph = percentage of variance inflation due to AC covariance for MZ twins, DZ s²_Ph = percentage of variance inflation due to AC covariance for DZ twins. We have marked instances of only one source of AC covariance in bold.

Effect of PGS Predictive Power

We looked at the effect of different PGS predictive power settings on statistical power to detect AC covariance (Fig. 3). Here, the focus lies on power for the regressions that only account for one source of AC covariance (see Table 1, bold values). The values show that larger AC covariance effect sizes and higher PGS predictive power increase the chance of a false positive (Type I error). This necessitates using the combined model – while larger AC covariance effects and a higher PGS increases the chance of true findings, they also increase the chance of finding an effect for the incorrect source of AC covariance if the applied model does not match the data-generating model.

DZ-only (Full Siblings) Analyses

Finally, we present OpenMx results of the analyses given DZ-only twin (i.e., full sibling) data. We set the sample size to n_DZ = 4000 pairs, and again used exact data simulation. Compared to the MZ& DZ sample, patterns are very similar (see Table 3 & 4). The power to detect cultural transmission is higher than the power to detect sibling interaction, for both samples. However, the discrepancy between the power to detect cultural transmission and the power to detect sibling interaction is larger in the DZ-only sample. Notably, the inflation of the overall variance due to AC covariance is lower than in the mixed sample given the presence of sibling interaction.

Table 3 Power Estimates in Factorial Design With two Factors (CT, SI) Each With 3 Levels (9 Cells) in a DZ-only Sample.^c

CT	SI	CT (m1)	SI (m2)	CT (m3)	SI (m3)	DZ s²_Ph
0	0	0.05	0.05	0.05	0.05	0.00%
0.05	0	0.152	0.05	0.261	0.156	6.82%
0.1	0	0.424	0.05	0.724	0.453	14.65%
0	0.05	0.05	0.171	0.179	0.302	3.41%
0	0.1	0.05	0.506	0.55	0.819	7.32%
0.05	0.05	0.146	0.163	0.64	0.651	10.74%
0.1	0.05	0.404	0.156	0.94	0.894	19.06%
0.05	0.1	0.14	0.48	0.915	0.962	15.15%
0.1	0.1	0.384	0.454	0.994	0.995	23.97%

^c CT = cultural transmission, SI = sibling interaction, m1 = model only assuming CT, m2 = model only assuming SI, m3 = model assuming both CT and SI, MZ s²_Ph = percentage of variance inflation due to AC covariance for MZ twins, DZ s²_Ph = percentage of variance inflation due to AC covariance for DZ twins. We have marked instances of only one source of AC covariance in bold.

Table 4 Parameter Estimates in Factorial Design With two Factors (CT, SI) Each With 3 Levels (9 Cells) in a DZ-only Sample.^d

CT	SI	CT (m1)	SI (m2)	CT (m3)	SI (m3)	MZ s²_Ph	DZ s²_Ph
0	0	0	0	0	0	0.00%	0.00%
0.05	0	0.05	0.033	0.05	0	6.82%	6.82%
0.1	0	0.1	0.067	0.1	0	14.65%	14.65%
0	0.05	0.039	0.05	0	0.05	6.57%	3.41%
0	0.1	0.079	0.1	0	0.1	13.65%	7.32%
0.05	0.05	0.089	0.083	0.05	0.05	13.90%	10.74%
0.1	0.05	0.14	0.117	0.1	0.05	22.22%	19.06%
0.05	0.1	0.131	0.133	0.05	0.1	21.47%	15.15%
0.1	0.1	0.182	0.167	0.1	0.1	30.30%	23.97%

^d CT = cultural transmission, SI = sibling interaction, m1 = model only assuming CT, m2 = model only assuming SI, m3 = model assuming both CT and SI, MZ s²_Ph = percentage of variance inflation due to AC covariance for MZ twins, DZ s²_Ph = percentage of variance inflation due to AC covariance for DZ twins. We have marked instances of only one source of AC covariance in bold.

We present the comparison between the MZ & DZ sample and the DZ-only sample for the OpenMx analyses in Fig. 4. The estimates are relatively similar, but the measures for the DZ-only sample differ more from each other and the power to detect sibling interaction is slightly lower than the power to detect cultural transmission.

In this simulation study, we studied the power to disentangle two sources of AC covariance, namely cultural transmission and sibling interaction, using PGSs of genotyped twins and their parents. We demonstrated that a model that only assumes one source of AC covariance might be used to detect the presence of AC covariance but cannot unambiguously identify the process underlying the AC covariance. The model that incorporates both sources of AC covariance offers the advantage, in principle, of correctly identifying both sources.

Strengths and Caveats

Our analyses provide a reference for effect sizes needed to detect cultural transmission and sibling interaction, which allows researchers to make informed decisions with respect to their research design, sample size, and assessment of whether the available PGSs of the phenotype, and the heritability (see Appendix A), are sufficient to detect an effect. We considered both MZ and DZ twin data, and DZ-only data. DZ-only data mimics full-sibling data (assuming AC covariance processes are identical in DZ twins and full siblings). The R package accompanying this paper (named gnomesims) can be used to do power analyses using both SEM and GEE-based GLM. We demonstrate that the results for GEE are similar to those obtained using SEM, which makes GEE a suitable alternative to carry out power analyses and to fit the models. Therefore, our models may be of use to researchers from different fields and applied to a variety of phenotypes. Compared to GEE, we note that the random regressor implementation in SEM has the considerable advantage that one can obtain maximum likelihood estimates and likelihood ratio tests in the presence of missing data, using full information maximum likelihood estimation (Arbuckle 1996).

An important caveat of the current study is that the measures of cultural transmission and sibling interaction, respectively, might be biased by other sources of AC covariance. When the combined model is applied to real data, effects may therefore still reflect other sources of AC covariance. This means that we may not be able to unambiguously interpret the results of the combined model either. There are many examples of other sources – among others, dynastic social processes in the extended family (Nivard et al 2024), geographical clustering (Tamimy et al 2021), or assortative mating (Sunde et al 2024).

Implications & Future Directions

This paper builds on existing studies on genetic confounding in family-based designs, which have frequently used the phenotype educational attainment. Indirect genetic effects are relevant to educational attainment (Young et al 2022), but the exact nature of these effects remains unclear. While one study comparing individuals with and without siblings (Howe et al 2022) found evidence for indirect genetic effects via the sibling, another study found evidence for cultural transmission, but no evidence for sibling interaction (Demange et al 2022). Our study offers some perspective on how to contextualize these contradictory findings; It may be that sibling effects are present but too small to detect due to insufficient power.

The cultural transmission-only and the sibling interaction-only regression models can be used to detect AC covariance, but the subsequent interpretation of AC covariance in terms of a given process may not be correct. We recommend, if data resources permit it, the use of the combined model (Young et al 2022) to resolve sibling interaction and cultural transmission. For future research, we emphasize the need for further methodological studies and the application of existing models. First, other simulation studies should continue to address methodological issues, such as age-moderation of AC covariance. Second, future research should focus on the application of these models to different phenotypes and include other AC covariance mechanisms. Different mechanisms could include other passive sources of AC covariance, such as effects of the extended family, as well as evocative and active AC covariance. By accounting for mechanisms other than cultural transmission and sibling interaction we might gain a better understanding of which AC covariance processes are relevant to which phenotype. Understanding the relevant processes also relates to characterizing the mechanisms that we refer to as AC covariance, so that they can be clearly interpreted. Cultural transmission, for example, could refer to what is traditionally thought of as rearing effects, but also to the prenatal environment, the socio-economic status of the family, or other factors. A more precise characterization would help to move from mere detection to formulating recommendations and interventions for a given phenotype. One way to achieve this would be to incorporate ‘environmental’ measures - such as specific markers of the home environment, schooling, or parenting - into genetically informative designs, bearing in mind that such measures may be subject to genetic influences (Tucker-Drob & Bates 2016).

In summary, we present a simulation study and accompanying R package to address the power of detecting cultural transmission and sibling interaction in genomic family data. We emphasize the importance of using a model that matches the true AC covariance processes at play. Additionally, conclusions about the true effects need to be formulated more cautiously. Further research is required to identify practical issues in applying this model and incorporate other sources of confounding. We anticipate models that account for several sources of AC covariance to become more accessible and widely used in the future.

Author Contribution

C.V.D. conceptualised the work, J.B.B. and C.V.D. conducted the analyses, J.B.B. drafted the manuscript, E.v.B. and C.V.D. provided supervision, E.v.B. provided funding, all authors provided critical feedback and approved the final version.

Acknowledgements

This study was funded by the European Union (European Research Council, Project “InterGen”, PI: van Bergen, 101076726). Elsje van Bergen is furthermore funded by the Research Council of Norway (GenEd grant 335634), the Dutch Research Council (VIDI Talent Grant, https://doi.org/10.61686/VI.Vidi.221G.007), and she is a Jacobs Foundation Research Fellow. The position of Charlotte K. L. Pahnke is funded by the German Research Foundation (funding no. 428902522).

Data Availability

A script to replicate the data simulations conducted in this paper and load data files can be found on GitHub under https://github.com/josefinabernardo/covAC-simulation-study-2024/blob/main/data-submission.R.

Arbuckle JL (1996). Full Information Estimation in the Presence of Incomplete Data. In: Advanced Structural Equation Modeling. Psychology Press.
Balbona JV, Kim Y, Keller MC (2021). Estimation of Parental Effects Using Polygenic Scores. Behav Genet 51(3):264–278. https://doi.org/10.1007/s10519-020-10032-w
Bates TC, Maher BS, Medland SE, McAloney K, Wright MJ, Hansell NK, Kendler KS, Martin NG, Gillespie NA (2018). The nature of nurture: Using a virtual-parent design to test parenting effects on children’s educational attainment in genotyped families. Twin Research and Human Genetics 21(2):73–83. https://doi.org/10.1017/thg.2018.11
Bernardo JB (2024). gnomesims: Power calculations for genotyped family data (R package version 0.0.0.9000) [R]. https://github.com/josefinabernardo/gnomesims (Original work published 2024)
Carey G (1986). Sibling imitation and contrast effects. Behav Genet 16(3):319–341. https://doi.org/10.1007/BF01071314
de Zeeuw EL, van Beijsterveldt CEM, Glasner TJ, de Geus EJC, Boomsma DI (2016). Arithmetic, reading and writing performance has a strong genetic component: A study in primary school children. Learning and Individual Differences 47:156–166. https://doi.org/10.1016/j.lindif.2016.01.009
Demange PA, Hottenga JJ, Abdellaoui A, Eilertsen EM, Malanchini M, Domingue BW, Armstrong-Carter E, de Zeeuw EL, Rimfeld K, Boomsma DI, Van Bergen E, Breen G, Nivard MG, Cheesman R (2022). Estimating effects of parents’ cognitive and non-cognitive skills on offspring education using polygenic scores. Nat Commun 13(1):4801. https://doi.org/10.1038/s41467-022-32003-x
Dolan CV, De Kort JM, Van Beijsterveldt TCEM, Bartels M, Boomsma DI (2014). GE covariance through phenotype to environment transmission: An assessment in longitudinal twin data and application to childhood anxiety. Behav Genet 44(3):240–253. https://doi.org/10.1007/s10519-014-9659-5
Dolan CV, Huijskens RCA, Minică CC, Neale MC, Boomsma DI (2021). Incorporating polygenic risk scores in the ACE twin model to estimate A–C covariance. Behav Genet 51(3):237–249. https://doi.org/10.1007/s10519-020-10035-7
D’Onofrio BM, Turkheimer EN, Eaves LJ, Corey LA, Berg K, Solaas MH, Emery RE (2003). The role of the children of twins design in elucidating causal relations between parent characteristics and child outcomes. Journal of Child Psychology and Psychiatry 44(8):1130–1144. https://doi.org/10.1111/1469-7610.00196
Eaves LJ (1976). A model for sibling effects in man. Heredity 36(2):205–214. https://doi.org/10.1038/hdy.1976.25
Eaves LJ, Last KA, Young PA, Martin NG (1978). Model-fitting approaches to the analysis of human behaviour. Heredity 41(3):249–320. https://doi.org/10.1038/hdy.1978.101
Eaves LJ, Last KA, Martin NG, Jinks JL (1977). A progressive approach to non‐additivity and genotype‐environmental covariance in the analysis of human differences. British Journal of Mathematical and Statistical Psychology 30(1):1–42. https://doi.org/10.1111/j.2044-8317.1977.tb00722.x
Eaves LJ, Silberg JL, Maes HH (2005). Revisiting the children of twins: Can they be used to resolve the environmental effects of dyadic parental treatment on child behavior? Twin Research and Human Genetics 8(4):283–290. https://doi.org/10.1375/twin.8.4.283
Fulker DW (1988). Genetic and cultural transmission in human behavior. Proceedings of the Second International Conference on Quantitative Genetics, 318–340.
Hewitt JK, Heath AC (1988). A note on computing the chi-square noncentrality parameter for power analyses. Behav Genet 18(1):105–108. https://doi.org/10.1007/BF01067079
Højsgaard S, Halekoh U, Yan J, Ekstrøm CT (2022). Package “geepack”: Cran-R Project [Computer software]. https://vps.fmvz.usp.br/CRAN/web/packages/geepack/geepack.pdf
Howe LJ, Evans DM, Hemani G, Smith GD, Davies NM (2022). Evaluating indirect genetic effects of siblings using singletons. PLOS Genetics 18(7). https://doi.org/10.1371/journal.pgen.1010247
Keller MC, Medland SE, Duncan LE, Hatemi PK, Neale MC, Maes HHM, Eaves LJ (2009). Modeling extended twin family data I: Description of the cascade model. Twin Research and Human Genetics 12(1):8–18. https://doi.org/10.1375/twin.12.1.8
Kong A, Thorleifsson G, Frigge ML, Vilhjalmsson BJ, Young AI, Thorgeirsson TE, Benonisdottir S, Oddsson A, Halldorsson BV, Masson G, Gudbjartsson DF, Helgason A, Bjornsdottir G, Thorsteinsdottir U, Stefansson K (2018). The nature of nurture: Effects of parental genotypes. Science 359(6374):424–428. https://doi.org/10.1126/science.aan6877
Lee JJ, Wedow R, Okbay A, Kong E, Maghzian O, Zacher M, Nguyen-Viet TA, Bowers P, Sidorenko J, Karlsson Linnér R, Fontana MA, Kundu T, Lee C, Li H, Li R, Royer R, Timshel PN, Walters RK, Willoughby EA et al (2018). Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals. Nat Genet 50(8):1112–1121. https://doi.org/10.1038/s41588-018-0147-3
McAdams TA, Cheesman R, Ahmadzadeh YI (2023). Annual research review: Towards a deeper understanding of nature and nurture: Combining family-based quasi-experimental methods with genomic data. Journal of Child Psychology and Psychiatry 64(4):693–707. https://doi.org/10.1111/jcpp.13720
McAdams TA, Hannigan LJ, Eilertsen EM, Gjerde LC, Ystrom E, Rijsdijk FV (2018). Revisiting the children-of-twins design: Improving existing models for the exploration of intergenerational associations. Behav Genet 48(5):397–412. https://doi.org/10.1007/s10519-018-9912-4
Neale MC, Hunter MD, Pritikin JN, Zahery M, Brick TR, Kirkpatrick RM, Estabrook R, Bates TC, Maes HH, Boker SM (2016). OpenMx 2.0: Extended structural equation and statistical modeling. Psychometrika 81(2):535–549. https://doi.org/10.1007/s11336-014-9435-8
Nivard MG, Belsky DW, Harden KP, Baier T, Andreassen OA, Ystrøm E, van Bergen E, Lyngstad TH (2024). More than nature and nurture, indirect genetic effects on children’s academic achievement are consequences of dynastic social processes. Nat Hum Behav 8(4):771–778. https://doi.org/10.1038/s41562-023-01796-2
Okbay A, Wu Y, Wang N, Jayashankar H, Bennett M, Nehzati SM, Sidorenko J, Kweon H, Goldman G, Gjorgjieva T, Jiang Y, Hicks B, Tian C, Hinds DA, Ahlskog R, Magnusson PKE, Oskarsson S, Hayward C, Campbell A et al (2022). Polygenic prediction of educational attainment within and between families from genome-wide association analyses in 3 million individuals. Nat Genet 54(4): Article 4. https://doi.org/10.1038/s41588-022-01016-z
Plomin R, von Stumm S (2018) The new genetics of intelligence. Nat Rev Genet 19(3):148–159. https://doi.org/10.1038/nrg.2017.104
Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D (2006) Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet 38(8):904–909. https://doi.org/10.1038/ng1847
Purcell S (2002) Variance components models for gene–environment interaction in twin analysis. Twin Res 5(6):554–571. https://doi.org/10.1375/136905202762342026
R Development Core Team (2023) R (Version 4.3.1): A Language and Environment for Statistical Computing. (Version 4.3.1) [Computer software].
Ripley B, Venables B, Bates DM, Hornik K, Gebhardt A, Firth D (2023) Package “mass”: Cran-R Project [Computer software]. http://cran-r.c3sl.ufpr.br/web/packages/MASS/MASS.pdf
RStudio Team (2023) RStudio: Integrated Development Environment for R (Version 2023.09.0+463) [Computer software].
Sauce B, Matzel L (2018) The Paradox of Intelligence: Heritability and Malleability Coexist in Hidden Gene-Environment Interplay. Psychol Bull 144. https://doi.org/10.1037/bul0000131
Selzam S, Ritchie SJ, Pingault J-B, Reynolds CA, O’Reilly PF, Plomin R (2019) Comparing Within- and Between-Family Polygenic Score Prediction. Am J Hum Genet 105(2):351–363. https://doi.org/10.1016/j.ajhg.2019.06.006
Sunde HF, Eilertsen EM, Torvik FA (2024) Understanding Indirect Assortative Mating and Its Intergenerational Consequences (p. 2024.06.21.600029). bioRxiv. https://doi.org/10.1101/2024.06.21.600029
Tamimy Z, Kevenaar ST, Hottenga JJ, Hunter MD, De Zeeuw EL, Neale MC, Van Beijsterveldt CE, Dolan CV, Van Bergen E, Boomsma DI (2021) Multilevel Twin Models: Geographical Region as a Third Level Variable. Behav Genet 51(3):319–330. https://doi.org/10.1007/s10519-021-10047-x
Tubbs JD, Porsch RM, Cherny SS, Sham PC (2020) The Genes We Inherit and Those We Don’t: Maternal Genetic Nurture and Child BMI Trajectories. Behav Genet 50(5):310–319. https://doi.org/10.1007/s10519-020-10008-w
Tucker-Drob EM, Bates TC (2016) Large Cross-National Differences in Gene × Socioeconomic Status Interaction on Intelligence. Psychol Sci 27(2):138–149. https://doi.org/10.1177/0956797615612727
van der Sluis S, Dolan CV, Neale MC, Posthuma D (2008) Power Calculations Using Exact Data Simulation: A Useful Tool for Genetic Study Designs. Behav Genet 38(2):202–211. https://doi.org/10.1007/s10519-007-9184-x
Verhulst B, Maes HH, Neale MC (2017) GW-SEM: A Statistical Package to Conduct Genome-Wide Structural Equation Modeling. Behav Genet 47(3):345–359. https://doi.org/10.1007/s10519-017-9842-6
Yang J, Lee SH, Goddard ME, Visscher PM (2011) GCTA: A Tool for Genome-Wide Complex Trait Analysis. Am J Hum Genet 88(1):76–82. https://doi.org/10.1016/j.ajhg.2010.11.011
Young AI, Nehzati SM, Benonisdottir S, Okbay A, Jayashankar H, Lee C, Cesarini D, Benjamin DJ, Turley P, Kong A (2022) Mendelian Imputation of Parental Genotypes Improves Estimates of Direct Genetic Effects. Nat Genet 54(6):897–905. https://doi.org/10.1038/s41588-022-01085-0

No competing interests reported.

APPENDIXThePowertoResolveCulturalTransmissionandSiblingInteractionUsingPolygenicScores.docx

The Power to Resolve Cultural Transmission and Sibling Interaction Using Polygenic Scores

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Method

Procedure

Materials

Results

Discussion

Conclusion

Declarations

Author Contribution

Acknowledgements

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1