Introduction

The classical twin design (CTD; Eaves et al. 1978; Jinks and Fulker 1970) has been one of the most productive genetically informative designs in the study of human traits (Polderman et al. 2015). Twin studies have contributed greatly to our knowledge concerning genetic and environmental contributions to individual differences in psychological and medical traits, disease phenotypes and 'omics' variables (van Dongen et al. 2012). Multivariate and longitudinal extensions of the CTD have provided insights into the etiology of comorbidity and stability of traits and disorders. It is well understood that the correct interpretation of results based on the CTD depend on the tenability of the model assumptions (Eaves et al. 1977; Jinks and Fulker 1970; Plomin et al. 2016). The main assumptions of the CTD concern genotype-environment covariance (assumed to be absent), genotype-environment interaction (assumed to be absent), the equal environment assumption (environment does not cause larger resemblance in MZ than in DZ twins), and parental mating (assumed to be random, or that parental resemblance is due to social homogamy rather than phenotypic assortment). Given these assumptions, the results from the CTD can provide unbiased estimates of additive genetic (A), unshared environmental (E), and common environmental (C) and dominance (D) variance components. The effect of violations of these assumptions are well understood (Verhulst and Hatemi 2013; Purcell 2002; Keller et al. 2010), so that estimates of variance components obtained in the twin model may be interpreted in the light of possible model violations.

Many papers have been devoted to the detection and accommodation of model violations, either within the CTD (e.g., Purcell 2002; Molenaar et al. 2012; Eaves and Erkanli 2003; Carey 1986; Dolan et al. 2014; Beam and Turkheimer 2013), or in extended designs (e.g., Plomin et al. 1985; Narusyte et al. 2008; Neale and Fulker 1984; Fulker 1988; D’Onofrio et al. 2003; Keller et al. 2009; Heath et al. 1985; Maes et al. 2006). The aim of the present paper is to demonstrate that the incorporation of polygenic risk scores (PRSs) in the classical twin design allows one to estimate the covariance between A and C (σAC). In the study of childhood intelligence, σAC > 0 is considered plausible, stemming from a process of cultural transmission (Keller et al. 2009; Fulker 1988), which gives rise to passive genotype-environment covariance in children (Plomin et al. 1977; Scarr and McCartney 1983, Knafo and Jaffee 2013; Kendler 2011; Rutter and Silberg 2002).

Measured genetic variables have been incorporated in genetically informative designs with various aims, apart from the gene finding of traditional linkage or combined linkage-association analysis (e.g., Fulker et al. 1999; Neale 2000). For instance, van den Oord and Snieder (2002) presented an extended twin model with measured genetic variables to test association in the presence of population stratification and to test causal relationships. Neale et al. (2000) partitioned variation in serum APOE levels into that associated with the APOE locus and residual genetic variance. In a study of attention problems, van Beijsterveldt et al. (2011) incorporated measured candidate gene information on SNPs in the serotonergic, dopaminergic system and the BDNF gene. The effect of SNPs was tested on a latent factor that summarized multiple assessments of attention problems across childhood. Minică et al. (2018) presented an integration of the CTD and Mendelian randomization method, in which PRSs feature as genetic instruments.

The use of measured genetic information specifically to study genotype-environment covariance is relatively new. Bates et al. (2018) and Kong et al. (2018) proposed the use of polygenic scores based on transmitted and non-transmitted alleles from parents to offspring to detect the effects of non-transmitted alleles on phenotype outcomes in their children. Warrington et al. (2018) used structural equation modeling to determine the fetal and maternal effects of measured genetic variants on birthweight, thus revealing genotype-environment covariance. Cheesman et al. (2020) applied PRSs in an adoption design to detect (passive) gene-environment correlation in educational attainment. Selzam et al. (2019) used PRSs measured in DZ twin pairs to demonstrate the presence of gene-environment correlation for cognitive abilities, and the mediating role of social economic status therein. Wertz et al. (2018) used PRSs in a parent and offspring design to demonstrate the gene-environment correlation originating in parental behavior.

The present aim is to incorporate PRSs in the twin design with the aim of estimating A–C covariance. This approach allows us to determine the presence of A–C covariance, but sheds no light on the process that gave rise to the A–C covariance. For instance, the A–C covariance may be due to active (e.g., niche picking), passive (e.g., cultural transmission) or evocative processes (Plomin et al. 1977). The outline of this paper is as follows. First, we present the classical twin model, and the model extended with PRSs. Second, given the model for PRSs in MZ and DZ twins, we address the issues of identification and statistical power. Third, we present the results of a small simulation to determine the effects of using estimated weights in calculating PRSs (i.e., the standard procedure) in comparison to exact known weights.

