Control for Population Structure and Relatedness for Binary Traits in Genetic Association Studies via Logistic Mixed Models

Logistic Mixed Models and Score Tests

For a single-variant test, we consider the following logistic mixed model:

where ${\pi }_i=P\left(y_i=1|{\mathbf{X}}_i,G_i,b_i\right)$ is the probability of a binary phenotype (e.g., disease status) for subject i, conditional on their covariates, genotype, and random effects b_i, X_i is a 1 × p row vector of covariates for subject i, α is a p × 1 column vector of fixed covariate effects including an intercept, G_i is the genotype of a genetic variant for subject i, and β is the genotype effect. We assume that $\mathbf{b}\sim N\left(0,{\sum }_{k=1}^K{\tau }_k{\mathbf{V}}_k\right)$ is an n × 1 column vector of random effects, where ${\tau }_k$ are the variance component parameters and V_k are known n × n relatedness matrices. When the number of variance components K = 1, V₁ is usually the genetic relationship matrix estimated from a large number of genetic variants. We allow for multiple random effects to account for complex sampling designs, e.g., hierarchical designs, and environmental covariance structure. The binary phenotypes y_i are assumed to be independent conditional on the random effects b.

To perform the score test for the null hypothesis H₀: β = 0, we need to fit the null logistic mixed model, which is the same for all genetic variants, as

where ${\pi }_{i0}=P\left(y_i=1|{\mathbf{X}}_i,b_i\right)$. We fit Equation 1 using the penalized quasi-likelihood (PQL) method.²² Specifically, let $\mathbf{W}=diag\left\{v_{i0}\right\}$, where $v_{i0}={\pi }_{i0}\left(1-{\pi }_{i0}\right)$, $\mathbf{X}={\left(\begin{array}{cccc}{\mathbf{X}}_1^T& {\mathbf{X}}_2^T& \cdots & {\mathbf{X}}_n^T\end{array}\right)}^T$ be an n × p covariate matrix including an intercept, and $\tilde{\mathbf{Y}}$ be the “working vector” with components ${\tilde{Y}}_i={\mathbf{X}}_i\mathbf{\alpha }+b_i+v_{i0}^{-1}\left(y_i-{\pi }_{i0}\right)$. Under the null hypothesis H₀: β = 0, we iteratively fit the working linear mixed model $\tilde{\mathbf{Y}}=\mathbf{X}\mathbf{\alpha }+\mathbf{b}+\mathbf{\epsilon }$ where $\mathbf{\epsilon }\sim N\left(0,{\mathbf{W}}^{-1}\right)$. We use the computationally efficient AI-REML algorithm^{6, 23} to estimate ${\hat{\tau }}_k$. After obtaining the variance component estimates, the estimates of the fixed effects $\hat{\mathbf{\alpha }}$ and random effects $\hat{\mathbf{b}}$ can be updated, followed by the working vector $\tilde{\mathbf{Y}}$. The process continues until convergence.

At each iteration, we perform a matrix inversion based on Cholesky decomposition with complexity of O(n³) and matrix multiplications with complexity of O(pn² + K²n² + p²n), where n is the sample size, p is the number of covariates, and K is the number of variance components. Thus, the total complexity of fitting a logistic mixed model with K known relatedness matrices is O(in³ + ipn² + iK²n² + ip²n), where i is the number of iterations required to reach convergence.

