1

1

Causality, \(G \times E\) interaction, Missing heritability, Partial identification, Twin studies

1 Introduction↩︎

Nature versus nurture is a long-standing debate in biological and social sciences. Central to this discussion is the concept of heritability, which, roughly speaking, measures the proportion of phenotypic variation due to genetic variation. There are several related yet distinct notions of heritability [1]–[5]. The most commonly used notions are the broad-sense heritability \(H^2\), which refers to the proportion of phenotypic variance that is attributable to total genetic variation; the narrow-sense heritability \(h^2\), which refers to the proportion attributable to additive genetic effects in a statistical genetics model; and the SNP-based heritability \(h_{\text{SNP}}^2\), which refers to narrow-sense heritability based on single-nucleotide polymorphisms (SNP) only. In twin studies, heritability is typically measured by Falconer’s formula (twice the difference of the correlation between identical twins and the correlation between non-identical twins). Closely related are heritability estimates from sibling regression using the identity-by-descent (IBD) kinship matrix [6]. Heritability is also used by plant and animal breeders to measure the expected response to selection of a trait. There, heritability is usually defined on an entry-mean basis (averaging over experimental replications) [7], [8]. Although existing notions of genetic heritability can shed light on the mechanisms that influence human physiology and behavior, they have two serious drawbacks. First and foremost, it has been observed since the beginning of the GWAS era that, for many human traits, there is a large gap between heritability estimates from twin studies and from unrelated individuals [9], [10]. Many theories and explanations for this “missing heritability problem” have been proposed [11]–[13], and this debate has also motivated a flurry of new methods for heritability estimation [10], [14]–[16]. One contributing factor to the “missing heritability problem” is that these heritability estimates use different methods and may simply approximate different population quantities. For example, the twin-based estimate from the ACE model (see Example 2 below) is usually believed to approximate the narrow-sense heritability \(h^2\), and the GWAS-based estimates are thought to approximate the SNP-based heritability \(h_{\text{SNP}}^2\) (this is basically how \(h_{\text{SNP}}^2\) is “defined”). However, it is not easy to communicate the subtle differences between different concepts of heritability to applied scientists and to the public. Fundamentally, the difficulty is that these notions of heritability are based on specific statistical models that make different assumptions for different designs. This means that they are largely ill-defined when the modeling assumptions are not satisfied. This leads us to the second drawback of existing approaches: although heritability is usually described as the fraction of phenotypic variation attributable to genetic variation and thus is an inherently causal concept, most existing notions of heritability only model statistical relations (e.g., associations and conditional associations) instead of causal relations. A notable exception is the ACE model for twin studies [17], [18]. This is a special case of linear structural equation models that were originally introduced by [19] and are routinely used in behavioural genetics [20] to estimate heritability. However, the ACE model still suffers from a lack of causal consideration: the genetic and environmental effects, which are assumed to be independent in the ACE model, are generally confounded due to parental genetic effects. This observation has recently regenerated interests in using family-based designs in human genetics, but mainly for giving causal interpretations to GWAS results [21]–[24] or reducing the bias of Mendelian randomization studies [25], [26]. In sum, although heritability is often thought as a form of causal attribution, the causal perspective is largely missing in the mathematical definitions and models used in the literature.

In this article, we fill this gap with a new, counterfactual definition of heritability. The key idea is to measure heritability by the average magnitude of the squared difference between an individual and its counterfactual—a hypothetical “non-identical twin” that goes through exactly the same environment. The precise construction of this counterfactual will depend on the specific context and more details can be found below. We rigorously define this natural conception of heritability using the Neyman-Rubin potential outcomes model [27] and the nonparametric structural equation model [28], [29] that generalizes Wright’s linear structural equation model. Because our definition of heritability is inherently counterfactual—belonging to the third layer of Pearl’s ladder of causality [30], [31], it is generally only partially identifiable from empirical data. To this end, we provide lower and upper bounds of counterfactual heritability that can be computed from observational data by estimating the conditional distribution of the phenotype given genetic and environmental factors. We further compare counterfactual heritability and its bounds with several existing notions of heritability.

2 Counterfactual Heritability↩︎

2.1 Definition of Counterfactual Heritability↩︎

Let the random variables \(G \in \mathcal{G} \subseteq \mathbb{R}^d\) and \(Y \in \mathcal{Y} \subseteq \mathbb{R}\) denote, respectively, the genotype and phenotype of an individual. Suppose one can construct, at least hypothetically in a thought experiment, the potential outcomes \(Y(g)\) in an experiment that sets \(G\) to a fixed value \(g \in \mathcal{G}\). For example, if \(Y\) represents the weight of a potato, then \(Y(g)\) is the weight of the potato when the genotype is set as \(g\) by CRISPR gene editing.

Our definition of counterfactual heritability rests on comparing the potential outcomes of two “non-identical twins” under exactly the same environment. More specifically, consider the following setting: \[G \perp \!\!\! \perp G'\perp \!\!\! \perp\{Y(g): g \in \mathcal{G}\}\quad and\quad G'\overset{d}=G.\label{orth-uncond}\tag{1}\] Here we assume \(G\) and \(G'\) are two independent and identically distributed (IID) genotypes that are independent of the potential outcomes. In this case, we define counterfactual heritability as \[\label{eq:def-h2} \xi:=\frac{\mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G'))}{2\mathop{\mathrm{\mathrm{Var}}}(Y(G))},\tag{2}\] where \(\xi\) means “explanability”. That is, counterfactual heritability is the average magnitude of the squared difference between the potential outcomes of two “non-identical twins” (with genotype \(G\) and \(G'\)), normalized by dividing twice the phenotypic variance.

To make this definition more concrete, suppose the potential outcomes are determined by the structural equation \[Y(g) = f(g, E),~g \in \mathcal{G},\] where \(f\) is some unknown function and \(E\) contains all environmental factors (possibly unmeasured and including random variations ariging in the organism itself). Condition [1 ] essentially assumes \(G\), \(G'\) and \(E\) are independent, and the counterfactual heritability in [2 ] can be written as \[\xi = \frac{\mathop{\mathrm{\mathrm{Var}}}(f(G,E) - f(G',E))}{2 \mathop{\mathrm{\mathrm{Var}}}(f(G,E))}.\] It is clear from this expression that the “non-identical twins” are assumed to experience exactly the same environment \(E\). As an even more concrete example, if \[Y(g) = \beta_g g + \beta_e E + \beta_{ge} g E,\] it is not difficult to verify that \[\xi = \frac{\beta_g^2 \mathop{\mathrm{\mathrm{Var}}}(G) + \beta_{ge}^2 \mathop{\mathrm{\mathrm{Var}}}(G) \mathop{\mathrm{\mathrm{Var}}}(E)}{\beta_g^2 \mathop{\mathrm{\mathrm{Var}}}(G) + \beta_e^2 \mathop{\mathrm{\mathrm{Var}}}(E) + \beta_{ge}^2 \mathop{\mathrm{\mathrm{Var}}}(G) \mathop{\mathrm{\mathrm{Var}}}(E)}.\] Thus, the counterfactual heritability \(\xi\) includes the variation from gene-environment interaction measured by \(\beta_{ge}^2 \mathop{\mathrm{\mathrm{Var}}}(G) \mathop{\mathrm{\mathrm{Var}}}(E)\). This way of “accounting” is different from most existing notions of heritability that only include the variation of the genetic “main effects” \(\beta_g^2 \mathop{\mathrm{\mathrm{Var}}}(G)\) in this simple model.

