Mixed Model with Correction for Case-Control Ascertainment Increases Association Power - PMC

Overview of Method

We improve upon standard mixed-model methods¹¹ by using a retrospective association-score statistic (LTMLM) computed from PMLs under the liability-threshold model. The improvement over previous approaches comes from appropriate modeling of case-control ascertainment. We consider all individuals simultaneously and incorporate prevalence information.

Our method consists of three steps. First, the GRM is calculated, and a corresponding heritability parameter is estimated, modeling the phenotype covariance of all individuals (see “Estimation of Heritability Parameter”). The heritability parameter is estimated with Haseman-Elston (H-E) regression on the observed scale and then transformed to the liability scale. Second, PMLs are estimated with a truncated multivariate normal Gibbs sampler (see “PMLs”). The PML of each individual is conditional on that individual’s case-control status, on every other individual’s case-control status, and on disease prevalence and liability-scale phenotypic covariance. Third, a χ² (1 degree of freedom) association-score statistic is computed on the basis of the association between the candidate SNP and the PML (see “LTMLM Association Statistic”).

The toy example in Figure 1 illustrates how genetic relatedness to a disease-affected case subject can increase an individual’s PML. In Figures 1A and 1B, we plot the distribution of liabilities in 10,000 unrelated individuals who were randomly ascertained and ascertained according to case-control status (for a disease with a prevalence of 0.1%), respectively. In Figures 1C and 1D, we plot the same distributions but condition on an individual’s having a genetic relatedness of 0.5 with a disease-affected case subject while assuming a liability-scale heritability of 1.0. In each case, the posterior distribution of liabilities (and hence the PML) is shifted upward. (The magnitude and direction of this effect would be different for an individual with a genetic relatedness of 0.5 with a control subject). Our main focus below is on much lower levels of genetic relatedness (identity by state) among many unrelated samples, but the same principles apply.

Estimation of Heritability Parameter

Mixed-model association statistics rely on the estimation of a heritability parameter. We note that this heritability parameter, which Kang et al.⁴ refer to as “pseudo-heritability,” is generally lower than the total narrow-sense heritability (h²) in datasets not dominated by family relatedness but can be larger than the heritability explained by genotyped SNPs (h_g²)¹⁷ in datasets with population structure or family relatedness. However, for ease of notation, we use the symbol h² to represent this heritability parameter. A list of all notation used below is provided in Table S1.

The goal is to test for association between a candidate SNP and a phenotype. We first consider a quantitative trait:

The phenotypic data (transformed to have mean 0 and variance 1) can be represented as a vector (φ) with values for each individual (i). Genotype values of candidate SNPs are transformed to a vector (x) with mean 0, variance 1, and effect size β. The quantitative trait value depends on the fixed effect of the candidate SNP (βx), the genetic random effect excluding the candidate SNP (u), and the environmental component (e). We extend to case-control traits via the liability-threshold model, in which each individual has an underlying, unobserved, normally distributed trait called the liability. An individual is affected by a disease if the liability exceeds a specified threshold (t), corresponding to disease prevalence¹⁶ (Figure S1).

Standard mixed-model association methods generally estimate h² from a GRM and phenotypes by using restricted maximum likelihood (REML).^4,11 Genotypic data (excluding the candidate SNP¹¹) are used for building a GRM:

where X is a matrix of non-candidate SNPs normalized to mean 0 and variance 1 and M is the number of SNPs. We estimate h² by using H-E regression followed by a transformation to the liability scale. We obtain the H-E regression estimate by regressing the product of the case-control phenotypes on the off-diagonal terms of the GRM:^18,19

where π_i denotes the case-control status of individual i and $\hat{\Theta }$ _ik is the genetic relatedness of individuals i and k. This gives an estimate on the observed scale, and then the estimate is transformed to the liability scale:²⁰

where z is the height of the standard normal density ($\left(1/2\pi \right)e^{-t^2/2}$) at the liability threshold t, K is disease prevalence, and P is the proportion of case subjects in the sample.²⁰

Then, the variance between the individuals is modeled as the phenotypic covariance,

where $\hat{\Theta }$ is the N × N GRM, V is the phenotypic covariance, h² is the heritability parameter, and I is the identity matrix.

