Using the phenotype differences model to identify genetic effects in samples of partially genotyped sibling pairs

The phenotype differences model recovers estimates of genetic effects that are free of environmental confounding when two key assumptions hold: random assignment of genotype within families, which is an assumption also required by fixed effects, and equal population variance of the genetic predictor in the genotyped and ungenotyped siblings, an assumption unique to phenotype differences.

Let $g_{1j}$ have population variance ${\sigma }_1^2$, and let $g_{2j}$ have population variance ${\sigma }_2^2$. Thus, the assumption of phenotype differences is that ${\sigma }_1^2={\sigma }_2^2$. Because we do not directly observe $g_{2j}$, this assumption is inherently untestable. However, if we know that we observe the genotype of a random sibling (i.e., random sampling of genotype within families), then this assumption is trivially met.

When $g_{\mathit{ij}}$ is a normally distributed polygenic index, there exists no mechanical mean-variance dependence. In this case, average differences in genetic characteristics between the genotyped sibling and the ungenotyped siblings are not inherently a problem; this is because the existence of genetic characteristics that linearly increase or decrease an individual’s likelihood of being the genotyped (versus ungenotyped) sibling would not distort the variance of $g_{1j}$ (compared to $g_{2j}$) and therefore will not violate our assumption. However, nonlinear and/or nonmonotonic selection into genotyping may impact the variance of $g_{1j}$ (relative to $g_{2j}$) and therefore induce violations.

When $g_{\mathit{ij}}$ is the number of major alleles at a single-nucleotide polymorphism (i.e., taking the value of 0, 1, or 2), such as in GWAS, mean differences likely entail variance differences; therefore, systematic differences in the allele frequencies across the genotyped and ungenotyped sibling will most likely induce violations of the equal variance assumption. When this key assumption is not met, estimates of the true genetic effect become biased as a function of the degree of variance discordance. Beginning from SI Appendix, Eq. S18 of the phenotype differences consistency proof, we observe that:

When $\rho =\frac{1}{2}$, this reduces to:

We empirically test the equal variance assumption in two different ways. First, while it is not possible to test whether a difference in genetic predictor variance exists between individuals who survived until the year that genotyping began (2006) and those who did not, we can artificially “move back” the genotyping date. Specifically, we test whether the 47 polygenic indices have equal variance among genotyped individuals who survived until 2020 ($N=$ 5,406) and among genotyped individuals who died prior to 2017 ($N=$ 2,302) using a Brown–Forsythe test (46, 47). Across the 47 separate tests, we observe no $P<0.05$, consistent with the lack of systematic differences in genetic predictor variance between genotyped individuals who did versus did not survive until 2020.

Second, while we definitionally do not observe the genetic predictor of both siblings in our One Genotype sample, we often do observe phenotype variables for both siblings. SI Appendix, Table S2 displays results from Brown–Forsythe equal variance tests using phenotypic data. We focus on 7 phenotypes that are largely continuously distributed to reduce the presence of mean-variance dependencies (which do not exist for normally distributed polygenic indices). Although one phenotype (neuroticism) has $P<0.05$, the magnitude of any phenotypic variance discordance that we observe between the genotyped and ungenotyped siblings is generally small and insignificant.

As we show in SI Appendix, section F, when the effects of $g_{\mathit{ij}}$ are small, fixed effects estimates will have asymptotically identical SEs to phenotype differences estimates derived from the same number of genotypes (although, when using phenotype differences, one typically has half as many genotypes per family). However, as the fraction of within-family outcome variation explained by $g_{\mathit{ij}}$ grows, phenotype difference provides less precise effect estimates than fixed effects per genotype. Specifically, for models fit on the same number of sibling pairs—that is, $N=M$, implying twice as many genotypes used for fixed effects than for phenotype differences—this decrease in comparative precision is governed by the following formula:

Here, $\varphi$ is the fraction of within-family variation in the $y_{\mathit{ij}}$ that is explained by within-family variation in $g_{\mathit{ij}}$:

The within-$R^2$ of the fixed effects model can provide a useful estimate of $\varphi$. For currently available genetic predictors, this reduction in precision due to $\varphi >0$ is relatively modest. For example, when $g_{\mathit{ij}}$ is the polygenic index for height—one of the most predictive scores in the SSGAC repository—$\widehat{\varphi }=0.11$ in the WLS Two Genotypes sample. Therefore, the comparative precision per genotype is $2\times \sqrt{\frac{1-0.11}{4-0.11}}=0.96$. That is, when fit on the same number of genotypes, fixed effects estimates of the effect of the height polygenic index on phenotypic height will have expected SEs that are 96% as large as the SEs of phenotype difference estimates. See SI Appendix, Fig. S2 for a graphical display of the relationship between $\varphi$ and the asymptotic ratio of SEs of the fixed effects and phenotype differences models. Note that, for ease of exposition, these derivations assume $\rho$ is known and equal to $\frac{1}{2}$.