It is often not realistic to assume that the gene and environment are independent factors unconditionally. Instead, there may be some environmental factors and/or background genetic factors \(X\) such that, given \(X\), the gene is unconfounded with the rest of the environment. Formally, this unconfoundedness assumption can be described using the following modification to condition [1 ]: \[G \perp \!\!\! \perp G'\perp \!\!\! \perp\{Y(g): g \in \mathcal{G}\} \mid X \quad and\quad G'\overset{d}=G \mid X\label{orth}\tag{3}\] The first part of [3 ] on the conditional independence of \(G\) and potential outcomes of \(Y\) is known as ignorability [27] or no unmeasured confounding [32] in the causal inference literature. Assuming a nonparametric structural equation model with independent errors, this can be read off from a causal diagram for \(G\), \(Y\), and background factors [28]. We will still define the counterfactual heritability as \(\xi\) in equation [2 ], although the distribution of \(G'\) is now different and the variance is also averaged over the distribution of \(X\). One can view [1 ] as a special case of [3 ] that takes \(X\) to be empty. All our theoretical results below will be stated in the more general case under condition [3 ].

Under condition [1 ] or [3 ], \(\xi\) is always bounded between \(0\) and \(1\). In fact, the counterfactual heritability in [2 ] can be written as \[\xi = 1 - \mathop{\mathrm{\mathrm{Cor}}}(Y(G), Y(G')),\] and it can be shown that \(\mathop{\mathrm{\mathrm{Cor}}}(Y(G), Y(G')) \geq 0\) (see Lemma 3). This justifies treating \(\xi\) as a notion of heritability or variable importance.

Note that the definition of counterfactual heritability depends on the distribution of \(G\). We view this as a feature of this definition because heritability is a population genetics concept and thus should depend on the population under investigation. We will only investigate mathematical properties of \(\xi\) in this Section but will return to its practical interpretation later.

2.2 Identifiable bounds on counterfactual heritability↩︎

As mentioned in the Introduction, our definition of counterfactual heritability relies on a hypothetical thought experiment, because no experiment in the real world can provide exactly the same environment for two individuals. Thus, \(\xi\) should be regarded as an idealized quantity that we should strive to approximate as accurately as possible.

From a statistical perspective, the difficulty of approximating \(\xi\) is that it is not possible to observe \(Y(G)\) and \(Y(G')\) at the same time. This is commonly referred to as the “fundamental problem of causal inference” [33]. In fact, in reality we can observe just one realization of the genotype \(G\) and the phenotype \(Y\) for any given individual. This gives rise the so-called causal identification problem: can we use the distribution of the observed variables (in our case, \(G\), \(Y\), and possibly \(X\)) to determine a counterfactual quantity (in our case, the counterfactual heritability \(\xi\))?

To make any progress, we first assume that the potential and factual outcomes are consistent in the sense that \(Y(G) = Y\) [32]; this is also referred to as the stable unit treatment value assumption (SUTVA) [34] or sometimes the no interference assumption in the causal inference literature. By combining this with [3 ], it is possible to identify the distribution of \(Y(g)\) or the average treatment effect \(\relax(Y(g) - Y(g'))\) for any fixed \(g, g'\) [35].

However, the counterfactual heritability \(\xi\) in [2 ] clearly depends on the joint distribution of the potential outcomes, which cannot be identified from empirical data. So exact identification of \(\xi\) is not possible in general. Nonetheless, it is possible to provide nontrivial bounds on \(\xi\) as we demonstrate below.

Bounding \(\xi\) is an instance of optimal transport problem. Note that bounding \(\xi\) is equivalent to bounding \(V = \mathop{\mathrm{\mathrm{Var}}}(Y(G) - Y(G'))\), because the denominator in \(\xi\) can be directly identified by \(2\mathop{\mathrm{\mathrm{Var}}}(Y(G)) = 2 \mathop{\mathrm{\mathrm{Var}}}(Y)\). When the genotype \(G\) is binary, the next result provides tight upper and lower bounds of \(V\) using the Fréchet-Hoeffding inequality [36], [37].

When \(G\) is not binary, \(\xi_l\) is still a tight lower bound for \(\xi\), but obtaining tight upper bounds for \(V\) (and thus \(\xi\)) is still an open problem [38]. To provide some intuition, it can be shown (see Lemma 3 in the Appendix) that if \(X\) is empty (so [3 ] reduces to [1 ]), we have \[\label{1-h2} 1 - \xi = \frac{\mathop{\mathrm{\mathrm{Var}}}\{ \relax(Y(G) \mid Y(\cdot))\}}{\mathop{\mathrm{\mathrm{Var}}}\{Y(G) \}},\tag{4}\] where \(Y(\cdot) = \{Y(g): g \in \mathcal{G}\}\) denotes the potential outcomes “schedule”. Thus, if there exists a joint distribution on \(Y(\cdot)\) (with constrained marginals given by the observed data) such that \(\relax(Y(G)\mid Y(\cdot))\) is a constant, the trivial upper bound \(\xi \leq 1\) is tight; see [39] for some examples. Heuristically, if we do not have any environmental factor and observe different phenotypes for two “non-identical twins”, we cannot rule out the possibility that such difference is entirely due to the difference in their genotypes.

The next result provides bounds in the non-binary case. All mathematical proofs can be found in the Appendix.

Theorem 1 provides two lower bounds and tow upper bounds on the counterfactual heritability \(\xi\) that can be computed using a sample of \((X,G,Y)\). The lower bound \(\xi_l\) is tight but requires modeling the conditional quantile function of \(Y\) given \(X\) and \(G\). The more relaxed lower bound \(\xi_{l}'\) could be easier to compute. On the other side, the new upper bound \(\xi_u'\) uses the expected conditional variance of \(Y\) given the observed environmental factors. Thus, including some environmental factors \(X\) that are independent of the genes but predictive of the phenotype can lead to non-trivial upper bounds of the counterfactual heritability. Intuitively, the genotypes can explain at most the variance not possibly explained by \(X\), so \[\frac{1}{2} \mathop{\mathrm{\mathrm{Var}}}(Y(G) - Y(G')) \leq \relax(\mathop{\mathrm{\mathrm{Var}}}(Y \mid X)).\] The inequality \(\xi \leq \xi_u\) in Proposition 1 still holds for non-binary \(G\), but it is generally not very useful. For example, when \(\mathbb{P}(G = g) = 1/3\) and \(Y(g) \sim \text{N}(0,1)\) for \(g \in \{0,1,2\}\), it is not difficult to show that \(\xi_u = 4/3 > 1\) and thus the upper bound \(\xi_u\) is trivial.

3 Connection with different notions of heritability↩︎

3.1 Comparison with broad-sense and narrow-sense heritability↩︎

Next, we compare counterfactual heritability \(\xi\) with other existing notions of heritability. Roughly speaking, “broad-sense heritability”, often denoted as \(H^2\), is the proportion of phenotypic variation statistically explained by total genetic variation, including dominance and epistasis. On the other hand, “narrow-sense heritability”, often denoted as \(h^2\), refers to phenotypic variation explained by additive genetic variation only. These quantities are not always precisely defined in the literature (their definition usually requires a specific statistical model). Consider the example below.

Next, we compare \(H^2\) and \(h^2\) with counterfactual heritability \(\xi\) and its bounds in Theorem 1. This is not straightforward because the precise definitions of \(H^2\) and \(h^2\) are often tied to specific parametric models of the phenotype. To make these quantities comparable to \(\xi\) when the parametric models are not correct, we will use the following nonparametric definition \[\label{def:H2} H^2 = \frac{\mathop{\mathrm{\mathrm{Var}}}(\relax(Y \mid G))}{\mathop{\mathrm{\mathrm{Var}}}(Y)}\tag{5}\] and \[\label{def:h2} h^2 = \frac{\mathop{\mathrm{\mathrm{Var}}}(G^T \beta_0)}{\mathop{\mathrm{\mathrm{Var}}}(Y)},\tag{6}\] where \(G^T \beta_0\) is the least-squares projection of \(\relax(Y \mid G)\) onto linear functions of \(G\), that is, \[\beta_0 = \arg\min_{\beta \in \mathbb{R}^d} \mathop{\mathrm{\mathrm{Var}}}(\relax(Y \mid G) - G^T \beta).\] It follows that \(0 \leq h^2 \leq H^2 \leq 1\). Here we deliberately exclude any measured environmental factors \(X\) even though they are often used in practical estimation methods. This is because the exact regression model often changes the definition of \(H^2\) and \(h^2\) and makes the comparison more difficult.

It can be shown that \(H^2\le \xi_l'\) if \(G \perp \!\!\! \perp X\) and \(H^2 = \xi_l'\) if \((Y,G) \perp \!\!\! \perp X\) (see Lemma 2). However, \(H^2\le \xi_l'\) is not always true; for example, when \(G = X\) with probability 1, we have \(\xi_l' = 0\) but \(H^2\) can still be positive.

max width=

Table [tab:heritability1] compares counterfactual heritability and its bounds with narrow-sense and broad-sense heritability in a number of examples. To simplify the expressions, we consider two independent, binary genotypes \(G_1,G_2 \sim \text{Bern}(0.5)\) and some IID normally distributed observed environmental factor \(X\) and unobserved factors \(E_1,E_2\) with mean \(0\) and variance \(1/4\). Table [tab:heritability2] further provides the numeric values of different notions of heritability in a variety of examples with two causal genetic variants \(G_1,G_2\) (the exact settings are described in the Online Supplement).

It is easy to see from Table [tab:heritability1] that all six notions of heritability are different. As another instructive example, in the two settings where \(Y(g) = \beta g_1 E_2 + E_1\) and \(Y(g) = \beta g_1 X + E_1\), the counterfactual heritability \(\xi\) is the same but the interval \([\xi_l, \xi_u']\) is tighter in the latter case because the environmental factor \(X\) is known.

In Table [tab:heritability2], we observe that when there is no interaction between gene and unknown environmental factors in generating \(Y\) (as in the examples (1),(2),(4),(5)), the true counterfactual heritability \(\xi\), its tight lower bound \(\xi_l'\), and its lower bound \(\xi_l\) are the same. This is because the outcomes \(Y(G)\) and \(Y(G')\) are comonotonic. Additionally, we observe that the narrow-sense heritability \(h^2\) and broad-sense heritability \(H^2\) are both no greater than \(\xi\). This is because \(h^2\) and \(H^2\) only consider genetic effect as an independent component in determining the outcome \(Y\), and do not account for the interaction effect between gene and environment. Furthermore, when the known covariate \(X\) is not included in the model, the upper bound \(\xi_u'\) becomes trivial (\(\xi_u'=1\)).

max width=

3.2 Counterfactual heritability of one generation and comparison with twin-study heritability↩︎

Our definition of counterfactual heritability requires a suitable choice of environmental confounders \(X\), so \(\xi\) can be understood as the proportion of phenotypic variation that is attributable (in a counterfactual sense) to genetic variation given \(X\). In a population that evolves naturally, such as humans, it is reasonable to let \(X\) collect all haplotypes of the parents inherited from grandparents, so \(\xi\) essentially represents the counterfactual heritability in one generation due to the randomness of meiosis.

The next example demonstrates, in a setting similar to Example 1, a close connection between this notion of one-generation counterfactual heritability and twin-study heritability that is often calculated by Falconer’s formula, \[\label{def:twin} h^2_{\text{twin}} = 2(\rho_{\text{MZ}}-\rho_{\text{DZ}}),\tag{7}\] where \(\rho_{\text{MZ}}\) and \(\rho_{\text{DZ}}\) are the correlations between monozygotic and dizygotic twins.

Example 2 shows that we can understand Falconer’s formula as (twice) the one-generation counterfactual heritability in twin studies under an additive model. In this sense, counterfactual heritability provides a natural generalization of Falconer’s formula to other study designs and non-additive models.

3.3 Distinguishing notions of heritability by the counterfactual comparison↩︎

To further illustrate the concept of counterfactual heritability of one generation, it is useful to consider the simple graphical model in Figure [fig:direct-indirect-a] which shows that parental haplotype \(X\) has a direct and indirect effect (via offspring genotypes) on the offspring phenotype \(Y\). The unconfoundedness assumption, that is \(G \perp \!\!\! \perp\{Y(g): g \in \mathcal{G}\} \mid X\) in [3 ], holds under a standard interpretation of Figure 1 as a nonparametric structural equation model with independent errors [28]. Let \(G'\) be an independent realization of \(G\) given \(X\), \(X'\) be an independent realization of \(X\), and \(G''\) be an independent realization given \(X'\). In this case, the one-generation counterfactual heritability is \[\xi_{\text{fraternal}} = \frac{\mathop{\mathrm{\mathrm{Var}}}\{Y(X,G)-Y(X,G')\}}{2\mathop{\mathrm{\mathrm{Var}}}(Y(X,G))},\] and the population counterfactual heritability is \[\xi_{\text{unrelated}} = \frac{\mathop{\mathrm{\mathrm{Var}}}\{Y(X,G)-Y(X',G'')\}}{2\mathop{\mathrm{\mathrm{Var}}}(Y(X,G))}.\] The data generating mechanisms for the counterfactuals in these definitions are shown in panels (b-d) of Figure 1. The one-generation \(\xi_{\text{fraternal}}\) compares an individual in panel (b) of Figure 1 to their within-family fraternal counterfactual in panel (c), while the population-based \(\xi_{\text{unrelated}}\) compares an individual in panel (b) to their population (unrelated) counterfactual in panel (d). In reality, because the haplotype \(X\) of the parents is further inherited from their parents, and so on, it is perhaps more appropriate to think about \(\xi_{\text{unrelated}}\) as “infinite-generation” heritability.

Table [tab:heritability3] provides a numeric comparison of different notions of heritability in a specific with-in family example. Since \(h^2\) and \(H^2\) do not account for family structure, their values remain the same as in Table [tab:heritability2]. For models without gene-environment interaction (as seen in rows analogous to those in Table [tab:heritability2]), the true counterfactual heritability \(\xi\), the tight lower bound \(\xi'_l\), and the lower bound \(\xi_l\) are again identical. Due to the additional information of parents’ genotypes coded in \(X\), the upper bound \(\xi'_u\) is strictly less than \(1\). Additionally, since now \(G\) and \(G'\) are conditionally independent given the parents’ genotype, the variance of \(Y(G) - Y(G')\) is smaller than the marginally independent case in Table [tab:heritability2]. Therefore, the values of \(\xi\) in Table [tab:heritability3] are smaller than in Table [tab:heritability2]. The exact settings that generated Table [tab:heritability3] can be found in the Online Supplement.

max width=

3.4 Comparison with relatedness disequilibrium regression↩︎

The relatedness disequilibrium regression (RDR) is a recent approach that defines heritability through random segregation and thus bears conceptual similarity to counterfactual heritability [16]. A key observation in [16] is that traditional heritability estimator can be biased by confounding environmental factors. In particular, even though traditional kinship methods based on twins and siblings are robust to additive family effects, they can still be biased by sibling interference effects, meaning one’s genotype may influence the phenotype of their sibling (this is called “reciprocal genetic effects” in [16]). Consider the following example.

The definition of \(h^2_{\text{RDR}}\) in this example motivates a counterfactual quantity \[\xi_{\text{adopted}} = \frac{\mathop{\mathrm{\mathrm{Var}}}\{Y(X,G)-Y(X,G'')\}}{2 \mathop{\mathrm{\mathrm{Var}}}(Y(X,G))},\] where \(G''\) is the genotype of the offspring in a unrelated family with parental haplotype \(X'\) independent of \(X\); see panels (d) and (e) in Figure 1. So in the definition of \(\xi_{\text{adopted}}\), an individual is compared with their adopted counterfactual. It is not difficult to show that \(h^2_{\text{RDR}} = \xi_{\text{adopted}}\) in the setting of Example 3. However, even though it may be sensible to compare to an adopted counterfactual in some scenarios, this definition seems to be flawed: in the setting of Example 3, \(h^2_{\text{RDR}}\) or \(\xi_{\text{adopted}}\) can be larger than \(1\) when \(f_{\text{direct}}(G_1)\) and \(f_{\text{indirect}}(G_2)\) are negatively correlated.

3.5 Comparison with plant breeding heritability↩︎

In the plant breeding literature, heritability is usually defined on an entry-mean basis because the phenotype of interest is usually an aggregated value over replicated genotypes in experiments [7], [8]. The definition of counterfactual heritability can be easily extended to this case by considering the counterfactual of the aggregated phenotype, as illustrated below.

Let \(Y(g, x)\) denote the potential outcome of a phenotype of the plant when the genotype is set to \(g \in \{1,\dots,n_g\}\) and environment is set to \(x \in \{1,\dots,n_x\}\). A plant breeding experiment often measures \(Y_i = Y(G_i, X_i), i=1,\dots,n\) under a factorial design such that \[\sum_{i=1}^n 1_{\{G_i = g, X_i = e\}} = n_r,~\text{for all}~g~\text{and}~x,\] where \(n_r\) is the number of replications (so \(n = n_g n_x n_r\)). Denote the mean phenotype for all plants with genotype \(g = 1,\dots,n_g\) and environment \(x = 1,\dots,n_x\) as \[\bar{Y}(g, x) = \frac{1}{n_r} \sum_{i=1}^n 1_{\{G_i = g, X_i = x\}} Y_i(g,x).\] One can think of \(\bar{Y}(g, x)\) as the potential outcome of the mean phenotype \(\bar{Y}\) when we set the genotype as \(g\) and environment as \(x\) in an experiment with \(n_r\) replications. Furthermore, denote the mean phenotype for all plants with genotype \(g = 1,\dots,n_g\) as \[\label{eq:y-bar-g} \bar{Y}(g) = \frac{1}{n_x} \sum_{x=1}^{n_x} \bar{Y}(g, x) = \frac{1}{n_r n_x} \sum_{i=1}^n 1_{\{G_i = g\}} Y_i,\tag{8}\] which can be thought as the potential outcome of the mean phenotype \(\bar{Y}\) when we set the genotypes as \(g\) in \(n_r n_x\) replications over \(n_x\) environmental conditions.

Let \(G\) and \(G'\) be two independent random variables that are both uniformly distributed over \(\{1,\dots,n_g\}\), and let \(X\) be uniformly distributed over \(\{1,\dots,n_x\}\). Following equation [2 ], the counterfactual heritability of the mean phenotype \(\bar{Y}\) is given by \[\xi = \frac{\mathop{\mathrm{\mathrm{Var}}}(\bar{Y}(G) - \bar{Y}(G'))}{2 \mathop{\mathrm{\mathrm{Var}}}(\bar{Y}(G))}.\] We next compare this with an existing definition of plant breeding heritability assuming \(G\) and \(X\) are independent by design.

The above example shows that existing definitions of heritability in breeding experiments based on variance components may not have a clear causal interpretation even under relatively simple models. Adopting a counterfactual perspective can help researchers state their assumptions (e.g.random effects versus fixed effects) more explictly and interpret their estimand more easily.

4 Discussion↩︎

This article provides a unifying, counterfactual perspective on heritability. As a major advantage, our notion of counterfactual heritability does not depend on specific statistical modeling assumptions. Our results and comparisons with existing notions of heritability highlight the importance of clarifying the causal structural assumptions and counterfactual comparisons in reasoning with heritability. This is in line with in the ongoing “causal inference revolution” in statistics and other disciplines, where the focus of scientific research is shifted from statistical models and estimators to study designs and estimands/parameters.

A key requirement in defining counterfactual heritability is to find background environmental or genetic factors \(X\) such that the genotype \(G\) is randomized conditional on \(X\). The appropriate choice of \(X\) depends on the design of the study. In plan breeding experiments, the natural choice of \(X\) is different and may include control variables used in the design. In observational human genetics study, the most natural choice of \(X\) is the parental haplotypes, leading to the notion of one-generation counterfactual heritability discussed in Section 3.2. Reassuringly, in twin studies this basically reduces to familiar Falconer’s formula under additive models. As discussed earlier, the notion of one-generation heritability is closely related to recent discussion in human genetics that highlight the distinction between (direct) genetic effects and (indirect) nurturing effects [16], [21]. Our framework can be broadly viewed as a nonparametric, counterfactual extension to those efforts.

Counterfactual heritability can be further extended to a notion of counterfactual explainability of other factors. For example, the same definition can be used to define the counterfactual explainability of specific genetic factors and can be further extended to define the explainability of their interactions; similarly, it can be used to define the heritability of environmental factors and extended to define the explainability of gene-environment interaction. This is explored by [42] beyond the genetic heritability problem and is closely related to the functional analysis of variance [43], [44].

Our notion of counterfactual heritability is naturally defined for continuous (“quantitative”) and binary (“qualitative”) traits. For binary traits, it is also possible to define its counterfactual heritability through a latent, continuous liability index. Moreover, one can naturally generalize [2 ] to define counterfactual the genetic correlation between two phenotypes \(Y\) and \(Z\) as \[\rho = \frac{\mathop{\mathrm{\mathrm{Cov}}}(Y(G) - Y(G'), Z(G)- Z(G'))}{2 \sqrt{\mathop{\mathrm{\mathrm{Var}}}(Y(G)) \mathop{\mathrm{\mathrm{Var}}}(Z(G))}}.\] We leave the investigation of these quantities (e.g., their bounds and interpretations) for future work.

When discussing heritability, it is important to bear in mind that nature vs.nurture is a complex question, and heritability is just one number and provides inherently limited information. The fact that counterfactual heritability is only partially identified (can only be bounded using empirical data) reflects our inherent uncertainty about the nature vs.nurture decomposition, due to the impossibility of ascertaining the interaction strength between genes and unobserved environmental factors. We believe any notion of heritability should be interpreted against this backdrop.

Acknowledgment↩︎

Qingyuan Zhao was supported in part by the Engineering and Physical Sciences Research Council (EP/V049968/1). Hongyuan Cao is supported in part by NSF DMS 2311249. The authors thank Yunpeng Wang and Matthew Stephens for their comments on oral presentations of this work.

5 Proofs↩︎

5.1 Proofs with counterfactual heritability within population↩︎

5.2 Derivation of the result in Table [tab:heritability1]↩︎

5.3 Relationship between twin heritability and the counterfactual heritability↩︎

We assume that \(g_m,g_f\in \{0,1,2\}\) are the mother and father’s genes which are generated from the population. We assume that \(Y_1(g)=f(g)+\epsilon_1\), and \(Y_2(g)=f(g)+\epsilon_2\), where \[f(g)=\begin{cases}-a&g=0\\d&g=1\\a&g=2\end{cases}\] and \(\epsilon_1,\epsilon_2\) are i.i.d. The children’s genes are generated based on the parents’ genes. The detail can be found in Table 1. Denote \(G_1\) and \(G_2\) as the genotypes of participants. For MZ twin \(G_1=G_2=G\) a.s., and for DZ twin \(G_1\perp G_2\mid g_m,g_f\). It can be shown that \[\mathop{\mathrm{\mathrm{Var}}}(f(G))=\sigma_a^2+\sigma_d^2,\] and \[\label{eq:full95sib} \mathop{\mathrm{\mathrm{Cov}}}(f(G_1),f(G_2))=\frac{1}{2}\sigma_a^2+\frac{1}{4}\sigma_d^2.\tag{10}\]

The covariance between MZ twins is \[\begin{align} \rho_{MZ}&=\frac{\mathop{\mathrm{\mathrm{Cov}}}(Y_1(G),Y_2(G))}{\mathop{\mathrm{\mathrm{Var}}}(Y)}\\&=\frac{\mathop{\mathrm{\mathrm{Var}}}(f(G))}{\mathop{\mathrm{\mathrm{Var}}}(Y)}\\&=\frac{\sigma_a^2+\sigma_d^2}{\mathop{\mathrm{\mathrm{Var}}}(Y)}, \end{align}\] The covariance between DZ twins is \[\begin{align} \rho_{DZ}&=\frac{\mathop{\mathrm{\mathrm{Cov}}}(Y_1(G_1),Y_2(G_2))}{\mathop{\mathrm{\mathrm{Var}}}(Y)}\\&=\frac{\mathop{\mathrm{\mathrm{Cov}}}(f(G_1),f(G_2))}{\mathop{\mathrm{\mathrm{Var}}}(Y)}\\&=\frac{\frac{1}{2}\sigma_a^2+\frac{1}{4}\sigma_d^2}{\mathop{\mathrm{\mathrm{Var}}}(Y)}. \end{align}\]

The twin heritability is \[h^2_{\text{twin}}=2(\rho_{MZ}-\rho_{DZ})=(\sigma_a^2+\frac{3}{2}\sigma_d^2)/\mathop{\mathrm{\mathrm{Var}}}(Y).\]

The counterfactual heritability is defined as \[\begin{align} \xi&=\frac{\mathop{\mathrm{\mathrm{Var}}}(Y_1(G_1)-Y_1(G_2))}{2\mathop{\mathrm{\mathrm{Var}}}(Y)}\\&=\frac{\mathop{\mathrm{\mathrm{Var}}}(f(G_1)-f(G_2))}{2\mathop{\mathrm{\mathrm{Var}}}(Y)}\\&=\frac{\mathop{\mathrm{\mathrm{Var}}}\{f(G_1)\}-\mathop{\mathrm{\mathrm{Cov}}}\{f(G_1),f(G_2)\}}{\mathop{\mathrm{\mathrm{Var}}}(Y)}\\&=\frac{\frac{1}{2}\sigma_a^2+\frac{3}{4}\sigma_d^2}{\mathop{\mathrm{\mathrm{Var}}}(Y)}. \end{align}\]

In fact, the following lemma holds for the relationship between twin heritability \(h^2_{\text{twin}}\) and \(\xi_l'\).

Here is the derivation of the equation 10 . Note \(G_m\) and \(G_f\) are the mother’s and father’s genotype, \(G_1\) and \(G_2\) are the genotypes of two siblings. Define \(g(G)=f(G)-\relax(f(G))\) to be the centralized genetic effect. We have \[p^2g(AA)+2pqg(Aa)+q^2g(aa)=0.\] Note that \(G_1\perp G_2\mid G_m,G_f\) and \(G_1\overset{d}{=}G_2\mid G_m,G_f\). We obtain that \[\begin{align} &\mathop{\mathrm{\mathrm{Cov}}}(f(G_1),f(G_2))\\=&\relax(g(G_1)g(G_2)) \\=&\relax[\relax\{g(G_1)g(G_2)\mid G_m,G_f\}] \\=&\relax[\relax\{g(G_1)\mid G_m,G_f\}\relax\{g(G_2)\mid G_m,G_f\}] \\=&p^4 \{g(AA)\}^2+4p^3q \{\frac{1}{2} g({AA})+\frac{1}{2}g({Aa})\}^2+2p^2q^2\{g({Aa})\}^2\\&\quad+4p^2q^2\{\frac{1}{4} g({AA})+\frac{1}{2}g({Aa})+\frac{1}{4} g({aa})\}^2\\&\quad+4pq^3\{\frac{1}{2}g({Aa})+\frac{1}{2} g({aa})\}^2+q^4\{g({aa})\}^2 \\&\quad -\{p^2g(AA)+2pqg(Aa)+q^2g(aa)\}^2 \\=&p^3q\{g({AA})-g({Aa})\}^2+pq^3\{g({aa})-g({Aa})\}^2\\&\quad+p^2q^2\{\frac{\{g(AA)+2g(Aa)+g(aa)\}^2}{4}-2g(Aa)^2-2g(AA)g(aa)\} \\=&\quad p^2q^2\left\{\frac{\{x+y+4g(Aa)\}^2}{4}-2g^2(Aa)-2\{x+g(Aa)\}\{y+g(Aa)\}\right\}\\&\quad +pq(p^2x^2+q^2y^2) \\=&pq[p^2x^2+q^2y^2+pq\{\frac{(x+y)^2}{4}-2xy\}] \end{align}\] where \(x=g(AA)-g(Aa)=a-d\), and \(y=g(aa)-g(Aa)=-a-d\). It can be simplified by \[\begin{align} &\mathop{\mathrm{\mathrm{Cov}}}(f(G_1),f(G_2))\\=&pq\{p^2(a-d)^2+q^2(-a-d)^2\}\\&\quad+p^2q^2\{\frac{(a-d-a-d)^2}{4}-2(a-d)(-a-d)\} \\=&pq\{p^2(a-d)^2+q^2(a+d)^2+pq(2a^2-d^2)\} \\=&pq\{(p^2+q^2+2pq)a^2+2(q^2-p^2)ad+(p^2+q^2-pq)d^2\} \\=&pq\{a^2+2(q-p)ad+(p^2+q^2-pq)d^2\} \\=&pq\{a+(q-p)d\}^2+p^2q^2d^2 \\=&\sigma_a^2/2+\sigma_d^2/4. \end{align}\]

5.4 Simulation details for Table [tab:heritability2] and [tab:heritability3]↩︎

In this section, we provide details for deriving results in Table [tab:heritability2] and [tab:heritability3]. We first introduce some notations. Consider two independent genes with minor allele frequency (MAF) \(0.5\). Denote parents’ genotypes as \(G_p=(G_{1f},G_{1m},G_{2f},G_{2m})\), where the subscripts \(f\) and \(m\) indicate “father” and “mother”, respectively. Each element of \(G_p\) is i.i.d following distribution \(\mathcal{P}_G\), with probabilities \(0.25\), \(0.5\), and \(0.25\) for values \(0\), \(1\), and \(2\), respectively. Conditional on parents’ genotype, we can obtain children’s genotype \(G=(G_1,G_2)\) from Table 1. Averaging over the population, marginally, \(G_1\sim \mathcal{P}_G\) and \(G_2\sim \mathcal{P}_G.\) Let \(X\) be the covariate, which contains environmental variables in Table [tab:heritability2] or parents’ phenotypes in Table [tab:heritability3]. In Table [tab:heritability2], \(X\sim P_X\) independently of \(G\) as an external environmental effect, where \(\mathcal{P}_X(0),\mathcal{P}_X(0.5),\mathcal{P}_X(1),\mathcal{P}_X(1.5)\) and \(\mathcal{P}_X(2)\) are \(1/16,4/16,6/16,4/16\), and \(1/16,\) respectively. The error terms, \(E_1\sim N(\mu=0,\sigma^2=0.5)\) and \(E_2\sim P_X\). Moreover, \(E_1\) and \(E_2\) are independent of \(G\) and \(G_p,\) and \(G'=(G'_1,G'_2)\) is an i.i.d copy of \(G\). We calculate heritability based on different definitions (e.g. narrow, broad, counterfactual) and obtain the lower and upper bounds of the counterfactual heritability. The results are summarized in Table [tab:heritability2].

In Table [tab:heritability3], the setting is exactly the same as that of Table [tab:heritability2], except that the covariate \(X\) is defined as \(X=(G_{1f}+G_{1m})/2\), representing the average parental genotype (i.e., a proxy for parental phenotype). This introduces a correlation between children’s genotype and environmental factors. Additionally, \(G'\) is an i.i.d copy of \(G\) conditional on \(G_p\), mimicking heritability within one generation. We also calculate \(h^2_{\text{twin}}\) defined in (7 ). All results are presented in Table [tab:heritability3].

In the simulation, we know the data-generating process, which allows us to compute the true values of heritabilities in all cases except \(\xi_l\). When we calculate \(\xi\), \(\xi_l'\), \(\xi_u'\) and \(h^2_{\text{twin}}\), which are defined in (2 ),(?? ),(?? ) and (7 ), it suffices to evaluate \(\relax(Y\mid G,G_p)\), \(\mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)\), \(\relax(Y(G)-Y(G')\mid G,G',G_p)\) and \(\mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)\). The results of these expressions are provided in subsection 5.4.1. The calculations of narrow-sense heritability \(h^2,\) defined in (6 ), and broad-sense heritability \(H^2\), defined in (5 ), involve calculating \(\relax(Y\mid G)\), which can be obtained from \(\relax(Y\mid G, G_p)\) and the posterior probability \(P(G_p\mid G)\).

We next elaborate on the calculation of \(\xi_l\), defined in (?? ). In the binary case, the quantiles \(F^{-1}_{G,X}(U)\) and \(F^{-1}_{G',X}(U)\) (as defined in Proposition 1) can be calculated directly by \(P(Y=1\mid G,X)\). In the continuous case, when the tight lower bound \(\xi_l\) does not have the same value as \(\xi_l'\) (e.g. scenarios \((3)\) of Table [tab:heritability3]), the calculation of \(\xi_l\) involves the quantile functions \(F^{-1}_{G,X}(U)\) and \(F^{-1}_{G',X}(U)\). These are the quantiles of a mixture of normal distributions, and we resort to a numerical approximation. Specifically, for each fixed combination of the parents’ genotype (\(G_p\)), we can generate two sets of potential outcomes: \(Y_{1i}(G)\mid G_p\), and \(Y_{2i}(G')\mid G_p\), for \(i=1,\cdots ,N\), where \(N=1,000,000\). We then sort each sample to obtain the order statistics \(Y_{1(i)}\), and \(Y_{2(i)}\). By law of large numbers, \[\frac{1}{N}\sum\limits_{i=1}^N (Y_{1(i)}-Y_{2(i)})^2\to E\big[\big(F^{-1}_{G_p,G}(U)-F^{-1}_{G_p,G'}(U)\big)^2\mid G_p,G,G'\big].\] The tight lower bound \(\xi_l\) can then be estimated from (?? ).

5.4.1 Conditional expectations and conditional variances↩︎

Here are the results of \(\relax(Y\mid G,G_p)\), \(\mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)\), \(\relax(Y(G)-Y(G')\mid G,G',G_p)\) and \(\mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)\) in different settings of Table [tab:heritability3]. Denote \(\Phi(\cdot)\) as the cdf of a standard normal distribution.

When \(Y=G_1+G_2+E_1,\) \[\begin{align} \relax(Y\mid G,G_p)&=G_1+G_2, \\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\mathop{\mathrm{\mathrm{Var}}}(E_1)=0.5,\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=(G_1+G_2)-(G_1'+G_2'),\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=0. \end{align}\]

When \(Y=G_1G_2+E_1,\) \[\begin{align} \relax(Y\mid G,G_p)&=G_1G_2, \\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\mathop{\mathrm{\mathrm{Var}}}(E_1)=0.5,\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=G_1G_2-G_1'G_2',\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=0. \end{align}\]

When \(Y=G_1E_2+E_1,\) \[\begin{align} \relax(Y\mid G,G_p)&=G_1\relax(E_2)=G_1,\\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\mathop{\mathrm{\mathrm{Var}}}(E_1)+G_1^2\mathop{\mathrm{\mathrm{Var}}}(E_2)=0.5+0.25G_1^2,\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=(G_1-G_1')\relax(E_2) = G_1 - G_1',\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=(G_1-G_1')^2\mathop{\mathrm{\mathrm{Var}}}(E_2) = 0.25(G_1 - G_1')^2. \end{align}\]

When \(Y=G_1+X+E_1,\) \[\begin{align} \relax(Y\mid G,G_p)&=G_1+X=G_1+\frac{G_{1f}+G_{2f}}{2},\\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\mathop{\mathrm{\mathrm{Var}}}(E_1)=0.5,\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=G_1-G_1',\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=0. \end{align}\]

When \(Y=G_1X+E_1,\) \[\begin{align} \relax(Y\mid G,G_p)&=G_1X=G_1\frac{G_{1f}+G_{1m}}{2},\\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\mathop{\mathrm{\mathrm{Var}}}(E_1)=0.5,\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=(G_1-G_1')X=(G_1-G_1')\frac{G_{1f}+G_{1m}}{2},\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=0. \end{align}\]

When \(Y=1\{G_1+G_2+E_1>0\}\), \[\begin{align} \relax(Y\mid G,G_p)&=\Phi(\sqrt2 (G_1+G_2)),\\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\relax(Y\mid G,G_p)(1-\relax(Y\mid G,G_p)),\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=\Phi(\sqrt2 (G_1+G_2))-\Phi(\sqrt2 (G_1'+G_2')),\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}\\&\quad\cdot(1-abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}). \end{align}\]

When \(Y=1\{G_1G_2+E_1>0\}\), \[\begin{align} \relax(Y\mid G,G_p)&=\Phi(\sqrt2 (G_1G_2)),\\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\relax(Y\mid G,G_p)(1-\relax(Y\mid G,G_p)),\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=\Phi(\sqrt2 (G_1G_2))-\Phi(\sqrt2 (G_1'G_2')),\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}\\&\quad\cdot(1-abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}). \end{align}\]

When \(Y=1\{G_1E_2+E_1>0\}\), \[\begin{align} \relax(Y\mid G,G_p)&=\relax_{E_2}\{\Phi(\sqrt2 G_1E_2)\} \\&=\frac{1}{16}\Phi(0)+\frac{4}{16}\Phi(\frac{\sqrt2}{2} G_1)+\frac{6}{16}\Phi(G_1)\\&\quad+\frac{4}{16}\Phi(\frac{3\sqrt2}{2} G_1)+\frac{1}{16}\Phi(2 G_1), \\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\relax(Y\mid G,G_p)(1-\relax(Y\mid G,G_p)),\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=\relax_{E_2}\{\Phi(\sqrt2 G_1E_2)\}-\relax_{E_2}\{\Phi(\sqrt2 G_1'E_2)\}\\ &=\frac{4}{16}T(1/2)+\frac{6}{16}T(1)+\frac{4}{16}T(3/2)+\frac{1}{16}T(2),\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=*, \end{align}\] where \[\begin{align} *=&\mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)\\=& \mathop{\mathrm{\mathrm{Var}}}\{\relax(Y(G)-Y(G')\mid G,G',G_p,E_2)\mid G,G',G_p\} \\&\quad+ \relax\{\mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p,E_2)\mid G,G',G_p\} \\=& \mathop{\mathrm{\mathrm{Var}}}_{E_2}(T\mid G,G',G_p) + \relax_{E_2}\{|T|(1-|T|)\mid G,G',G_p\}, \end{align}\] and \(T(E_2)=\Phi(\sqrt{2}G_1E_2)-\Phi(\sqrt{2}G_1'E_2).\) Here, the direct calculation shows that \[\begin{align} &\mathop{\mathrm{\mathrm{Var}}}_{E_2}(T\mid G,G',G_p) \\&=\frac{4}{16}T^2(1/2)+\frac{6}{16}T^2(1)+\frac{4}{16}T^2(3/2)+\frac{1}{16}T^2(2)\\&\quad-\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}^2, \end{align}\] and \[\begin{align} &\relax_{E_2}\{|T|(1-|T|)\mid G,G',G_p\} \\&=\frac{4}{16}|T(1/2)|\{1-|T(1/2)|\}+\frac{6}{16}|T(1)|\{1-|T(1)|\}\\&\quad+\frac{4}{16}|T(3/2)|\{1-|T(3/2)|\}+\frac{1}{16}|T(2)|\{1-|T(2)|\}. \end{align}\]

When \(Y=1\{G_1+X+E_1>0\}\), \[\begin{align} \relax(Y\mid G,G_p)&=\Phi(\sqrt2 (G_1+X))=\Phi(\sqrt2 (G_1+\frac{G_{1f}+G_{1m}}{2})),\\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\relax(Y\mid G,G_p)(1-\relax(Y\mid G,G_p)),\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=\Phi(\sqrt2 (G_1+X))-\Phi(\sqrt2 (G_1'+X))\\ &=\Phi\{\sqrt{2}(G_1+\frac{G_{1f}+G_{1m}}{2})\}-\Phi\{\sqrt{2}(G_1'+\frac{G_{1f}+G_{1m}}{2})\},\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}\\&\quad\cdot(1-abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}). \end{align}\]

When \(Y=1\{G_1X+E_1>0\}\), \[\begin{align} \relax(Y\mid G,G_p)&=\Phi(\sqrt2 G_1X)=\Phi(\sqrt2 G_1\frac{G_{1f}+G_{1m}}{2}),\\ \mathop{\mathrm{\mathrm{Var}}}(Y\mid G,G_p)&=\relax(Y\mid G,G_p)(1-\relax(Y\mid G,G_p)),\\ \relax(Y(G)-Y(G')\mid G,G',G_p)&=\Phi(\sqrt2 G_1X)-\Phi(\sqrt2 G_1'X)\\ &=\Phi(\sqrt{2}G_1\frac{G_{1f}+G_{1m}}{2})-\Phi(\sqrt{2}G_1'\frac{G_{1f}+G_{1m}}{2}),\\ \mathop{\mathrm{\mathrm{Var}}}(Y(G)-Y(G')\mid G,G',G_p)&=abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}\\&\quad\cdot(1-abs\{\relax(Y(G)-Y(G')\mid G,G',G_p)\}). \end{align}\]

References↩︎