Using the phenotypic covariance matrix V, the liability is modeled as a multivariate normal distribution:

We note that we observe the case-control phenotypes of the individuals and not the continuous liabilities.

PMLs

We first consider the univariate PML (PML_uni), constructed independently for each individual; we generalize to the multivariate setting below. As described in Equations 11 and 12 in Lee et al.,²⁰ these correspond to the expected value of the liability and are conditional on the case-control status:

These values are calculated analytically in the univariate setting, and depending on case or control status, they can be thought of as the mean of a truncated normal distribution above or below, respectively, the liability threshold t.²⁰

We next consider the multivariate PML (PML_multi), estimated jointly across individuals. The PML_multi for each individual is conditional on that individual’s case-control status, on every other individual’s case-control status, and on their phenotypic covariance. The PML_multi is estimated with a Gibbs sampler, as in previous work¹⁵ (which focused on family relatedness and did not consider association statistics). The Gibbs sampler is an iterative algorithm that generates random variables from conditional distributions in order to avoid the difficult task of explicitly calculating the marginal density for each random variable.

For each individual in turn, the conditional distribution of the liability is calculated on the basis of all of the other individuals, and a new value is generated. The algorithm is as follows:

Initialization: for each individual j,

For each Markov Chain Monte Carlo (MCMC ) iteration n

• For each individual

  i

  • Sample φ_i from the constrained conditional univariate normal distribution

  • L(φ_i) ∼exp(–φ^TV⁻¹φ/2) and constraint φ_i ≥ t if π_i = 1, φ_i < t if π_i = 0

  • (where φ_≠i are fixed)

We use 100 burn-in iterations followed by 1,000 additional MCMC iterations. We estimate the PML_multi by averaging over MCMC iterations. We reduce the number of MCMC iterations needed via Rao-Blackwellization, which averages (across n iterations) the posterior means of the distributions from which each φ_i is sampled.

LTMLM Association Statistic

The LTMLM association statistic is calculated with the PML_multi. For simplicity, we first consider the case where the liability is known. We jointly model the liability and the genotypes by using a retrospective model, enabling appropriate treatment of sample ascertainment. We concatenate the two vectors (φ, x) and derive the joint likelihood for these combined terms. The covariance of φ and x between individuals i and k is

where Θ is the true underlying GRM from which genotypes are sampled. (We note that Θ, which is unobserved, is different from the GRM $\hat{\Theta }$ estimated from the data). The variance of (φ, x) as a function of effect size $\beta$ is

Thus,

where both of these matrices are 2N × 2N. (We note that the product of the matrices in Equations 10 and 11 is $\left(\begin{array}{cc}I+O\left({\beta }^2\right)& 0\\ 0& I+O\left({\beta }^2\right)\end{array}\right),$ whose difference from the identity contains only O(β²) terms.)

The joint likelihood of the liability and genotypes is distributed as a multivariate normal N(0,C(β)), and thus

Taking the derivative of the log-likelihood results in the score equation. The determinant of the matrix V does not have any terms linear in β, so the terms with V alone drop out when we take the derivative:

The marginal score statistic tests the null hypothesis, which is that the fixed effect of the candidate SNP is 0 (H₀: β = 0), against the alternative hypothesis (H_A: β ≠ 0). The denominator of the score statistic is the variance of the score evaluated under the null:

This leads to the score statistic

where Θ, the true underlying genetic relatedness of the individuals, can be approximated by the identity matrix in datasets of unrelated individuals.

In Equations 9, 10, 11, 12, 13, 14, and 15, the liability was assumed to be known for simplicity. We now consider a case-control trait, with unobserved liability, and derive the score function by using the observed case-control status, π, of each individual. When we return to the score function and conditioning on case-control status,

When we introduce the unobserved quantitative liability, φ, the score function can be rewritten in terms of the probability density of the liability:

where P(φ) is the probability density of the liability and E[φ| π] is the PML. It follows that an appropriate LTMLM score statistic is

Again, Θ can be approximated by the identity matrix in datasets of unrelated individuals; we note that this choice affects only a constant calibration factor (because the denominator is the same for each candidate SNP) and that other calibration options are available (see below). As with other association statistics, the LTMLM score statistic generalizes to non-normally distributed genotypes.^21–23 The overall computational cost of computing the LTMLM statistic is O(MN²) when M > N > number of iterations (Table S2). We have fixed the number of iterations at 100 burn-in iterations followed by 1,000 additional iterations.

