Testing a Large Number of Composite Null Hypotheses Using Conditionally Symmetric Multidimensional Gaussian Mixtures in Genome-Wide Studies

In a genetic mediation study, suppose for each of $i=1,...,n$ subjects we observe a single scalar genotype $G_{i}$, an outcome $Y_{i}$, the $J$ mediators $M_{ij}$ $(j=1,...,J)$, and a vector of $p$ additional covariates $\mathbf{X}_{i}=(X_{i1},...,X_{ip})^{T}$. In lung cancer examples, $G_{i}$ is a SNP of interest, $Y_{i}$ is the binary lung cancer status, and the $M_{ij}$ are lung tissue gene expression values for $J$ genes across the genome. We want to perform mediation analysis for each mediator separately by fitting the $J$ pairs of regression models (Dai et al., 2022; Liu et al., 2022)

where $\mu_{ij}=E(Y_{i}|\mathbf{X}_{i},M_{ij},G_{i})$. Let $Z_{j1}$ be the result of a score test for $H_{0,j,1}:\alpha_{j}=0$, and let $Z_{j2}$ be the result of a score test for $H_{0,j,2}:\beta_{j}=0$. It can be shown that $Z_{j1}$ and $Z_{j2}$ are independent (Liu et al., 2022). Under some common assumptions, the SNP $G_{i}$ possesses an effect on lung cancer mediated through expression of gene $j$ if we reject the composite null $H^{med}_{0,j}:\alpha_{j}=0\cup\beta_{j}=0$ (Huang, 2019). This basic two-regression framework can be adapted in many ways, for instance Equations (1) and (2) could be applied in separate datasets if the data are not all available in one package. However, it is clear to see that the form of the composite null $H^{med}_{0,j}$ would be the same.

In a pleiotropy study, suppose that we have two GWAS datasets A and B. For instance, A may be a GWAS for lung cancer in the International Lung Cancer Consortium (ILCCO) dataset (McKay et al., 2017), and B may be a GWAS for coronary artery disease (CAD) in the UK Biobank (UKB) dataset (Bycroft et al., 2018). In cohort A, suppose that for each of $i=1,...,n_{A}$ subjects we observe genotype data $\mathbf{G}_{i}^{A}=(G_{i1}^{A},...,G_{iJ}^{A})^{T}$ for $J$ total SNPs as well as the outcome $Y_{i}^{A}$ and a vector of $p$ additional covariates $\mathbf{X}_{i}^{A}=(X_{i1}^{A},...,X_{ip}^{A})^{T}$. Let $\mathbf{G}_{i^{\prime}}^{B}$, $Y_{i^{\prime}}^{B}$, and $\mathbf{X}_{i^{\prime}}^{B}$ be defined similarly for the $i^{\prime}=1,...,n_{B}$ subjects in cohort B, where the $j$th genotype always indexes the same SNP in both datasets. A pleiotropy analysis fits

with $\mu_{ij}^{A}=E(Y_{i}^{A}|\mathbf{X}_{i}^{A},G_{ij}^{A})$ and $\mu_{i^{\prime}j}^{B}=E(Y_{i^{\prime}}^{B}|\mathbf{X}_{i^{\prime}}^{B},G_{i^{\prime}j}^{B})$. Let $Z_{j1}$ be the result of a score test for $H_{0,j,1}:\alpha_{j}=0$ and $Z_{j2}$ the result of a score test for $H_{0,j,2}:\beta_{j}=0$. SNP $j$ is a pleiotropic variant associated with both lung cancer and CAD if we reject the composite null $H^{pleio}_{0,j}:\alpha_{j}=0\cup\beta_{j}=0$. In massive genetic compendiums such as the UKB, A and B may also be the same dataset. In such cases, $Z_{j1}$ and $Z_{j2}$ would generally be correlated.

