A Note on Jöreskog’s ACDE Twin Model: A Solution, Not the Solution

Abstract

The classical twin design is used to decompose phenotypic variance into genetic and environmental variance components. In practice, for reasons of identification, either an ADE or an ACE twin model is fitted, i.e., models with three variance components. Jöreskog proposed to use the Moore-Penrose inverse to estimate the four genetic and environmental variance components in the ACDE twin model. We demonstrate that the variance components thus obtained do not equal the true genetic and environmental variance components. Given a ACDE model, resorting to an ADE or ACE model or applying Jöreskog’s method will not produce correct results. However, it is possible that Jöreskog’s method will produce a better approximation to the true variance components under certain configurations of A, D, C, and E variance components.

Keywords:

• Identification
• Moore–Penrose inverse
• twin model
• variance components

The univariate classical twin design is a means to the end of decomposing the variance of a given phenotype into the genetic and environmental variance components, based on the phenotypic (2 × 2) monozygotic (MZ) and dizygotic (DZ) covariance matrices, Σ_MZ and Σ_DZ. On the environmental side, shared or common (C) environmental and unshared (E) environmental variances are distinguished: σ_C² and σ_E², where C is defined as all non-genetic effects that increase resemblance between relatives (twins). On the genetic side, additive genetic (A) and dominance genetic (D) variance components are distinguished: σ_A² and σ_D². Subject to many assumptions (Eaves et al., 1978), including the independence of the genetic and environmental variables and random mating, the phenotypic variance equals σ_Ph² = σ_A² + σ_D²+ σ_C² + σ_E². The observed MZ and DZ (2 × 2) covariance matrices contain three statistics that have distinct model expectations: the phenotypic variance and the MZ and DZ twin covariances. These three statistics are insufficient to resolve all four variance components. As a solution, one may choose, as a point of departure, either the ACE model (σ_Ph² = σ_A² + σ_C² + σ_E²), or the ADE model (σ_Ph² = σ_A² + σ_D² + σ_E²). The choice is based on the following rule of thumb: ACE if r_MZ <2*r_DZ, and ADE if r_MZ >2*r_DZ, where r_MZ and r_DZ are phenotypic correlations. Another solution is to estimate an unsigned variance component for C, and to note that a significant and substantial negative estimate of σ_C² implies that variance due to genetic non-additivity outweighs variation due to common environment, assortative mating, or age effects (Verhulst et al., 2019). Structural equation modelling (SEM) is used for parameter estimation and model testing, using programs like LISREL (Cardon et al., 1991; Jöreskog & Sörbom, 2018) or OpenMx (Neale et al., 2016). So, from the point of view of SEM, the full ACDE model is not identified. This can be formalized in SEM terms in terms the rank of the information matrix or the rank of the Jacobian (Bollen & Bauldry, 2010). Equivalently, it can be formalized terms of the set of equations that related the unknown variance components to the observed information (Jöreskog, 2021). In the standardized case (σ_Ph² = 1), we have $${\sigma }_{\text{Ph}}{}^{\text{2}}=\text{1}={\sigma }_{\text{A}}{}^{\text{2}}+{\sigma }_{\text{D}}{}^{\text{2}}+{\sigma }_{\text{C}}{}^{\text{2}}+{\sigma }_{\text{E}}{}^{\text{2}}$$ $${\text{r}}_{\text{MZ}}={\sigma }_{\text{A}}{}^{\text{2}}+{\sigma }_{\text{D}}{}^{\text{2}}+{\sigma }_{\text{C}}{}^{\text{2}}$$ $${\text{r}}_{\text{DZ}}=\mathrm{½}\sigma {}_{\text{A}}{}^{\text{2}}+\mathrm{¼}{\sigma }_{\text{D}}{}^{\text{2}}+{\sigma }_{\text{C}}{}^{\text{2}}$$ where the coefficients ½ and ¼ are based on quantitative genetic theory (Mather & Jinks, 1971). This may be expressed as y = Bx, where the column vector x^t= [σ_A² σ_D² σ_C² σ_E²], and the column vector y^t = [1, r_MZ, r_DZ]. The solution involves the inverse is B⁻¹y = x. In the ACDE model the (3 × 4) matrix B has rows [1,1,1,1], [1,1,1,0] and [½, ¼, 1, 0], and the inverse B⁻¹ does not exist. As noted, setting either σ_C² to zero (r_MZ <2*r_DZ) or σ_D² to zero (r_MZ >2*r_DZ) is the common solution to this problem.

Jöreskog (2021) recently proposed an alternative method, namely by using the Moore–Penrose inverse to obtain an estimate of x: x=(B^t(BB^t)⁻¹)y. The Moore–Penrose provides the minimum norm solution, i.e., it minimizes √(x^tx) given y = Bx. This method can be implemented in programs like LISREL and OpenMx.