The Twin Model with Polygenic Risk Scores: A–C Covariance

Let Ph denote the phenotype of interest, and let GVk denote the k-th genetic variant (GV) contributing to the variance of Ph, where k = 1…K, and K is the number of GVs. We limit our presentation to diallelic GVs (e.g., SNPs) with additive effects (additively coded, e.g., 0, 1, or 2). The phenotype Ph is modeled as follows:

$${\text{Ph}} = {\text{b}}_{0} + \sum\nolimits_{k = 1}^{K} {b_{k} } {\text{GV}}_{{{\text{ki}}}} + {\text{C}}_{{\text{i}}} + {\text{E}}_{{\text{i}}} ,$$

where b0 is the intercept, bk is the k-th regression coefficient, subscript i denotes person, and E and C represent unshared and shared environmental factor scores of individual i. In the classical twin model, under the assumptions mentioned in the introduction, the variance of Ph is decomposed, in the ACE model, into the components σ2A, σ2C, and σ2E:

$${\upsigma }^{2}_{{{\text{Ph}}}} = {\upsigma }^{2}_{{\text{A}}} + {\upsigma }^{2}_{{\text{C}}} + {\upsigma }^{2}_{{\text{E}}} .$$

The additive genetic variance equals the sum of the contributions of the individual GVs and their covariances (σGVk,GVl) attributable to linkage disequilibrium:

$${\upsigma }^{2}_{{\text{A}}} = \sum\nolimits_{{{\text{k}} = 1}}^{{\text{K}}} {{\text{b}}_{{\text{k}}}^{2} {\upsigma }^{2}_{{{\text{GVk}}}} } + \sum\nolimits_{{{\text{k}} = 1}}^{{\text{K}}} {\sum\nolimits_{{{\text{l}} = 1,{\text{l}} \ne {\text{k}}}}^{{\text{K}}} {{\text{b}}_{{\text{k}}} {\text{b}}_{{\text{l}}} } } \,\,{\upsigma }_{{{\text{GVk}},{\text{GVl}}}}$$

Given σAC ≠ 0, we have (assuming no AE covariance)

$${\upsigma }^{2}_{{{\text{Ph}}}} = {\upsigma }^{2}_{{\text{A}}} + {\upsigma }^{2}_{{\text{C}}} + {\upsigma }^{2}_{{\text{E}}} + 2{\upsigma }_{{{\text{AC}}}} ,$$

where σAC is the covariance of A and C. In this model, the parameter σAC is not identified. If we assume σAC = 0, while in truth σAC > 0, the variance σ2C is biased in the twin model, as σAC acts as C, thus inflating the estimate of σ2C in the standard ACE model (see Purcell 2002; Verhulst and Hatemi 2013). If in truth, σAC < 0, the MZ and DZ twin correlations suggest the presence of dominance variance (D).

Given estimates of the regression coefficients (bk) obtained in independent genome-wide association studies (GWASs), the PRS can be calculated \({\sum }_{\mathrm{l}=1}^{\mathrm{L}}{\mathrm{b}}_{\mathrm{l}}\) GVli (Purcell et al. 2009; Evans et al. 2009; Dudbridge 2013), where the set of L GVs is a subset of the K GVs. The set of L GVs may be chosen on the basis of the p value of the individual GVs or other considerations. Let the PRS equal p*Ap, with variance p2*σ2Ap. The scaling parameter p accommodates the fact that the phenotype Ph and the PRS are not measured on the same scale. Let Aq denote the residual additive genetic variable. The model is now:

$${\text{PRS}}_{{\text{i}}} = {\text{p}}\,{\text{A}}_{{{\text{pi}}}}$$
$${\text{P}}_{{{\text{Phi}}}} = {\text{b}}_{0} + {\text{A}}_{{{\text{pi}}}} + {\text{A}}_{{{\text{qi}}}} + {\text{C}}_{{\text{i}}} + {\text{E}}_{{\text{i}}} ,$$

The variance decomposition of the additive genetic variable A and the phenotype Ph are:

$${\upsigma }^{{2}}_{{{\text{PRS}}}} = {\text{p}}^{{2}} {\upsigma }^{{2}}_{{{\text{Ap}}}}$$
$${\upsigma }^{{2}}_{{\text{A}}} = {\upsigma }^{{2}}_{{{\text{Ap}}}} + {\upsigma }^{{2}}_{{{\text{Aq}}}} + 2{\upsigma }_{{{\text{ApAq}}}}$$
$${\upsigma }^{{2}}_{{{\text{Ph}}}} = {\upsigma }^{{2}}_{{{\text{Ap}}}} + {\upsigma }^{{2}}_{{{\text{Aq}}}} + {\upsigma }^{{2}}_{{\text{C}}} + {\upsigma }^{{2}}_{{\text{E}}} + 2{\upsigma }_{{{\text{AC}}}} ,$$