When meaningful assortative mating exists, maternal and paternal genetic predictors become correlated with one another and $\rho \ne \frac{1}{2}$. In such a case, $\widehat{\rho }$ can be empirically estimated from a sample of fully genotyped sibling pairs, as we do with the WLS Two Genotypes Sample. It is important to note that while $\rho$ is fixed, $\widehat{\rho }$ is derived from a finite sample of $M$ fully genotyped sibling pairs and is therefore random. Thus, when fitting phenotype differences models with an estimated within-family correlation of genetic predictor, the uncertainty from $\widehat{\rho }$ propagates to $\widehat{\beta }$ and must be accounted for when constructing SEs. In SI Appendix, section E, we use the multivariate delta method (28, 29, 48) to show that the variance of ${\widehat{\beta }}^{\mathit{PD}}$ in this case is approximately:

where $N$ is the number of sibling pairs in the phenotype differences sample and $M$ is the number of sibling pairs used to estimate $\widehat{\rho }$. As can be seen in SI Appendix, Fig. S3, when $M>1,000$, this SE adjustment ends up being relatively small.

Occasionally, sibling pairs believe themselves to be full biological siblings but are, in actuality, only half-siblings. Because half-siblings have an expected within-family correlation of genetic predictor of $\frac{1}{4}$ (under random mating), sufficient prevalence of nonpaternity events can produce $\rho <\frac{1}{2}$. In such a case, phenotype differences estimates become biased away from 0. In the WLS Two Genotypes sample, approximately one-in-fifty sibling pairs who self-report being full biological siblings are actually half-siblings. If the variance of the overall distribution of $g_{\mathit{ij}}$ is identical for full and half-siblings, we would expect this frequency of nonpaternity events to induce $\rho \approx 0.495$. Such a small departure from $\rho =\frac{1}{2}$ will have only a trivial impact on phenotype differences estimates. Importantly, because misreporting individuals are unaware that they are half siblings, it is unlikely that there exist meaningful environmental differences between such half siblings that are correlated with paternal genotype, so confounding from population stratification and/or dynastic effects is unlikely.

The WLS is a survey based on a $\frac{1}{3}$ sample of all 1957 Wisconsin high school graduates and a randomly selected sibling of these graduates (22). The graduate respondents were originally empaneled with an in-person questionnaire at age 18 in 1957, which was followed by data collection at ages 25, 36, 54, 65, and finally 72 in 2012. The WLS includes a wide range of administrative and prospectively collected survey data from early life, adolescence, and adulthood. Genetic samples were assayed from saliva for a subsample of consenting WLS graduates and siblings. Genotyping was performed using the Illumina HumanOmniExpress 24 BeadChip arrays (Version 1/1.1; Illumina). We restrict our analytic sample to individuals of European ancestries because only these respondents have nonmissing polygenic index data in the SSGAC repository. However, due to the ancestral homogeneity of the WLS sample, this restriction is relatively inconsequential in practice; of the 9,012 genotyped WLS respondents, 8,927 were identified as European ancestries and have valid polygenic index data. While within-family models are robust to biases stemming from ancestry-related environmental confounding, including diverse ancestries together in a single phenotype differences model could produce inflated values of $\widehat{\rho }$ [due to ancestral assortative mating (49, 50)]. For the sake of open science, our main text analyses utilize the public-use release of the WLS data in which it is impossible to identify the 48 half siblings in the Two Genotypes Sample who believe themselves to be full siblings; however, in SI Appendix, Tables S3–S6, we rerun our analyses using the restricted-use data with these 48 cases removes and our findings our substantively identical.

The polygenic indices used in this study are drawn from Version 1.1 of the Social Science Genomics Association Consortium Polygenic Index Repository (23). All polygenic indices are standardized over the full sample of genotyped WLS graduates. A list of which polygenic indices are utilized in each analysis presented in this study can be found in SI Appendix, Table S1.

We align our phenotypic variables with those used in as Becker et al. (23); see supplementary table 12 of that paper for a list of the specific WLS survey items used. Though the repository contains 47 distinct polygenic indices, in the WLS, there only exist phenotype data as a subset of 30 traits. When multiple measurements were available, for variables such as depression, we first standardize the variable within each wave and then take the average for an individual across waves. For variables like educational attainment with multiple measurements, we take the maximum value across waves. Finally, all phenotypes are standardized over the full sample of genotyped WLS graduates.

The mortality data used in this study is derived from the National Death Index (51)—importantly, such data are available for all members of the WLS, regardless of how long they remained empaneled in the study. Mortality data are right censored in 2020; that is, we observe all deaths through December of 2019. Our key mortality variable is lifespan in years, which is calculated as the difference between death year and birth year (for individuals who are still alive, we use 2020 as their death year). Although this lifespan variable is right censored, because genotype is randomly assigned within families, we would not expect siblings differences in polygenic indices to be correlated with siblings differences in birth year. In SI Appendix, Fig. S1 and Dataset S4, we repeat our mortality analyses using a dichotomous indicator for survival to age 75. Sibling differences for these two mortality outcomes were residualized on a second-order polynomial of sibling difference in birth year.