The imposition of this mathematical constraint certainly does produce a solution to y = Bx, but we question whether it is the correct quantitative genetic solution. To illustrate this, we consider a phenotype, which is subject to the effect of a single biallelic locus, with two alleles denoted M and N. Let p denote the frequency of allele M (in a well-defined population), and q = 1-p denote the frequency of allele N. Assuming Hardy Weinberg equilibrium, the genotype frequencies are p² (genotype MM), 2pq (genotype MN) and q² (genotype NN). The conditional phenotypic means in the subpopulation of individual with genotype MM, MN and NN are µ+a, µ+d, and µ-a, respectively, and the genetic contributions of the locus to the phenotypic variance are (Mather & Jinks, 1971; Posthuma et al., 2003). $${\sigma }_{\mathrm{A}}{}^2=\mathrm{ }2\text{pq}{\left(\mathrm{a}+\mathrm{d}\left(\mathrm{q}-\mathrm{p}\right)\right)}^2$$ $${\sigma }_{\mathrm{D}}{}^2={\left(2\text{pqd}\right)}^2$$

For example, given p = 0.5, a = 0.5 and d = a (complete dominance), we have σ_A² = 0.125, and σ_D² = 0.0625. Suppose furthermore that σ_C² = 0.09375 and σ_E² = 0.09375, so that σ_Ph² = 0.375. The standardized values are σ_Ast² = 4/12, σ_Cst² = 3/12, σ_Dst² = 2/12, and σ_Est² = 3/12. The MZ and DZ covariance matrices are shown in Table 1.

Table 1. MZ and DZ covariance matrices.

	Twin 1	Twin 2
ΣMZ
Twin 1	0.37500	0.28125
Twin 2	0.28125	0.37500
ΣDZ
Twin 1	0.37500	0.171875
Twin 2	0.171875	0.37500

The correlations are r_MZ = 0.75 and r_DZ = 0.4583, which, by the aforementioned rule of thumb, suggests, incorrectly, an ACE model (r_MZ <2*r_DZ), with estimates, based of Falconer’s equations (Falconer & Mackay, 1996), of 2*(r_MZ-r_DZ) = 0.58333 (σ_Ast²) and 2*r_DZ-r_MZ = 0.1666 (σ_Cst²). Applying Jöreskog’s method, as implemented in OpenMx, we obtained the estimates shown in Table 2.

Table 2. True variance components and estimated variance components based on the Moore–Penrose inverse.

	A	C	D	E
True raw	0.125	0.09375	0.0625	0.09375
True standardized	0.3333	0.250	∼0.1667	0.250
ACDE raw	∼0.0913	∼0.1047	∼0.0846	0.0936
ACDE standardized	∼0.2440	∼0.2798	∼0.2262	0.250

We note that the estimates based on the Moore-Penrose inverse do not equal the true variance components. For instance, the true narrow-sense (σ_A²/σ_Ph²) and broad-sense ((σ_A² +σ_D²)/σ_Ph²) heritabilities are ∼0.3333 (0.125/0.375) and 0.5 ((0.125 + 0.0625)/0.375), respectively. However, the heritabilities based on the results obtained with the Moore-Penrose inverse equal ∼0.2440 (narrow sense) and ∼0.4702 (broad sense).

We conclude that the Moore-Penrose inverse provides a solution for the ACDE model by the imposition of a mathematical constraint. However, it is clear from our results that it does not provide the quantitative genetic solution, i.e., the variance components associated with the genetic and environmental influences. In this connection, we also note that Jöreskog’s Eq. 23, which provides the correct narrow-sense heritabilities in the ACE and ADE model (i.e., σ_A²), does not do so in the ACDE model. Specifically, according to Eq. 23, and given a standardized phenotype (σ_Ph² = 1), the narrow sense heritability equals (1/2)*rMZ - (2/7)*rDZ = (5/14)*σ_Ast²+ (3/14)*σ_Cst² + (6/14)*σ_Dst², whereas it should equal σ_Ast².

While we conclude that the solution provided by the Moore–Penrose inverse not the correct (i.e., quantitative genetic) solution, we quickly acknowledge that the common practice of resorting to an ACE or ADE model will not produce the correct results either, if the true model is the ACDE model. In that situation, both solutions produce incorrect results, which may be expressed in terms of an error of approximation. It is likely that, under certain configurations of A, D, C, and E variance components, Jöreskog’s approach has a smaller error of approximation. In closing, we point out that other genetically informative designs have been proposed, in which the four variance components are identified (e.g., Keller et al., 2009, 2010; Truett et al., 1994).