We calculate the GRM via leave one chromosome out (LOCO) analysis, i.e., for each candidate SNP on a given chromosome, we calculate GRM by using all of the other chromosomes. This prevents deflation due to double counting of the candidate SNP as both a fixed effect and a random effect in the mixed model.^4,6,11

Simulated Genotypes and Simulated Phenotypes

We performed simulations both by using simulated genotypes and simulated phenotypes and by using real genotypes and simulated phenotypes (see below). Quantitative liabilities for each individual were generated from SNP effects and an environmental component. All simulations included M candidate SNPs and an independent set of M SNPs used for calculating the GRM (so that candidate SNPs were not included in the GRM). For each scenario, a set of 100 simulations were run. We set ten candidate SNPs and ten GRM SNPs to be causal in simulations with N = 1,000 samples and set 50 candidate SNPs and 50 GRM SNPs to be causal in simulations with N = 5,000 samples to ensure that causal SNPs had similar average χ² statistics independent of M and N. We then dichotomized the resulting quantitative liabilities on the basis of the liability threshold to categorize each individual as a case or control subject. Case-control ascertainment was performed and simulated 50% case and 50% control subjects. We compared Armitage trend test (ATT), logistic regression (LogR), mixed linear model (MLM), and LTMLM statistics (see Table 1). MLM statistics were computed with the genome-wide complex trait analysis (GCTA)-LOCO statistic described in Yang et al.,¹¹ and the heritability parameters were estimated with the GCTA software.²⁴ We evaluated performance by using average χ² statistics at causal, null, and all markers, λ_GC at all markers (median χ² divided by 0.455),²⁵ and the proportion of causal and null markers that were significant at the p value thresholds of 0.05, 0.001, 1 × 10⁻⁶, and 5 × 10⁻⁸.

In the primary analyses, we simulated individuals without population structure or linkage disequilibrium (LD) with N = 1,000 or 5,000 samples; M = 1,000, 5,000, or 50,000 SNPs; and prevalence K = 50%, 10%, 1%, or 0.1%. Genotypes were sampled from independent binomials with allele frequencies sampled from a uniform distribution over the values [0.1, 0.9]. In secondary analyses, we simulated population structure by simulating two populations with F_ST = 0.01 and allele frequencies drawn from beta distributions with parameters p(1 – F_ST)/F_ST and (1 – p)(1 – F_ST)/F_ST, based on ancestral allele frequency p sampled from a uniform distribution over the values [0.1, 0.9].

To test the impact of the generative distribution, we simulated the underlying distribution by using a logit model instead of a liability-threshold model. The ten causal candidate SNPs were simulated with alternating fixed-effect sizes of β = 0.4 or β = −0.4. Then, case-control phenotypes were generated from a binomial distribution where the probability of being a case subject was $\left(\text{case}\right)=1/\left(1+\text{exp}\left(-\left[c+\mathrm{\beta x}\right]\right)\right)$, shifted by the affine term (c) and based on the desired disease prevalence.

WTCCC2 Genotypes and Simulated Phenotypes

We also conducted simulations by using real genotypes from WTCCC2 to incorporate LD and realistic population structure. The WTCCC2 data contained 360,557 SNPs and 15,633 samples, as described previously.¹¹ Because the goal of the power study is to demonstrate a comparison of the statistics under case-control ascertainment, we used N = 1,000 samples (500 case and 500 control subjects) and simulated phenotypes with a prevalence of 50%, 25%, and 10%. The prevalence was restricted to a lower bound of 10% because of the limitation of having only 15,633 WTCCC2 samples for simulating case-control ascertainment. We computed ATT, LogR, MLM, and LTMLM statistics as described above.

WTCCC2 Genotypes and MS Phenotypes