where σAC = σApC + σAqC. The parameters σApC and σAqC are the covariances of Ap and C and of Aq and C, respectively. The parameter σApAq is the covariance of the additive variables Ap and Aq. We parameterize the covariance terms σApC and σAqC as a function of the single covariance term σAC as follows. We derive the coefficient γp by tracing from C to A (C ↔ A with coefficient σAC), and then from A to Ap (A → Ap) where γp is the regression coefficient in the regression of Ap on A. We do the same with Aq using γq. The path diagram is shown in Fig. 1.

Fig. 1
figure 1

The covariance between C and Ap and Aq are derived as σACγp and σACγq, respectively, where γp = σ2Ap/σ2A and γq = σ2Aq/σ2A

Full size image

We thus obtain the constraints:

$${\upsigma }_{{{\text{ApC}}}} = {\upgamma }_{{\text{p}}} {\upsigma }_{{{\text{AC}}}} ,$$
$${\upsigma }_{{{\text{AqC}}}} = {\upgamma }_{{\text{q}}} {\upsigma }_{{{\text{AC}}}} ,$$

where

$${\upgamma }_{{\text{p}}} = {{\left\{ {{\upsigma }^{{2}}_{{{\text{Ap}}}} + {\upsigma }_{{{\text{ApAq}}}} } \right\}} \mathord{\left/ {\vphantom {{\left\{ {{\upsigma }^{{2}}_{{{\text{Ap}}}} + {\upsigma }_{{{\text{ApAq}}}} } \right\}} {{\upsigma }^{{2}}_{{\text{A}}} }}} \right. \kern-\nulldelimiterspace} {{\upsigma }^{{2}}_{{\text{A}}} }},$$
(1)
$${\upgamma }_{{\text{q}}} = {{\left\{ {{\upsigma }^{{2}}_{{{\text{Aq}}}} + {\upsigma }_{{{\text{ApAq}}}} } \right\}} \mathord{\left/ {\vphantom {{\left\{ {{\upsigma }^{{2}}_{{{\text{Ap}}}} + {\upsigma }_{{{\text{ApAq}}}} } \right\}} {{\upsigma }^{{2}}_{{\text{A}}} }}} \right. \kern-\nulldelimiterspace} {{\upsigma }^{{2}}_{{\text{A}}} }}.$$
(2)

Note that σAC = σApC + σAqC, as γ1 + γ2 = 1. At this point two comments are in order. First, it is not possible to estimate both the scaling parameter p and the parameter σApAq. We therefore set σApAq to equal zero. Given this identifying constraint, γp = σ2Ap/σ2A and γq = σ2Aq/σ2A. We demonstrate below that the constraint σApAq = 0 has no bearing on the likelihood ratio test of σAC = 0, or on the maximum likelihood estimates of σAC, σ2A, and σ2C. Second, we recognize that if σAC ≠ 0, the PRS weights (bl) obtained in meta-analyses of the results of independent GWASs, which are used to calculate the PRS, will be upwardly (downwardly) biased given σAC > 0 (σAC < 0). This raises the question of whether this has any effect on the estimate of σAC. We address this question below. The model is depicted in Fig. 2. We have parameterized this model in terms of variance components, i.e., we fixed the path coefficients terminating in the phenotype to one, and estimated the variance components (σ2Ap, σ2Aq, σ2C, σ2E, along with the parameter p). One may also consider fixing the variance components to one, and estimating the path coefficients. However, this parameterization may complicate statistical tests of the variance components (see Verhulst et al. 2019, for details).

Fig. 2
figure 2

ACE twin model with PRSs, including A–C covariances σApC and σAqC (dashed double headed arrows). This is the model in DZ twins (i.e., rz = 0.5). The covariance between Ap and Aq is fixed to zero, but, as demonstrated in the text, this has no bearing on the derived estimate of the total A, C covarianc (σA,C)

Full size image

Simulation I: Power