A replication study is also performed with Equations (3) and (4) except the outcomes $Y_{i}^{A}$ and $Y_{i^{\prime}}^{B}$ are lung cancer in both models. However, in a replication study we only want to reject the null when the parameters possess the same sign, which creates an even larger null space than the mediation or pleiotropy composite null. We can formalize the replication null as $H^{rep}_{0,j}:\alpha_{j}=0\cup\beta_{j}=0\cup\text{sign}(\alpha_{j})\neq\text{sign}(\beta_{j})$. Additionally, replication and pleiotropy studies may be conducted in more than two datasets. For instance, we may possess a third lung cancer cohort C with outcomes $Y_{i^{\prime\prime}}^{C}$ for $i^{\prime\prime}=1,...,n_{C}$ subjects and data on the same SNPs $\mathbf{G}_{i}^{C}=(G_{i^{\prime\prime}1}^{C},...,G_{i^{\prime\prime}J}^{C})^{T}$ and other covariates $\mathbf{X}_{i^{\prime\prime}}^{C}=(X_{i^{\prime\prime}1}^{C},...,X_{i^{\prime\prime}p}^{C})^{T}$. We can then fit $\mbox{logit}(\mu_{i^{\prime\prime}}^{C})=\gamma_{j0}+\boldsymbol{\gamma}_{jX}^{T}\mathbf{X}_{i^{\prime\prime}}^{C}+\gamma_{j}G_{i^{\prime\prime}j}^{C}$ and let $Z_{j3}$ be the result of a score test for $H_{0,j,3}:\gamma_{j}=0$.

For each of the mediation, pleiotropy, and replication settings, computing all the test statistics results in $J$ sets of the form $\mathbf{Z}_{j}=(Z_{j1},...,Z_{jK})^{T}$, where each element $Z_{jk}\sim N(0,1)$ if the $k$th null in the $j$th set is true, $k=1,...,K$. For example, in mediation analysis $K=2$, and in three-way pleiotropy, $K=3$. We refer to $K$ as the dimension of the composite null set. Let the observed value of $\mathbf{Z}_{j}$ be given by $\mathbf{z}_{j}=(z_{j1},...,z_{jK})^{T}$. These $\mathbf{z}_{j}$ are also called summary statistics and are often publicly available, e.g. in the GWAS catalog.

We briefly review the multidimensional two-group model to help introduce the csmGmm. Suppose we are interested in performing large-scale testing on $J$ sets of $K$ parameters $\boldsymbol{\theta}_{j}=(\theta_{j1},....,\theta_{jK})^{T}$, $j=1,...,J$. Let the general composite null for each set be $H_{0,j}^{comp,K}:\bigcup_{k=1}^{K}\theta_{jk}=0$, and suppose the available data consists of $K$-length test statistic vectors $\mathbf{z}_{1},...,\mathbf{z}_{J}$ that have been computed as described above or provided as summary statistics. We assume that the distribution of $\mathbf{Z}_{j}$ depends on the unobserved true association configuration $\mathbf{H}_{j}=(H_{j1},...,H_{jK})^{T}$ (Efron, 2008; Heller and Yekutieli, 2014). Here $H_{jk}\in\{-1,0,1\}\forall j,k$, where $H_{jk}=0$ denotes that the $k$th individual null in the $j$th set is true, i.e. $\theta_{jk}=0$. Let $H_{jk}=-1$ denote that the $k$th individual parameter in the $j$th set is negative, i.e. $\theta_{jk}<0$, and let $H_{jk}=1$ denote that the $k$th individual parameter in the $j$th set is positive, i.e. $\theta_{jk}>0$. While there are technically three groups, we use the two-group name to follow historical precedence.

Let $\mathbf{h}^{l}=(h^{l}_{1},...,h^{l}_{K})^{T}\in\mathcal{H}$ denote possible values of $\mathbf{H}_{j}$, where $l=0,...,L=3^{K}-1$, and each $l$ uniquely identifies one possible association configuration. We use $\mathbf{h}^{0}=(0,...,0)^{T}$ to denote the global null configuration always. The general composite null space includes all configurations with at least one 0, $\mathcal{H}_{0}=\{\mathbf{h}^{l}\in\mathcal{H}:\sum_{k=1}^{K}|h^{l}_{k}|<K\}$ (replication is a special case that will be discussed later). The local false discovery rate values are then given by $lfdr(\mathbf{z}_{j}):=\text{Pr}(\mathbf{H}_{j}\in\mathcal{H}_{0}|\mathbf{Z}_{j}=\mathbf{z}_{j})$ (Efron, 2008), which is simply

It is clear that the performance of the lfdr approach depends heavily on the specifications of the densities $f(\mathbf{z}_{j}|\mathbf{H}_{j}=\mathbf{h}^{l})$. Many options are available (Efron, 2008; Heller and Yekutieli, 2014; Wang and Wei, 2020) but may not (a) be flexible enough to unify the mediation, pleiotropy, and replication settings, (b) possess good operating characteristics in each application, or (c) provide interpretable results that do not contradict conventional frequentist significance rankings, which is especially important in the z-statistic-dominated field of genetics research. The csmGmm addresses these issues by building upon a Gaussian mixture model, which is known (Chen and Li, 2009) for its ability to approximate densities arbitrarily well, helping to address (a). We then modify the mixture with conditional symmetry constraints that greatly decrease the number of parameters, helping to address (b), and allow us to prove that (c) holds. More specific comparisons with other density models are provided in the Supplementary Materials.