Finally, we analyzed WTCCC2 individuals with ascertained case-control phenotypes for MS,¹¹ a disease with a prevalence of around 0.1%. As in previous work, we assume that the disease prevalence is known on the basis of external epidemiological literature.^20,26,27 For the WTCCC2 MS data, we used a threshold of 3.0, corresponding to a disease prevalence of 0.1%.²⁶ We computed ATT, LogR, MLM, and LTMLM statistics as described above. Although the underlying MS study was appropriately matched for ancestry,²⁸ the data made available to researchers included only pan-European case and UK control subjects. Thus, the WTCCC2 dataset shows a severe mismatch in ancestry of case and control subjects; this severe mismatch between case and control subjects is not representative of a typical GWAS. We thus restricted our primary analysis to 10,034 samples with only a moderate mismatch in ancestry, but we also performed analyses of unmatched and stringently matched datasets (Figure S2). The unmatched dataset contained 10,204 case and 5,429 control subjects. We performed matching by first calculating 20 principal components (PCs) in the full cohort and weighing the contribution of each PC on the basis of how much phenotype variance it explained in a multiple regression. A Euclidean distance over these 20 weighted dimensions was then computed for all pairs of individuals, and each case subject was greedily assigned the nearest unmatched control subject until no matched case-control pairs could be identified. Finally, any matched case-control pairs who were not within 6 SDs of the mean pairwise distance were removed as outliers, yielding the 5,017 case and 5,017 matched control subjects used in our primary analysis. We performed stringent matching by additionally removing any matched case-control pairs who were not within 2 SDs of the mean pairwise distance, yielding the 4,094 case and 4,094 matched control subjects used in our stringently matched analysis.

We compared association statistics at 75 published SNPs associated with MS.¹¹ We used a jackknife approach to assess the statistical significance of differences in association statistics by excluding each of the 75 published SNPs in turn.

Simulations: Simulated Genotypes and Simulated Phenotypes

We first conducted simulations by using simulated genotypes and simulated ascertained case-control phenotypes (see Material and Methods). Our main simulations involved unrelated individuals with no population structure, but the impact of population structure is explored below. We evaluated the power of ATT, LogR, MLM, and LTMLM. We report average χ² statistics at causal, null, and all markers and λ_GC at all markers in Table 2, and we report the proportion of causal and null markers that were significant at various p value thresholds in Table S3. For diseases with low prevalence, the LTMLM statistic outperformed the ATT, LogR, and MLM statistics. Improvements in average χ² statistics at causal markers (which are naturally interpreted as the increase in effective sample size) were larger than improvements in power to detect an association at a given p value threshold, most likely as a result of the variable (normally distributed) effect sizes in this simulation (see below). For LTMLM and MLM at disease prevalences of 0.1%, 26%, and 5%, improvements in average χ² statistics at causal markers were observed in simulations with 5,000 SNPs and 50,000 SNPs, respectively. Smaller improvements were observed at higher disease prevalences. Test statistics were well calibrated at null markers. Simulations at other values of M and N indicated that the magnitude of the improvement depends on the value of N/M (Tables S3 and S4). Simulations with population structure demonstrated similar results but also showed inflation in the ATT statistic as expected (Tables S5 and S6).

The MLM statistics were calculated with an h² parameter estimated via REML,⁴ but the LTMLM statistics were calculated with an h² parameter first estimated via H-E regression on case-control phenotypes and then transformed to the liability scale^18,20 (see Material and Methods). As case-control ascertainment became more severe, the H-E regression estimate of h² remained unbiased, whereas the variance-component estimate was severely downwardly biased even after transformation to the liability scale (Table 3 and Table S7), consistent with previous work (see Golan et al.²⁹ and Table S9 in Yang et al.¹¹). Population structure resulted in bias for both REML and H-E-regression estimates of h², although the bias for the REML estimates was higher (Table S8). These biases did not inflate LTMLM or MLM statistics under the null hypothesis (Tables S5 and S6). We note that previous work has shown that running MLM with the correct h² parameter does not ameliorate the loss in power for MLM.¹¹

We also evaluated performance in settings where the liability threshold (equivalently, the disease prevalence) was misspecified (Tables S9 and S10). The LTMLM statistic remained properly calibrated under the null and continued to outperform the MLM statistic given that the impact of misspecifying the liability threshold was small. Misspecifying the liability threshold led to bias in liability-scale heritability estimates as a result of the inaccurate conversion from the observed scale to the liability scale (Table S11).