The score for H₀: β = 0 is $T={\mathbf{G}}^T\left(\mathbf{y}-{\hat{\mathbf{\pi }}}_0\right)$, where $\mathbf{G}={\left(\begin{array}{cccc}{G}_1& G_2& \cdots & G_n\end{array}\right)}^T$ is the n × 1 column vector of genotypes, $\mathbf{y}={\left(\begin{array}{cccc}{y}_1& y_2& \cdots & y_n\end{array}\right)}^T$ is the n × 1 column vector of outcomes, and ${\hat{\mathbf{\pi }}}_0$ is a vector of fitted values under H₀, which is the same for all SNPs. The estimated variance of the score is $Var\left(T|H_0\right)={\mathbf{G}}^T\hat{\mathbf{P}}\mathbf{G}$ under the null hypothesis, where $\hat{\mathbf{P}}={\hat{\mathbf{\Sigma }}}^{-1}-{\hat{\mathbf{\Sigma }}}^{-1}\mathbf{X}{\left({\mathbf{X}}^T{\hat{\mathbf{\Sigma }}}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T{\hat{\mathbf{\Sigma }}}^{-1}$, and $\hat{\mathbf{\Sigma }}={\hat{\mathbf{W}}}^{-1}+{\sum }_{k=1}^K{\hat{\tau }}_k{\mathbf{V}}_k$. The test for each genetic variant involves a vector-matrix-vector multiplication and an inner product calculation for two vectors, and thus the score test step has complexity of O(qn² + qn), where q is the total number of genetic variants tested. Also, the score test step can be easily parallelized if necessary. We use the C++ library Armadillo to perform matrix and vector calculations. More details about fitting the logistic mixed model and performing the score test are given in Appendix A.

HCHS/SOL Genotyping and Statistical Analysis

HCHS/SOL subjects who consented to genetic studies had DNA extracted from blood. These DNA samples were genotyped on the SOL HCHS Custom 15041502 B3 array (custom content designed and developed by Papanicolaou, Rotter, and Taylor) at Illumina Microarray Services. This array comprised the Illumina Omni 2.5M array (HumanOmni2.5-8v1-1) and additional custom content selected for HCHS/SOL, including ancestry-informative markers, variants characteristic of Amerindians, and known GWAS hits.²⁵ Illumina Microarray Services, LA Biomed, and the SOL Genetic Analysis Center (GAC) performed quality control to generate recommended SNP- and sample-level quality filters. Samples were checked for annotated versus genetic sex, gross chromosomal anomalies, relatedness and population structure, missing call rates, batch effects, duplicate sample discordance, and Mendelian errors. At the SNP level, checks were performed for Hardy-Weinberg equilibrium, minor allele frequency (MAF), duplicate probe discordance, and missing call rate.

Study participants were recruited via a multi-stage survey sampling procedure, for which individuals were sampled within households that were sampled within block groups.²⁶ The study includes genetic data from 12,803 individuals of Hispanic/Latino origin, belonging to six different Hispanic/Latino background groups.²⁵ The HCHS/SOL study was approved by institutional review boards at participating institutions, and written informed consent was obtained from all participants. Standard quality control, similar to that described by Laurie et al.,²⁷ was performed to filter SNPs and samples of poor quality. Additionally, samples missing information on asthma diagnosis and 56 samples identified as ancestry outliers from principal-component analysis were removed, leaving a sample size of 11,222 for analysis. We also filtered SNPs with a MAF less than 5% in this overall sample, resulting in a set of 1,299,221 autosomal SNPs to be analyzed. The proportions of asthma cases in analyzed participants within each of the genetic background groups of Central American, Cuban, Dominican, Mexican, Puerto Rican, and South American descent are 4.5%, 9.6%, 9.6%, 3.9%, 25.6%, and 4.1%, respectively (Table 1). The Puerto Rican group is a clear outlier from the rest, so we compared it to the collection of all other background groups, which has a combined proportion of 5.8%. To account for potential shared environment effects resulting from the sampling design, relatedness matrices representing household and block group membership were included in addition to the standardized genetic relationship matrix, totaling three random effects for the linear and logistic mixed models. The top five PCs, as well as Hispanic/Latino background group indicators, were used to adjust for ancestry in both models. Additional fixed effects covariates included field center, sex, age, cigarette use, cigarette pack years, and sampling weight (i.e., weights to account for disproportionate selection of the sample according to the sample survey design²⁶). Treatment of sampling weight as a fixed effect in this way can effectively adjust for the marginal effect of design variables.²⁸

Simulation Studies

We perform coalescent simulations²⁹ to generate genotypes for a total of 8,000 founders with 1,000,000 independent SNPs from a 20 × 20 grid (Figure S1) to mimic spatially continuous populations (20 individuals per cell) with migration rate between adjacent cells M = 10 to represent population structure within Europe.^{30, 31} To simulate genotypes of an offspring cohort, we first sampled (without replacement) 10 pairs of parents for each cell of the 20 × 20 grid from the original cohort and then simulated two children for each family using the gene dropping algorithm,³² thus obtaining 8,000 individuals in the offspring cohort. We combined the two cohorts to get 16,000 individuals with both population structure and relatedness. For individual j in family i, the probability of being a case ${\pi }_{ij}$ was calculated from

where Z_i = 1 if family i was from a 10 × 10 grid in the top left, and Z_i = 0 otherwise. The parameters ${\alpha }_0$ and ${\alpha }_1$ were chosen such that the disease prevalence was 0.28 in the high-risk population group in the top left given no random effects and 0.05 in the low-risk population group given no random effects. The random effects were simulated as

where the variance component parameter ${\sigma }^2$ was set to 2. We randomly sampled 10,000 individuals from the combined original and offspring cohorts to form a simulated cohort study, and their standardized genetic relationship matrix was calculated using 625,583 genetic variants with MAF greater than 5% in the founder population. We removed family indicators and compared linear and logistic mixed models using the genetic relationship matrix. We included the top ten ancestry PCs as covariates in both models. We analyzed common genetic variants with MAF greater than 5%. In the simulated case-control study for rare disease, we chose ${\alpha }_0$ and ${\alpha }_1$ such that the disease prevalence was 0.045 in the high-risk population group in the top left given no random effects and 0.005 in the low-risk population group given no random effects. We randomly sampled 1,667 case subjects and 8,333 control subjects from the combined original and offspring cohorts to form a case-control study with a total sample size of 10,000 and a case-control ratio of 1:5. We performed the same analysis as described above.

Analysis of HCHS/SOL Asthma GWAS Data

We compared the results from LMMs and logistic mixed models using GMMAT for an analysis of physician-diagnosed asthma in the population-based HCHS/SOL cohort study.

Figure 1A shows the behavior of the overall quantile-quantile (QQ) plots from each method in the asthma analysis. The inflated results from the LMM are caused by violation of its constant residual variance (homoscedasticity) assumption. For binary traits, population stratification affects both population-specific means (disease prevalences in cohort studies) and variances of the trait; population groups with disease prevalence closer to 0.5 have larger variances (Table 1 and Figure 2). The mean-variance relationship assumed in linear models and LMMs is misspecified for binary traits. Although confounding by population structure can be accounted for by adjusting for population groups or ancestry PCs, unequal binary trait variances across different populations are not appropriately modeled in linear models and LMMs. The variance of the standard LMM-based test statistic for no genetic association is underestimated when population groups with larger binary trait variances (i.e., higher disease risk) also have higher MAFs, thus larger genotypic variances. This leads to inflation in the test statistic, and vice versa.

To demonstrate this, we partitioned all common SNPs across the genome into three categories based on their genotypic variances in Puerto Ricans and non-Puerto Ricans: (1) the SNPs with the ratio of expected variances, that is Var(SNP) = $2MAF\left(1-MAF\right)$, in Puerto Ricans (high-risk group) over non-Puerto Ricans (low-risk group) less than 0.8; (2) the SNPs with the ratio of Var(SNP) in Puerto Ricans over non-Puerto Ricans between 0.8 and 1.25; and (3) the SNPs with the ratio of Var(SNP) in Puerto Ricans over non-Puerto Ricans greater than 1.25. In this classification, the category 1 SNPs (n = 144,815 [11%]) have appreciably lower MAFs in Puerto Ricans (high risk) than non-Puerto Ricans (low risk). The category 2 SNPs have similar MAFs in both groups (n = 982,805 [76%]). The category 3 SNPs (n = 171,601 [13%]) have appreciably higher MAFs in Puerto Ricans (high risk) than non-Puerto Ricans (low risk). Both LMMs and GMMAT perform well for category 2 SNPs (Figure 1C). The problem with LMMs is more apparent for category 1 and 3 SNPs: Figure 1B shows conservative p values for category 1 SNPs, and Figure 1D shows anti-conservative p values for category 3 SNPs, as expected. This indicates the p values calculated from LMMs are likely to be incorrect for at least 24% of SNPs in HCHS/SOL asthma GWASs. Meta-analysis by analyzing Puerto Ricans and non-Puerto Ricans separately after removing related individuals across the two groups improves the LMM performance, but the results are still not well calibrated for category 1 and 3 SNPs (Figures 1B and 1D), probably due to the heteroscedasticity issue caused by within-group population stratification. In contrast, GMMAT performs well for all categories of SNPs considered.

Simulation Studies

We also conducted extensive simulation studies under the null hypothesis to illustrate this issue in various GWAS case-control and cohort study designs, to exclude the possibility that LMM-inflated results in HCHS/SOL were caused by polygenic effects. We restricted analyses to common SNPs with a MAF greater than 5%. First, in a simulated cohort study of 10,000 individuals with cryptic relatedness, we simulated two population groups with disease prevalences of 28% (population 1) and 5% (population 2), respectively, from a map of spatially continuous populations. The disease prevalences were chosen to mimic the asthma disease proportions of Puerto Ricans and non-Puerto Ricans in HCHS/SOL (Table 1). In both LMM and GMMAT, we adjusted for the top ten ancestry PCs as fixed effects covariates. Figure 3A includes results from 3,200 null simulated datasets with 625,583 SNPs and 10,000 related subjects per dataset. LMM has a type I error rate of 1.26 × 10⁻⁷ at the nominal 5 × 10⁻⁸ level, compared to 5.0 × 10⁻⁸ for GMMAT. Note that because more than 2 billion p values are included, type I error rate estimates at this level are accurate with standard error 5 × 10⁻⁹.

Following the HCHS/SOL example, we partitioned all the common SNPs into three categories: (1) the SNPs with the ratio of expected variances in population 1 (high risk) over population 2 (low risk) less than 0.8; (2) the SNPs with the ratio of expected variances in population 1 over population 2 between 0.8 and 1.25; and (3) the SNPs with the ratio of expected variances in population 1 over population 2 greater than 1.25. Category 1 includes the SNPs that have appreciably lower MAFs in population 1 than population 2 (24%). Category 2 includes the SNPs that have comparable MAFs in the two populations (58%). Category 3 includes the SNPs that have appreciably higher MAFs in population 1 than population 2 (18%). For category 2 SNPs, both LMM and GMMAT properly control for type I error rates (Figure 3C). However, despite adjusting for ancestry PCs, LMMs have deflated type I error rates for category 1 SNPs (Figure 3B) and inflated type I error rates for category 3 SNPs (Figure 3D). LMMs often fail to control type I error rates in the presence of moderate or strong population stratification when differences in prevalence by population groups cause large differences in binary trait variances between populations, as demonstrated by our simulation studies that were motivated by the HCHS/SOL study. However, this issue is not always evident in the overall QQ plot using all the SNPs in a GWAS, because the number of SNPs in a typical GWAS is usually in the range of 10⁵ and 10⁶, confidence intervals are relatively wide, and inflation and deflation are likely to mask each other in the overall QQ plot.

To illustrate this, Figure S2 shows results from only one simulated dataset with 625,583 common SNPs and 10,000 subjects in the same setting as Figure 3. There is a slight indication of inflation for p values from LMMs in Figure S2A, but the results are conservative for category 1 SNPs (n = 151,206 [24%]) (Figure S2B) and anti-conservative for category 3 SNPs (Figure S2D) (n = 111,455 [18%]). These results show that the p values from LMMs might be incorrect for at least 42% of SNPs in this analysis, even though inflation and deflation might not be obvious in the overall QQ plot.

We next simulated a case-control study with a total sample size of 10,000 and a case-control ratio of 1:5 from the same two populations as in our simulated cohort study in the presence of cryptic relatedness. Compared to cohort studies, case-control studies oversample cases. The disease prevalence was 4.5% and 0.5% in population 1 (high risk) and population 2 (low risk), respectively. Figure S3A shows that GMMAT works well according to the QQ plot of p values. The genomic control inflation factors show that both LMM and GMMAT have appropriate median p values. However, type I error rates at the nominal level of 5 × 10⁻⁸ are 1.75 × 10⁻⁷ for LMM and 4.6 × 10⁻⁸ for GMMAT, respectively. Moreover, when we divided the SNPs into three categories as we previously described, we observed strong deflation for category 1 SNPs (Figure S3B) and strong inflation for category 3 SNPs (Figure S3D) from LMM analysis, as seen in the simulated cohort study. Additional simulation studies show that the same problem for LMMs also exists in cohort studies without cryptic relatedness in the presence of population stratification (Figure S4). Although confounding by discrete populations can be appropriately accounted for by analyzing each population separately via LMMs followed by meta-analysis (Figure S5), this strategy does not perform well in the presence of confounding by continuous population structure (Figure S6). Moreover, LMMs do not work well in case-control studies of unrelated individuals with moderate to strong population stratification (Figure S7). These suggest that despite its wide use in both population-based and family-based genetic association studies, LMMs are generally not appropriate for binary traits due to misspecified phenotype variance in the model, probably yielding incorrect p values (see Appendix B for details).

We also performed additional simulations to compare GMMAT with ROADTRIPS, which performs a retrospective test for association in case-control data,³³ and the recently developed liability estimator as a phenotype (LEAP) approach. LEAP fits a liability mixed model that accounts for case-control ascertainment.²⁰ In the same case-control setting as that described for Figure S7, ROADTRIPS was not well calibrated due to failure to fully account for population stratification. LEAP was found to have well-behaved overall QQ plots but inflated type I error rates for category 1 SNPs, regardless of whether top ten ancestry PCs were adjusted as covariates or not (Figure S8).

We next simulated case-control studies with the disease prevalence of 1% by varying the case-control sampling ratios in two population groups. We first considered balanced cases and controls in the two populations. Specifically, when the case-control ratio is 1:1 in both groups, there is no population stratification, and ROADTRIPS, LEAP, LMMs, and GMMAT all properly control for type I error rates (Figure S9). When the case-control ratios are different but flipped in two populations (4:1 in population 1 and 1:4 in population 2), there is population stratification. However, because the variance of the binary trait is the same in both population groups in this situation, LMMs and GMMAT both perform well, whereas ROADTRIPS has inflated type I error rates in the tail and LEAP has an appropriate overall QQ plots but inflated or deflated type I error rates for the SNPs whose MAFs are different in the two population groups (Figure S10).

We also considered a situation where the case-control ratios are different in two populations (25:2 in population 1 and 25:48 in population 2) in a way that led to a smaller variance of the binary trait in population 1 than in population 2. ROADTRIPS and LEAP adjusting for the top ten PCs show inflation in the overall QQ plot, and LEAP without covariates and LMMs have inflated type I error rates for category 1 SNPs and deflated type I error rates for category 3 SNPs (Figure S11). In terms of required computational resources, LEAP requires more than 20 times the amount of memory compared to GMMAT for an analysis of a case-control study of sample size 10,000.

We also performed simulation studies to evaluate the performance of GMMAT in estimating odds ratios of genetic variants. We found that GMMAT had minor bias in estimating odds ratios when fitting logistic mixed models under the null and alternative hypotheses, compared to the true parameter values in large samples (n = 10,000). In addition, the performance in terms of absolute bias appears to be similar to that using MACAU,³⁴ a recently developed Markov Chain Monte Carlo-based mixed model approach for binomial count data, while the GMMAT odds ratio estimates are less variable (Figure S12).

In the absence of population stratification, LMMs and GMMAT have comparable power, but they are both less powerful than logistic regression (Table S1). In the presence of population stratification, LMMs have less power than GMMAT for causal genetic variants with a lower MAF in the high-risk population than in the low-risk population (Table S1), due to its conservative type I error rate control for such variants, as shown in Figure S4B.

Computational Speed and Memory Usage

We benchmarked our GMMAT package against SAS PROC GLIMMIX regarding computational speed and memory usage. To fit a null model with sample size 2,000, GMMAT takes less than 1.5% of the time required by SAS PROC GLIMMIX when fitting a logistic mixed model with one variance component, and less than 0.6% of the time with three variance components (Table S2), yielding the same numerical results to at least the fourth decimal place. On average, with one variance component of random effects, as is commonly used to account for genetic relatedness in GWASs, SAS PROC GLIMMIX requires about 28 min to fit a null model on a single core of an Intel Xeon E5-2690 CPU (2.90 GHz), compared to about 22 s using GMMAT. With three variance components of random effects that account for complex sampling designs in addition to genetic relatedness, SAS PROC GLIMMIX requires about 1.2 hr, whereas GMMAT needs still about 22 s. It takes about 14 min for GMMAT to perform score tests for 1,000,000 genetic variants on the same core, without parallelization. In practice, score tests for different variants can be easily run in parallel in a computing cluster, and testing each genetic variant takes only about 0.8 ms.

GMMAT requires less than 1 GB memory in analyzing 2,000 individuals. With sample size 10,000, it takes about 18 min to fit the null model with one variance component and 34 min with three variance components using GMMAT, and about 3.6 hr to perform score tests for 1,000,000 genetic variants on a single core (about 13 ms for testing each genetic variant) (Table S2). SAS PROC GLIMMIX reports insufficient memory when 100 GB memory is specified to fit the null model for 10,000 individuals, whereas GMMAT requires less than 14 GB memory for one variance component and 21 GB for three variance components.

The Generalized Linear Mixed Model

The derivations below are based on generalized linear mixed models (GLMMs), and logistic mixed models are a special case of GLMMs when the link function is logit and the dispersion parameter is fixed at 1. In the context of single-variant test, we consider the following GLMM

where X_i is a 1 × p row vector of covariates for subject i, $\mathbf{\alpha }$ is a p × 1 column vector of fixed covariate effects including the intercept, G_i is the genotype of the genetic variant of interest for subject i, and β is the fixed genotype effect. We assume that $\mathbf{b}\sim N\left(0,{\sum }_{k=1}^K{\tau }_k{\mathbf{V}}_k\right)$ is an n × 1 column vector of random effects, ${\tau }_k$ are the variance component parameters, $\mathbf{\tau }$ is a K × 1 column vector of ${\tau }_k$, and V_k are known n × n matrices. We also assume that given the random effects b, the outcome y_i is conditionally independent with mean $E\left(y_i|\mathbf{b}\right)={\mu }_i$ and variance $Var\left(y_i|\mathbf{b}\right)=\varphi a_i^{-1}v\left({\mu }_i\right)$, where $\varphi$ is the dispersion parameter (for binary and Poisson data $\varphi =1$), a_i are known weights, and $v(\cdot )$ is the variance function. The linear predictor ${\eta }_i$ is a monotonous function of the conditional mean ${\mu }_i$ via the link function ${\eta }_i=g\left({\mu }_i\right)$. For binary traits y_i, ${\mu }_i={\pi }_i=P\left(y_i=1|{\mathbf{X}}_i,G_i,b_i\right)$ is the probability of the binary outcome (e.g., disease status) for subject i.

For subject i, the quasi-likelihood given random effects b is

The log integrated quasi-likelihood function of $\left(\mathbf{\alpha },\beta ,\varphi ,\mathbf{\tau }\right)$ is

Let

we can use Laplace method to approximate the n-dimensional integral

thus Equation A1 becomes

where

is the solution of $f^{\prime }\left(\mathbf{b}\right)=0$.

The first partial derivative of $q{l}_i\left(\mathbf{\alpha },\beta ;\mathbf{b}\right)$ with respect to b is

where Z_i is a 1 × n vector of indicators such that $b_i={\mathbf{Z}}_i\mathbf{b}$, ${\mathbf{I}}_n=\left(\begin{array}{cccc}{\mathbf{Z}}_1^T& {\mathbf{Z}}_2^T& \cdots & {\mathbf{Z}}_n^T\end{array}\right)$, and the second derivative is

For canonical link functions, $v\left({\mu }_i\right)g^{\prime }\left({\mu }_i\right)=1$, the last two terms become 0. Let

then Equation A2 becomes

We assume that the weight matrix W changes slowly with respect to the conditional means (following Breslow and Clayton²²), that is

then we take the derivatives of Equation A3:

Under the null hypothesis $H_0:\beta =0$, if $\varphi$ and $\mathbf{\tau }$ are known, we jointly choose $\hat{\mathbf{\alpha }}\left(\varphi ,\mathbf{\tau }\right)$ and $\hat{\mathbf{b}}\left(\varphi ,\mathbf{\tau }\right)$ to maximize Equation A3, then $\hat{\mathbf{b}}\left(\varphi ,\mathbf{\tau }\right)=\tilde{\mathbf{b}}\left(\hat{\alpha }\left(\varphi ,\mathbf{\tau }\right),\beta =0\right)$ because $\tilde{\mathbf{b}}$ maximizes f(b) for given $\left(\mathbf{\alpha },\beta \right)$. Defining the working vector $\tilde{\mathbf{Y}}$ with elements ${\tilde{Y}}_i={\eta }_i+g^{\prime }\left({\mu }_i\right)\left(y_i-{\mu }_i\right)$, the solution of

can be written as the solution to the system

Let $\mathbf{\Sigma }={\mathbf{W}}^{-1}+{\sum }_{k=1}^K{\tau }_k{\mathbf{V}}_k$, $\mathbf{P}={\mathbf{\Sigma }}^{-1}-{\mathbf{\Sigma }}^{-1}\mathbf{X}{\left({\mathbf{X}}^T{\mathbf{\Sigma }}^{-1}\mathbf{X}\right)}^{-1}{\mathbf{X}}^T{\mathbf{\Sigma }}^{-1}$, then

is the solution that maximizes Equation A3. We note that

Estimation of Variance Component Parameters

Following Breslow and Clayton,²² we ignore the dependence of W on $\mathbf{\tau }$ and use Pearson chi-square statistic to approximate the deviance

Then Equation A3 at the maximum becomes

Similarly, the restricted maximum likelihood (REML) version is

Let ${\mathbf{V}}_0=diag\left\{a_i^{-1}v\left({\mu }_i\right){\left[g^{\prime }\left({\mu }_i\right)\right]}^2\right\}={\varphi }^{-1}{\mathbf{W}}^{-1}$, then $\mathbf{\Sigma }=\varphi {\mathbf{V}}_0+{\sum }_{k=1}^K{\tau }_k{\mathbf{V}}_k$, and the first derivatives of $q{l}_R\left(\hat{\mathbf{\alpha }}\left(\varphi ,\mathbf{\tau }\right),\beta =0,\varphi ,\mathbf{\tau }\right)$ with respect to $\varphi$ and $\mathbf{\tau }$ are

We define the average information⁶, ²³ matrix AI with the following entries

Let $\mathbf{\theta }$ be the variance component parameters to estimate, when $\varphi \ne 1$, $\mathbf{\theta }=\left(\varphi ,\mathbf{\tau }\right)$, and AI is a (K + 1)×(K + 1) matrix. For binary and Poisson data, $\varphi =1$, $\mathbf{\theta }=\mathbf{\tau }$, and AI is a K × K matrix containing only $\mathbf{A}{\mathbf{I}}_{{\tau }_k{\tau }_l}$.

We use the following algorithm to fit the null GLMM:

• 1.

  Fit a generalized linear model with $\mathbf{\tau }=0$ and get ${\hat{\mathbf{\alpha }}}^{\left(0\right)}$ and working vector ${\tilde{\mathbf{Y}}}^{\left(0\right)}$;

• 2.

  Use ${\mathbf{\theta }}^{\left(0\right)}=Var\left({\tilde{\mathbf{Y}}}^{\left(0\right)}\right)/K$ (if $\varphi =1$) or ${\mathbf{\theta }}^{\left(0\right)}=Var\left({\tilde{\mathbf{Y}}}^{\left(0\right)}\right)/\left(K+1\right)$ (if $\varphi \ne 1$) as the initial value of $\mathbf{\theta }$;

• 3.

  For each $k=\mathrm{0,1},\dots ,K$, update $\mathbf{\theta }$ using ${\theta }_k^{\left(1\right)}={\theta }_k^{\left(0\right)}+2n^{-1}{\left\{{\theta }_k^{\left(0\right)}\right\}}^2\left(\partial q{l}_R\left({\mathbf{\theta }}^{\left(0\right)}\right)/\partial {\theta }_k\right)$;

• 4.

  Use ${\tilde{\mathbf{Y}}}^{\left(1\right)}={\tilde{\mathbf{Y}}}^{\left(0\right)}$ as $\tilde{\mathbf{Y}}$ and update ${\mathbf{\theta }}^{\left(2\right)}={\mathbf{\theta }}^{\left(1\right)}+{\left\{\mathbf{A}{\mathbf{I}}^{\left(1\right)}\right\}}^{-1}\left(\partial q{l}_R\left({\mathbf{\theta }}^{\left(1\right)}\right)/\partial \mathbf{\theta }\right)$;

• 5.

  Calculate ${\hat{\mathbf{\alpha }}}^{\left(2\right)}$ and ${\hat{\mathbf{b}}}^{\left(2\right)}$ using ${\tilde{\mathbf{Y}}}^{\left(1\right)}$ and ${\mathbf{\theta }}^{\left(2\right)}$;

• 6.

  Update ${\tilde{\mathbf{Y}}}^{\left(2\right)}$ using ${\hat{\mathbf{\alpha }}}^{\left(2\right)}$ and ${\hat{\mathbf{b}}}^{\left(2\right)}$;

• 7.

  Repeat steps 4–6, until $2\max \left\{\left|{\hat{\mathbf{\alpha }}}^{\left(i\right)}-{\hat{\mathbf{\alpha }}}^{(i-1)}\right|/(\left|{\hat{\mathbf{\alpha }}}^{\left(i\right)}\right|+\left|{\hat{\mathbf{\alpha }}}^{(i-1)}\right|),\left|{\hat{\mathbf{\theta }}}^{\left(i\right)}-{\hat{\mathbf{\theta }}}^{(i-1)}\right|/(\left|{\hat{\mathbf{\theta }}}^{\left(i\right)}\right|+\left|{\hat{\mathbf{\theta }}}^{(i-1)}\right|)\right\}\le$ tolerance.

The Score Test

Once $\left(\hat{\mathbf{\alpha }},\hat{\varphi },\hat{\mathbf{\tau }}\right)$ is estimated under the null hypothesis $H_0:\beta =0$, the score test can be constructed by evaluating Equation A4 at $\left(\hat{\mathbf{\alpha }},\beta =0,\hat{\varphi },\hat{\mathbf{\tau }}\right)$, that is

Its variance under the null hypothesis is

the last equality holds because $\hat{\mathbf{P}}\hat{\mathbf{\Sigma }}\hat{\mathbf{P}}=\hat{\mathbf{P}}$.

Unrelated Individuals with Population Stratification

We performed additional simulation studies to compare LMM and GMMAT in unrelated individuals in the presence of population stratification. We used the coalescent model²⁹ to simulate genotypes for a total of 16,000 unrelated individuals with 1,000,000 genetic variants from a 20 × 20 grid (Figure S1) of spatially continuous populations (40 individuals per cell) with migration rate between adjacent cells M = 10 to mimic population structure within Europe.³⁰, ³¹ For individual i, the probability of being a case ${\pi }_i$ was calculated from

where Z_i = 1 if individual i was from a 10 × 10 grid in the top left (population 1) and Z_i = 0 otherwise (population 2). The parameters ${\alpha }_0$ and ${\alpha }_1$ were chosen such that the disease prevalence was 0.28 in population 1 (high risk) and 0.05 in population 2 (low risk). Note that the mean model is not mis-specified when a linear link function is used, thus our simulation setting is not in favor of logistic models. We randomly sampled 10,000 individuals to form a simulated cohort study for common disease (Figure S4) and calculated their standardized genetic relationship matrix using 625,504 genetic variants with MAF greater than 5%. We adjusted for top ten PCs in both models and combined results from 3,200 null simulation replicates. We compared our method with the strategy in which population 1 and population 2 were analyzed separately using LMMs and the results were subsequently combined in a meta-analysis (which we referred to as LMM meta), using 100 null simulation replicates (Figure S5). Moreover, we considered continuous population stratification where Z_i for each cell is the minimum of its row and column coordinates in Figure S1, which ranges from 0 to 19, and ${\alpha }_0$ was chosen such that the disease prevalence in populations with Z_i = 0 was 0.02 and ${\alpha }_1=0.2$ (Figure S6).

We also simulated a case-control study for rare disease (Figure S7). We chose ${\alpha }_0$ and ${\alpha }_1$ such that the disease prevalence was 0.045 in population 1 (high risk) and 0.005 in population 2 (low risk). We randomly sampled 1,667 case subjects and 8,333 control subjects to form a case-control study with a total sample size of 10,000 and a case-control ratio of 1:5. We performed the same analysis as described above and combined results from 3,200 null simulation replicates.

Comparison with the Existing Methods

We also compared GMMAT with ROADTRIPS³³ and the recently developed liability (probit) mixed models.²⁰, ²¹ ROADTRIPS does not allow for covariates. Because LTMLM²¹ does not allow covariate adjustment, we used LEAP²⁰ and set the disease prevalence to 0.015, which is the pooled prevalence in both high- and low-risk population groups. We compared ROADTRIPS and LEAP without covariates as well as LEAP adjusting for top ten PCs as covariates with LMM and GMMAT. The simulation settings were the same as in Figure S7, but p values from only one simulation replicate were shown (Figure S8).

Moreover, we simulated three case-control settings with the same prevalence but different case-control sampling schemes in two population groups in Figure S1. The disease prevalence was set to 0.01 in both groups. We compared ROADTRIPS and LEAP with and without adjusting for top ten PCs, with LMM and GMMAT.

In the first setting, we randomly sampled 1,250 case subjects and 1,250 control subjects from population 1 and 3,750 case subjects and 3,750 control subjects from population 2. This was a balanced case-control study with balanced designs in both population groups. Because the case-control ratio was the same, there was no population stratification (Figure S9).

In the second setting, we randomly sampled 2,000 case subjects and 500 control subjects from population 1 and 1,500 case subjects and 6,000 control subjects from population 2. This was an unbalanced case-control study with unbalanced designs but equal binary trait variances in two population groups. Population stratification was created by different case-control ratios in two population groups (4:1 in population 1 and 1:4 in population 2), instead of different disease prevalence (Figure S10).

In the third setting, we randomly sampled 2,500 case subjects and 200 control subjects from population 1 and 2,500 case subjects and 4,800 control subjects from population 2. This was a balanced case-control study with unbalanced designs and unequal binary trait variances in two population groups. Because the disease prevalence was 0.01 in both groups, there were no high-risk or low-risk groups, but population 1 was the low binary trait variance group with variance 0.0686 and population 2 was the high binary trait variance group with variance 0.2252 (Figure S11).

Simulations with Genetic Effects

We first conducted simulation studies to evaluate the performance of GMMAT for estimating the odds ratios of genetic variants in the presence of population stratification. We used the same genotype data as in our null simulations of unrelated individuals with population stratification. For individual i, the probability of being a case ${\pi }_i$ was calculated from

where Z_i = 1 if individual i was from population 1, Z_i = 0 otherwise (population 2). G_i is the additively coded genotype for the causal genetic variant. The parameters ${\alpha }_0$ and ${\alpha }_1$ were chosen such that the disease prevalence was 0.28 in population 1 and 0.05 in population 2 for individuals with G_i = 0. β was chosen such that the odds ratio varied from 1.0, 1.1, 1.25, 1.5, to 2.0. The total sample size was 10,000. We compared the odds ratio estimates and the p values calculated by GMMAT and MACAU (Figure S12). We then compared the power of LMM and GMMAT for identifying causal variants that had lower MAFs in population 1 than population 2 (Table S1).

We also simulated a case-control study with 1,667 case subjects and 8,333 control subjects with no population stratification. We assumed the disease prevalence was 0.01 for individuals with G_i = 0 and used the same sampling scheme in both population 1 and population 2. We compared the powers of logistic regression, LMM, and GMMAT for identifying causal variants (Table S1).