In this model, the observed statistics are the 3 × 3 MZ (PRS, phenotype twin 1, phenotype twin 2) and the 4 × 4 DZ (PRS1, PRS2, phenotype twin 1, phenotype twin 2) covariance matrices (ΣMZ and ΣDZ, respectively), and the 3- and 4-dimensional mean vectors. The MZ covariance matrix ΣMZ is 3 × 3, as MZ twins, being genetically identical, have identical polygenic scores. We do not consider the mean structure of the phenotype data beyond noting that we adopt the standard (testable) assumptions that the means are equal over twins within a pair and over zygosity. Let the vector θ contain the six parameters of the covariance structure model:

$${\mathbf{\theta = }}\left[ {{\text{p}}\quad {\upsigma }^{2}_{{{\text{Ap}}}} \quad {\upsigma }^{2}_{{{\text{Aq}}}} \quad {\upsigma }^{2}_{{\text{C}}} \quad {\upsigma }^{2}_{{\text{E}}} \quad {\upsigma }^{2}_{{{\text{AC}}}} } \right].$$

The vector θ does not include σApC and σAqC explicitly, because these parameters depend on the parameters σ2Ap, σ2Aq, and the total covariance σAC, as shown above. We evaluated local identification numerically using the OpenMx function mxCheckIdentification (written by Michael Hunter). The model is locally identified if the Jacobian matrix, J(θ), is full column rank (Bekker et al. 1994). The Jacobian matrix contains the first-order derivatives of the non-redundant elements in the matrices ΣMZ(θ) and ΣDZ(θ) with respect to the parameters in θ. Given 6 + 10 elements in ΣMZ(θ) and ΣDZ(θ), and 6 parameters J(θ) is a 16 × 6 matrix. The mxCheckIdentification function is convenient as it does the necessary calculations automatically, and can be applied directly to the OpenMx script that one uses to fit the model.

Having established local identification, we proceeded to address the question of resolution by considering the statistical power to reject σAC = 0 given various parameter settings. We used exact data simulation (van der Sluis, et al. 2008), which is equivalent to analyzing the expected (true) covariance matrices. We used normal theory maximum likelihood estimation throughout, and based our power calculations on the non-centrality parameter (NCP) associated with the (non-central) chi-square distribution (Martin et al. 1978). Given exact data simulation, the NCP equals the loglikelihood ratio (LLR) test of σAC = 0. We used the OpenMx library (Boker et al. 2011; Neale et al. 2016) in the R program (R Core Team 2018), and we used R for data simulation, and power calculations. In the power analyses, we set the MZ and DZ sample sizes equal (Nmz = 1000, Ndz = 1000), and we report the power given Nmz = Ndz = 1000, and the required sample sizes to achieve a power of 0.80, given α = 0.05.

Results I

The numerical check of model identification demonstrated that the model is locally identified, bearing in mind that we have set σApAq = 0. That is, given the 3 × 3 MZ and the 4 × 4 DZ phenotypic covariance matrices, ΣMZ(θ) and ΣDZ(θ), we can obtain unique estimates of the six parameters p, σ2Ap, σ2Aq, σ2C, σ2E, and σAC. From the perspective of the path model (Fig. 2), the key to the identification is the covariance between the phenotype and the PRS, which does not depend on zygosity. This covariance equals p*σAp2 + p*σApC or p*γ1*(σA2 + σAC), where p, γ1 and σA2 are identified based on the phenotypic and PRS MZ and DZ twin covariances. Table 1 contains the results of the power study. This table includes the 16 parameter settings and the power to reject σAC = 0.

Table 1 Statistical power to reject the null hypothesis that A–C covariance is zero (alpha = 0.05)
Full size table

Table 1 contains the proportions of genetic and phenotypic variance explained by the PRS (prPRS and prPh in Table 1). These range from 0.2 to 0.4 (prPRS) and 0.050 to 0.177 (prPh). The correlation between A and C (rAC) was chosen to equal 0.2 or 0.3. In addition to this correlation, we express the σAC effect size as the proportion (2σAC)/σ2Ph (i.e., pr2σAC in Table 1), where σ2Ph = σ2A + σ2C + σ2E + 2σAC. This proportion ranges from 0.088 to 0.189. The proportion of A variance is ~ 0.26 or ~ 0.40; the proportion of C variance is ~ 0.18 or ~ 0.25, and the proportion of E variance varies between 0.09 and 0.44). The settings, which are limited, were chosen merely to identify some circumstance in which the power to reject σAC = 0 is acceptable, given the present sample sizes.