Finally, we evaluated performance when phenotypes were generated with a logit model instead of a liability-threshold model; for simplicity, we used a fixed effect size for causal candidate SNPs (see Material and Methods). At low disease prevalence, we observed improvements for LTMLM both in average χ² statistics at causal markers (Table S12) and in power to detect an association at a given p value threshold (Table S13); improvements in power depend heavily on the distribution of causal effect sizes and are larger in simulations with fixed causal effect sizes than in simulations with variable causal effect sizes.

Simulations: WTCCC2 Genotypes and Simulated Phenotypes

We next conducted simulations by using real WTCCC2 genotypes and simulated phenotypes from ascertained case and control individuals (see Material and Methods).^11,28 For a given value of M (M SNPs for calculating the GRM plus M candidate SNPs for a total of 2M SNPs), we used the first M/2 SNPs from each of the first four chromosomes. The GRM was calculated with SNPs on chromosomes 3 and 4, and SNPs on chromosomes 1 and 2 were treated as the candidate SNPs. The simulated phenotypes were generated from chromosomes 1 and 3, where 1% of the SNPs were randomly selected as being causal. Results are reported for causal SNPs on chromosome 1 and for null SNPs on chromosome 2, the latter of which were not used for building the GRM.

Results for simulations including 1,000 and 10,000 SNPs (M) and sample sizes fixed at 500 case and 500 control subjects are displayed in Table 4 and Tables S14 and S15. Once again, the LTMLM statistic outperformed ATT and MLM as case-control ascertainment became more severe. (A limitation of these simulations is that performing case-control ascertainment on a fixed set of individuals limits case-control sample size; thus, these simulations were restricted to a disease prevalence of 10% or higher. It is reasonable to infer that for rarer diseases with more extreme case-control ascertainment the LTMLM statistic would achieve even higher power gains, as was demonstrated in simulations with simulated genotypes.)

The h² parameter estimates for simulations with real genotypes are displayed in Table 5. The H-E regression estimates are unbiased, but the REML estimates are again downwardly biased at lower prevalence and large N/M.

WTCCC2 MS Dataset

We analyzed the WTCCC2 genotypes together with MS case-control phenotypes: 5,172 MS case subjects and 5,172 control subjects were genotyped on Illumina chips^11,28(see Material and Methods). We compared ATT, ATT with five PCs (ATT + PCs),²² LogR, LogR with five PCs (LogR + PCs), MLM, and LTMLM. We evaluated calibration by using the average χ² and λ_GC over all SNPs; we note that the average χ² and λ_GC are expected to be greater than 1 as a result of polygenic effects.^11,30 We believe that LTMLM is effective in correcting for confounding and that a higher value of λ_GC for LTMLM than for MLM is most likely due to true polygenic signal and reflects the higher power of LTMLM.

We evaluated power by using the average χ² over the 75 published SNPs (Table 6) and the proportion of published SNPs that were significant at various p value thresholds (Table S16). The LTMLM method performed best: it had a 4.3% improvement in average χ² statistics scaled by λ_GC over MLM (jackknife p = 0.005; see Material and Methods) and an even larger improvement over ATT and ATT + PCs, consistent with simulations (Table 2 and Table S3). LTMLM also detected 56/75 known associations as nominally significant (p < 0.05) after λ_GC correction, whereas MLM detected only 53/75, although this difference is not statistically significant. Similar results were obtained when association statistics were calibrated via LD-score regression³¹ (Table S17). A perfectly matched dataset with 4,094 MS case subjects and 4,094 control subjects yielded a similar improvement for LTMLM over MLM (Table S18). We also applied LTMLM to the full unmatched dataset of 10,204 MS subjects and 5,429 control subjects, where a severe mismatch in ancestry between case and control subjects was not representative of a typical GWAS. The LOCO estimates of h² demonstrated inflation before we controlled for population structure (Table S19). In this analysis, the H-E regression estimate of h² produced an unrealistic value of 7.3 on the observed scale (corresponding to 2.8 on the liability scale), which is outside the plausible 0–1 range, suggesting severe population stratification or other severe problems with the data. We do not recommend the use of LTMLM on unmatched samples when such severe problems are detected. For completeness, we report the results of running LTMLM, which results in a loss of power (Table S18).