We will first consider the mediation framework with $K=2$ as a representative example to motivate and propose the csmGmm. The general setting of $K\geq 2$ follows.

In a pleiotropy analysis conducted within a single dataset, e.g. the UK Biobank, test statistics in a set may be correlated. As the base csmGmm assumes independence between test statistics in the same set, it would be a poor fit for such correlated situations. Thus, we propose a correlated csmGmm (c-csmGmm) model for sets of dimension $K=2$. Suppose the correlation between the test statistics is known or can be well-estimated as $\rho$. For instance, for two genome-wide association studies, one can simply use the sample correlation of the test statistics (Liu and Lin, 2018) to estimate $\rho$. Then, conditional on the hypothesis configuration $\mathbf{h}^{l}$, the model for $\mathbf{Z}_{j}$ is

Model (9) assumes two-dimensional composite nulls and $M_{1}=M_{2}=M_{3}=1$. Here, similar to Section 3.2.2, we let $\mu_{b,k}=0$ when configurations with binary representation $b$ have a 0 in the $k$th dimension. We can then prove a result about incongruous results similar to Theorem 1 for the c-csmGmm as well.

In other words, the c-csmGmm prevents incongruous results in the tail. This behavior is desirable since the tail is generally the area of interest, where we would expect to find significant discoveries. Incongruous results for test statistics near 0 may possess little impact on the final conclusions of an analysis. The condition that $\rho\geq 0$ is also not restrictive as generally pleiotropy studies are conducted on diseases that are positively correlated with each other. We prove Theorem 2 in the Supplementary Materials. The approach is substantially more complicated than Theorem 1, as the correlation adds complexity to each joint density.

Replication analyses also require slightly different considerations because of the different form of the composite null. In a replication analysis, if we are again testing the parameters $\boldsymbol{\theta}_{j}$ with $\mathbf{z}_{j}$, the general replication null is $H_{0,j}^{rep,K}:\{\bigcup_{k=1}^{K}\theta_{jk}=0\}\cup\{\bigcup_{k\neq k^{\prime}}\mbox{sign}(\theta_{jk})\neq\mbox{sign}(\theta_{jk^{\prime}})\}$. Correspondingly, we need to modify the null and alternative spaces in the two-group model, with $\mathcal{H}^{rep}_{a}=\{(1,1,...,1)^{T}_{K\times 1},(-1,-1,...,-1)^{T}_{K\times 1}\}$ and $\mathcal{H}^{rep}_{0}=\mathcal{H}\backslash\mathcal{H}^{rep}_{a}$.

The elements of $\mathbf{z}_{j}$ will still be mutually independent, and only a slight modification of the base csmGmm is needed to construct the replication csmGmm (r-csmGmm). Specifically, the r-csmGmm is defined as model (7) with the additional restriction that $M_{2^{K}-1}=1$. In other words, conditional on alternative configurations with $b_{l}=2^{K}-1$, the density of each dimension must be a mixture of two Gaussian distributions with the same mean magnitude and opposite directions (see example of Section 3.2.1). This constraint is reasonable when test statistics under the alternative could plausibly exhibit such a bimodal symmetric shape. The assumption may not be appropriate if test statistics under the alternative are believed to possess a density with more than two peaks. We then have the following theorem.

Theorem 3 differs from Theorem 1 in that congruous results are only guaranteed when all observed test statistics possess the same sign. This condition is not restrictive as observed test statistics with different signs generally would not be considered significant in a replication analysis anyway. We prove Theorem 3 in the Supplementary Materials.

The full data likelihood of $\mathbf{z}_{1},...,\mathbf{z}_{J}$ is very complex in genetic association studies due to correlation between SNPs. We can instead consider the composite likelihood ignoring correlations among features $j^{\prime}\neq j$, which has been justified in local dependency settings (Lindsay, 1988; Heller and Yekutieli, 2014). This composite likelihood is

where $\mathbf{C}_{j}=(\mathbf{H}_{j}^{T},D_{j})^{T}$ identifies the assocation configuration and mean vector, $\mathbf{C}_{j}=\mathbf{c}^{l,m}$ denotes that $\mathbf{C}_{j}=((\mathbf{h}^{l})^{T},m)^{T}$, and $\mathcal{C}$ is the collection of all possible mixture components. The log-likelihood is then