The greatest power is obtained in settings 8, 14, and 16: 0.906 (8), 0.793 (14), and 0.950 (16). Here, the PRS accounts for 10.1% (8), 16.8% (14), and 16.2% (16) of the phenotypic variance, and 2σAC accounts for 15.3% (8), 15.9% (14), and 18.9% (16) of the phenotypic variance. We see the lowest power given settings 1 (0.216) and 9 (0.238). Unsurprisingly, these are associated with relative low values of pr2σAC (8.8% and 11.2%) and prPh (5.4% and 8.8%). The relative contributions of pr2σAC and prPh to the power are apparent in the correlations of these with the power: 0.71 and 0.53, respectively. Regressing the power on these parameters, we found that they explain 65% of the variance in power (β pr2σAC: 0.625 and β prPh 0.400). Both are important, but the contribution of pr2σAc to the power is greater. The comparison of the σ2A = 0.3 and the σ2A = 0.5 conditions show that the magnitude of σ2A does not greatly influence the power. In terms the ratio of the N (power = 0.80) are about 1.1–1.2 (favoring the σ2A = 0.5 conditions). Finally, to determine the effect of the sign of σAC, we repeated the power analysis with rAC set to equal − 0.2 or − 0.3, and all other parameters unchanged. Figure 3 displays the plot of the power given positive and negative rAC. We note that the difference in power is small suggesting that the sign of rAC is unimportant in calculating power.

Fig. 3
figure 3

Power of the LLR test to reject σAC = 0 given positive and negative σAC. The parameter settings are given in Table 1. The only difference is the sign of σAC. The power to reject σAC = 0 given positive σAC is given in Table 1

Full size image

Simulation II: Bias and Type I Error Rate

As noted above, the weights (bl) used to calculate the PRS are expected be biased upwards if σAC > 0. We investigated the effects of this in a simulation study. Specifically, we simulated data according to the nuclear twin family (NTF) design (Fulker 1988; Keller et al. 2009), which comprises MZ and DZ twins and their parents. This model includes cultural transmission, i.e., the direct contribution of the parental phenotype to the twin's environment (parameter m in the notation of Keller et al. 2009 notation). This contribution gives rise to a shared environmental variable (F in the notation of Keller et al. 2009), and to covariance between this variable (F) and the additive genetic factor (A), σAF. In addition to the shared variance due to F, the model includes a shared environmental variance term, due to shared influences other than cultural transmission. We denote this C* to distinguish it from the C in the standard ACE model. The original NTF model accommodates phenotypic assortative mating. However, here we assume that mating is random. The MZ and DZ expected phenotypic covariance matrices are shown in Table 2. It is not possible to resolve F and C* in the absence of parental phenotypes. Thus, in the standard ACE model and in the model with PRSs, we estimate a single shared environmental variance, σ2C which actually equals σ2C* + σ2F.

Table 2 The expected covariance matrices in simulations 1–3
Full size table

First, we considered the model without cultural transmission, i.e., cultural transmission parameter (m) was zero. This implies that σ2F and σAF are zero, as the parameter m is the source of the A–C covariance. However, we included σ2C* > 0. Thus, there is shared (by the twins) environmental variance, but it is not due to cultural transmission. This model allows us to determine the Type I error rate associated with the test of σAC = 0. As an aside, in this simulation, the inclusion of σ2C* > 0 also allows us to establish the power to detect C variance in a twin model with PRS, in addition to checking the Type I error rate in test of σAC = 0. Second, we considered a model with cultural transmission, with m > 0, so that σAF > 0 and σ2F > 0. In this model we set σ2C* to zero, so that there are no shared environmental influences other than those stemming from the cultural transmission. Third, we considered a model with cultural transmission (m > 0), so that σAF > 0 and σ2F > 0, and σ2C* > 0. As mentioned, we cannot resolve σ2F and σ2C*, so we fitted a single shared environmental factor, representing C and F. In summary, given σ2C = σ2F + σ2C*, we have the following settings. Simulation 1: m = 0, σ2F = 0, σ2C* > 0, σ2C = σ2C*, σAC = σAF = 0. Simulation 2: m > 0, σ2F > 0, σ2C* = 0, σ2C = σ2F, σAC = σAF and σAF > 0. Simulation 3: m > 0, σ2F > 0, σ2C* > 0, σ2C = σ2F + σ2C*, σAC = σAF and σAF > 0.

In each simulation study based on these three models, we carried out 500 replications. Each data set comprised genotypic and phenotypic data in parents and twins. The parental data were discarded, and the twin data were used to fit the model. The additive genetic variable comprised 100 uncorrelated diallelic genetic variants, of which 40 were used to calculate the PRS. We carried out the simulation twice: once with exact, unbiased PRS weights (i.e., the parameters b in the expression PRSi = \({\sum }_{\mathrm{l}=1}^{\mathrm{L}}{\mathrm{b}}_{\mathrm{l}}\) GVli), and once with estimated weights b. We estimated the weights in independent data (not the data used to fit the actual model) by regressing the phenotype on the genetic variants. Given the absence of A–C covariance, the estimated weights are unbiased. In each replication the sample sizes were N = 2000 to estimate the parameter weights b, and Nmz = 1000 and Ndz = 1000 (in total 2000 pairs) to fit the actual model. Parameter values and effect sizes are given in Table 3. In simulations 1–3, the true values of prPh are 0.2, 0.141, and 0.147, and the true values of pr2σA,C*+F are 0.0, 0.353, and 0.368, respectively.

Table 3 Means and standard deviation of parameter estimates in simulation 1–3 based on 500 replications (Nmz = 1000; Ndz = 1000)
Full size table

In summary, the aims were (1) to establish that the Type 1 error rate was correct and (2) to investigate the effects, if any, of biased weights on the Type I error rate and the parameter estimates (bias), (3) to determine whether the presence of PRS, given zero A–C covariance, increases the power to detect C variance.

Results II

We first discuss the parameter estimates and then the loglikelihood ratio (LLR) test statistics. Table 3 contains the true parameter values and the mean and standard deviation of the parameter estimates based on the 500 replications. In simulation 1, given exact PRS weights, the parameter estimates are unbiased, as expected. Given estimated PRS weights, the estimates of σ2Ap and σ2Aq are biased: the mean values are 0.184 (underestimated; true 0.2) and 0.316 (overestimated; true: 0.3), respectively. The estimate of the covariance term σAC is unbiased: the mean value is 0.001 (true: 0.0). In simulation two, given exact weights, the parameter estimates are again unbiased, and given estimated PRS weights, the estimates of σ2Ap and σ2Aq are again biased: the mean values are 0.189 (underestimated; true 0.2) and 0.315 (overestimated; true: 0.3), respectively. The estimate of σAC is unbiased: mean value 0.124 (true 0.125). Simulation three produced the same results as simulation two, in terms of parameter bias stemming from using estimated PRS weight. The main finding is that using estimated PRS weights results in biased estimates of the variance components σ2Ap and σ2Aq, but has little effect on the estimate of the covariance term σAC.

Table 4 contains the results of the LLR tests. We tested the hypotheses σA,C = 0, and σ2C = 0 given σA,F+C = 0 in the twin model with the PRSs. In addition, we tested the hypothesis σC2 = 0 in the standard univariate ACE model. The test of σ2C = 0 is of interest in simulation 1, as the comparison of σ2C = 0 given σA,F+C = 0 (in the full model) and hypothesis σ2C = 0 in the standard univariate ACE model tell us whether the presence of PRSs helps to resolve C. The LLR statistic associated with the test of σA,F+C = 0 in simulation 1, where in truth σA,F+C = 0, should follow central chi2(1) distribution, which is characterized by a mean of 1 and a standard deviation of √(2) = 1.414. Given exact PRS weights, the mean and standard deviation of the LLR statistic are 0.970 and 1.423 (see Table 4). These do not differ from the expected values of 1 and √2 (LLR test: 0.27, df = 2, p = 0.87). The Type I error rate equaled 0.049 (CI95: 0.032–0.072). Given estimated weights, the values are 1.041 and 1.515. While these do not appear to deviate from the expected values (LLR test: 5.18, df = 2, p = 0.075), the variance is larger (1.515 vs. 1.414), as is the Type I error rate: 0.063 (CI95: 0.043–0.088). In terms of the mean LLR test, we note that the test of σC2 = 0 is more powerful in the full model (given σA,C = 0) than in the standard univariate ACE twin model. With exact PRS weights, the mean LLR test statistics equal 22.2 (with PRSs) and 16.06 (standard ACE model), and with estimated PRS weight, 21.5 and 15.6.

Table 4 Means and standard deviation of loglikelihood ratio tests simulation 1–3 based on 500 replications (Nmz = 1000; Ndz = 1000)
Full size table

The results of simulation 2 and 3 are comparable. The test of σAC = 0 suffers slightly given estimated PRS weights: the mean LLR test statistics equal 15.10 and 13.5 (simulation 2) and 13.48 and 11.96 (simulation 3). The effect of the test of σ2C = 0 given σAC = 0 is slight (in simulation 3: 78.8 vs 74.6). We note that the test of C in the standard twin model appears to be more powerful. However, given σAC > 0, this is a combined test of σC2 = 0 and σAC = 0.