where $\boldsymbol{\Lambda}_{l}:=\text{diag}\{\mathbf{h}^{l}\}$, and $\boldsymbol{\Sigma}$ is either the identity matrix or correlation matrix specified in the c-csmGmm, depending on the setting. We can then estimate the unknown parameters $\boldsymbol{\pi}$ and $\boldsymbol{\mu}$ using a standard expectation-maximization (EM) algorithm (Dempster et al., 1977).

Much literature has investigated the convergence of the EM algorithm and other technical issues in mixture distributions (Chen and Li, 2009). The full details of our implementation are provided in the Supplementary Materials, as we note that the focus of this work is performing inference on the $\textbf{z}_{j}$. We do highlight that effect directions may not appear symmetrically distributed around 0 for certain studies. For instance, some lipids summary statistics are oriented such that all effects are positive (Willer et al., 2013). However, using all positive orientations is clearly an artificial decision, and no information is lost by simply switching effect and non-effect alleles, which negates the test statistic. When test statistics appear skewed to one side, simulations show that randomly performing effect-reference switching to achieve rough symmetry about the origin works well (see Supplementary Materials).

To perform inference, we can apply the standard result that the mean lfdr-value in a rejection region is the probability of a false rejection (Efron and Tibshirani, 2002; Heller and Yekutieli, 2014). Specifically, we define the rejection region for a nominal false discovery rate $q$ to be $\mathcal{R}_{q}:=\{\mathbf{z}_{j}:\widehat{lfdr}(\mathbf{z}_{j})\leq\widehat{lfdr}_{(\hat{t}(q))}\}$. Here $\widehat{lfdr}_{(j)}$ denotes the order statistics of the $\widehat{lfdr}(\mathbf{z}_{j})$, $\widehat{lfdr}(\mathbf{z}_{j})$ are the lfdr-values computed under Equation (5) and a csmGmm model with estimates $(\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\pi}})$ in place of $(\boldsymbol{\mu},\boldsymbol{\pi})$, and $\widehat{t}(q)=\max_{r}\sum_{j=1}^{r}\widehat{lfdr}_{(j)}/r\leq q$ is the largest $r$ such that an average of the $r$-smallest lfdr-values is less than or equal to $q$.

Figure 1 considers the original csmGmm in the basic two-dimensional mediation setting. Each simulation iteration generates $j=1,...,10^{5}$ sets of two test statistics by simulating under $M_{ij}=\alpha_{0}+\alpha_{j}G_{ij}+\epsilon_{ij}$ and $\mbox{logit}(\mu_{ij})=\beta_{0}+\beta_{j}M_{ij}+\beta_{G}G_{ij}$. Here $\mbox{logit}(\mu_{ij})=E(Y_{ij}|M_{ij},G_{ij})$, $G_{ij}$ is a Binomial$(2,0.3)$ random variable, $(\alpha_{0},\beta_{0},\beta_{G})=(0,-1,0)$, $\epsilon_{ij}\overset{iid}{\sim}N(0,1)$, and $n=1000$. Standard score statistics for $(\alpha_{j},\beta_{j})$ are calculated by fitting the same regression models. We vary both the proportion of $(\alpha_{j},\beta_{j})$ that are non-zero and their effect sizes; effects are split evenly between positive and negative directions (see Supplementary Materials for all positive case). Figures 1A and 1B consider the composite null setting where there are no true mediated effects, while Figures 1C and 1D include a non-zero proportion of mediated effects (case 3), specifically $\tau_{2}=3\tau_{1}^{2}$. Here $\tau_{1}$ is the common proportion of each configuration with only one true association (cases 1 and 2).

These simulations demonstrate two major points about the mediation setting. First, when $\tau_{2}=0$ and there are no composite null sets generated under the alternative (Figures 1A and 1B), all methods can protect the false discovery rate well, and all are generally conservative. This finding confirms previous reports (Dai et al., 2022; Liu et al., 2022). Second, when some sets are generated under the alternative, all methods can protect the nominal error rate if $\tau_{1}$ (and thus $\tau_{2}$) is small, or if associations throughout the dataset are small. However if effects are generally large (Figure 1C) or $\tau_{1}$ is large (Figure 1D), some approaches may be anticonservative. In Figure 1C we see that both the 7 degrees of freedom spline and the DACT are anticonservative when associations generate test statistics with means in the range of 3 to 4 ($x\approx 0.15$). The DACT was developed for rare-weak mediation effects, and our simulations show that it performs well in that setting. Although the HDMT protects the nominal false discovery rate, it is much more conservative than other methods, including when $\tau_{1}$ is large. The csmGmm generally protects the nominal false discovery rate well even in more extreme settings and is less conservative than the HDMT.