In simulation 2, we set σ2C* = 0 and σ2F = 0.091 (6.4% of the phenotypic variance). This relatively low value does not rule out a considerable contribution of σAC to the phenotypic variance (35%). The true twin correlations in the simulation 2 (rMZ = 0.74, rDZ = 0.52) suggest a considerable contribution of C (~ 30%). This implies that (1) substantial C in the classical twin model may be mainly due to A–C covariance, while (2) the C variance, which in simulation 2 comprised only σ2F, may be small. Given that this variance may be small, its estimate may, given the present variance component parameterization, assume negative values. Indeed, in simulation 2, we encountered a negative variance component estimate in about 11% (exact PRS weights) and 13% (estimated weights) of the replications (we did not remove these cases in calculating the results in Tables 3 and 4 to avoid the bias caused by truncating the distribution of the parameters to the admissible solutions). The problem of negative variance can be solved by imposing the constraint that the 2 × 2 covariance matrix of Ap and C and the 2 × 2 covariance matrix of Aq and C be positive definite. This means that the eigenvalues of the covariance matrices are constrained to be larger than zero. This PD constraint is simple to implement in OpenMx using an mxAlgebra statement, as OpenMx includes a function to calculate eigenvalues (see the OpenMx script). Only constraining σ2C to be greater than zero is insufficient, as this by itself does not ensure that the covariance matrix of Ap (or Aq) and C is positive definite. In addition, if σ2C were to hit the lower bound, the parameter σAC would no longer be defined. We repeated simulation 2 with these PD constraints to gauge the effects on the parameters. The results, as obtained with estimated PRS weights, are included in Tables 3 and 4. These results are largely consistent with those obtained without the PD constraints. We do note that the PD constraints affect the distribution of the estimates of σ2C positively skewed, as shown in Fig. 4. Without the positive definiteness constraints, this distribution is normal. In contrast, the distribution of the estimates of σAC appear to be quite normal, which suggests that, at least in the present scenario, the LLR test of σAC = 0 will not be affected by the PD constraints.

Fig. 4
figure 4

Distribution of estimates of σ2F (left) and σAF (right) given positive definiteness constraints (simulation 2). The true values are σ2F = 0.091 and σAF = 0.125

Full size image

The Identifying Constraint σApAq = 0

As mentioned above, we have constrained σApAq to equal 0, because it is not possible to estimate the scaling parameter p (in σ2PRS = p2σ2Ap) and the covariance simultaneously. This raises the question how the results are affected if in fact σApAq > 0, as σApAq = 0 is implausible given linkage disequilibrium. In fact, the constraint σApAq = 0 has no effect of the estimates of σAC, σ2A, and σ2C. To demonstrate this, we used simulation exact data with σApAq > 0, and fitted the model twice. Once with σApAq fixed to its true value, and once σApAq fixed to zero. Specifically, we chose the parameter values shown in Table 5. Fitting the model with σApAq fixed to equal its true value (σApAq = 0.1224; correlation: ρApAq = 0.5), we recovered the parameter values, including σAC = 0.077, and the total additive genetic variance 0.2 + 0.3 + 2*0.1224 =  ~ 0.745 (i.e., σ2Ap + σ2Aq + 2σApAq). The power to reject σAC = 0 equals 0.742 (α = 0.05). Fixing σApAq = 0, we obtain identical results, except for the values of σ2Ap and σ2Aq. The total genetic variance is now composed as follows 0.5199 + 0.2250 = 0.745 (i.e., σ2Ap + σ2Aq). We checked this result with a wide variety of parameter values. So in principle, one can constrain σApAq to equal any sensible value. However, while the estimates and tests of σAC, σ2A, and σ2C are unaffected, we note that the values of σ2Ap and σ2Aq do depend on this sensible value.

Table 5 Results with σApAq fixed to its true value (row A), and σApAq fixed to equal zero (row B)
Full size table

Discussion

The present aim was to estimate A–C covariance in the classical twin model with PRS. To this end, we proposed the model depicted in Fig. 1, in which the covariance between Ap (PRS) and Aq and C are modeled as a function of the single covariance of A (Ap + Aq) and C. We found that the power to reject σAC = 0 depends mainly on the proportion of phenotypic variance due to the covariance term (σAC) and the PRS, where the former is more important than the latter. We investigated the influence of using estimated PRS weights. The use of estimated weights resulted in downwards bias of σAp and an upwards bias of σAq. However, the estimate of σAC was not affected. The use of estimated weights had a small effect of the Type I error rate in the test of σAC = 0.

In the most favorable settings qua power (8, 14, and 16 in Table 1), the proportions pr2σAC (phenotypic variance due to σAC) equaled 0.153, 0.159, and 0.189, and the proportions prPh (phenotype variance due to the PRS) equaled 0.101, 0.168 and 0.162. We consider these values (prPh) to be generally large by today's standards, but note that, while pr2σAC is given, the proportion prPh is likely to increase with the ongoing progress of GWAS meta analyses of many phenotypes. For instance, at present PRSs explain ~ 15% of the variance of educational attainment and ~ 11% of the variance of IQ (Allegrini et al. 2019).

The results of the power study and simulations shed some light on the viability of the model. But the results of simulation 2 also demonstrated that positive cultural transmission can result in large C in the standard ACE model, while most of this C variance is due to σAC. The actual C variance (without σAC) can be quite small. It is therefore advisable to fit the model with the positive definiteness constraint, outlined above. In this connection, we note that the finding that C variance in cognitive abilities is large in young children, but declines quickly in magnitude as children grow older (Haworth et al. 2010; Tucker-Drob and Bates 2016) may well be due to a decline in the magnitude of cultural transmission, in combination with an increase in genetic variance. This may be testing by extending the present model to include age as a moderator in the manner of Purcell (2002). This is relatively simple to do in OpenMx. In this connection, we also note that the estimate of σAC obtained using the present model may tell us that σAC is present. It does not, unlike other models, reveal the source of the σAC. For instance, in the NTF design, cultural transmission is the source the covariance between A and F, where the distinction is made between F (shared environmental effects due to cultural transmission) and residual C, which we denoted C* above.

We considered negative σAC in the power study, and found that the power to reject σAC = 0 was about the same regardless of the sign of σAC. We note that negative σAC (e.g., originating in negative cultural transmission) tends to produce twin correlations, which are suggestive of an ADE model (2*rDZ < rMZ). This is to be expected as − σAC lowers the MZ and DZ correlations to the same extent. Finally, the present results demonstrated that the addition of PRSs to the ACE model increases the power to detect C variance, assuming σAC = 0. This may be of interest, as in the classical twin model, the power to detect C variance is known to be generally poor (Visscher et al. 2008; Martin et al. 1978).

In closing, we note the following limitations. We have assumed that dominance variance (D) is absent, and acknowledge that the twin univariate design is limited to ACE or ADE. As demonstrated by Keller et al. (2010), a well fitting ACE model does not rule out the represent of D. It is possible that the addition of PRSs may aid in resolving D (in an ACDE) model, but we consider this beyond the present scope. Boomsma et al (2020; see this issue) showed that it is possible to estimate all four variance components (A, C, D, and E) in special cases of the multivariate twin model. Second as mentioned above, the settings of the power study and the simulations are limited in scope. In addition, we considered only equal MZ and DZ sample sizes (Nmz = Ndz = 1000). The ratio Nmz/Ndz has a general bearing on the power in the twin design (Visscher 2004). However, power calculations with unequal MZ and DZ sample sizes are simple to carry out. Third, the simulations that we carried out to gauge the effect of using estimated PRS weights, involved only a small number of associated GVs with relatively large effects. Simulation studies with more realistic designs will provide additional information concerning the effects of using estimated PRS weights. Fourth, we assumed that phenotypic mating is random. However, we note that the PRSs in the DZ twins offer the means to tests this, as the correlation between the additive genetic PRSs will equal 0.5 given phenotypic random mating. In addition, the present model may be extended to include parental data to accommodate phenotypic assortative mating as outlined in Keller et al. (2009). Fifth, we have not considered the effect of violations of other standard twin design assumptions on the estimate of σAC (Eaves et al. 1977; Purcell 2002; Keller et al 2010). Genotype—unshared environment covariance (σAE) is not identified in the present model. Unmodeled (positive) σAE and A × C interaction contribute to A. We do not see how either could result in spurious σAC. A × E and C × E interaction contribute to E, which has no bearing on A, C, or σAC. Sixth, we have not considered the possibility that the σAC is due to stratification. We know that spatial (geographical) allele frequency gradients may given rise to spurious C variance in the classical twin design (Tamimy et al. 2020; see this issue). A positive spatial correlation between C effects and allele frequencies, may given rise to A–C covariance. One way to detect this by including as fixed covariates principal components that reflect the allele frequency gradient (Price et al. 2006). If this kind of stratification is an issue, then the size of the C variance (see Tamimy et al. 2020) and the size of the A–C covariance should decline following the introduction of these covariates. Finally, we have assumed that the PRSs weights were obtained from GWASs of the phenotype of interest. Whether the present approach can be adapted to handle PRSs weights based on a genetically correlated phenotype (i.e., correlated with the phenotype of interest) remains to be seen.