Statistical Inference for Qualitative Interactions with Applications to Precision Medicine and Differential Network Analysis

As we begin to formalize the problem, we first introduce some notation. We consider two sub-populations, labeled by $g\in\{1,2\}$. Let $\theta_{g}\in\mathds{R}$ denote a measure of association in sub-population $g$. When convenient, we write $\theta=(\theta_{1},\theta_{2})$. We can consider various measures of association, such as: correlation coefficients, indicating the strength of linear relationship between two variables of interest; log odds ratios, describing the relationship between predictors and a binary outcome; and log hazard ratios, describing the association between predictors and a time-to-event outcome.

We assume that given i.i.d. samples of size $n_{1}$ and $n_{2}$ from each sub-population, $\sqrt{n_{g}}$-consistent and asymptotically normal estimates $\hat{\theta}_{g}$ of $\theta_{g}$ are available, i.e.,

with ${\sigma^{2}_{g}}>0$ denoting the asymptotic variance. For expositional simplicity, we assume balanced sample sizes $n_{1}=n_{2}=n$ (key results are stated more generally in the Appendix). We also assume ${\sigma^{2}_{g}}$ is known, though we can instead use a consistent estimate, as is commonly done in practice.

We now formally state the null hypotheses of no positive/negative interactions and no absence/presence interaction, labeled $H_{0}^{\text{P/N}}$ and $H_{0}^{\text{A/P}}$, respectively:

We let $\Theta_{0}^{\text{P/N}},\Theta_{0}^{\text{A/P}}$ and $\Theta_{1}^{\text{P/N}},\Theta_{1}^{\text{A/P}}$ denote the corresponding null and alternative regions of the parameter space, depicted in Figure 1. (Recall that the null region is the set of parameters such that the null hypothesis holds, and the alternative region is the complement of the null region.) The positive/negative null region is the union of the the non-negative and non-positive orthants, and the absence/presence null region is the union of all open orthants and the origin.

Our goal is to use the estimate $\hat{\theta}$ to perform tests of $H_{0}^{\text{P/N}}$ and $H_{0}^{\text{A/P}}$ such that the size is controlled asymptotically under mild assumptions. Recall that for a null hypothesis $H_{0}$ with accompanying null region $\Theta_{0}$, the size of a test is defined as

In words, the size is the largest possible type-I error rate that could be achieved under any probability distribution given that $\theta$ belongs to the null region.

Here, we describe our general approach for controlling the size at a pre-specified level $\alpha\in(0,1)$. We first define a test statistic $T$, a map from observable data to a real-valued number, with larger values of $T$ corresponding to more evidence against the null hypothesis. We then calculate the test statistic on the observed data, which we denote by $t$. We write $\mathbb{P}(T>t|\theta=\theta_{0})$ as the probability of observing a random test statistic at least as large as the observed value $t$, assuming $\theta=\theta_{0}$. We reject the null hypothesis if

One can think of $\mathbb{P}(T>t|\theta=\theta_{0})$ as the p-value under a specific null distribution $\theta=\theta_{0}$; $\rho(t)$ is then the largest of all such p-values. We reject $H_{0}$ when there is sufficient evidence to reject all hypotheses $\theta=\theta_{0}$. We can view $\rho(t)$ as a generalization of the usual p-value for simple null hypotheses to tests with composite null hypotheses, and will simply refer to $\rho(t)$ as “p-value”. Tests of the above form are guaranteed to control the size (Casella and Berger, 2002).

We now review existing approaches to test for qualitative interactions. We first discuss testing positive/negative interactions before moving to absence/presence interactions.

Gail and Simon (1985) developed the most widely used procedure to test for positive/negative interactions. Though Gail and Simon proposed a general $K$-sample test, we focus on the two-sample problem in this paper. We note that various $K$-sample tests for positive/negative interactions have been proposed (Piantadosi and Gail, 1993; Silvapulle, 2001; Li and Chan, 2006), but these procedures are essentially equivalent to the Gail-Simon test in the two-sample setting.

Gail and Simon’s approach is to perform a likelihood ratio test based on the asymptotic sampling distribution of $\hat{\theta}$. The likelihood ratio test rejects $H_{0}^{\text{P/N}}$ for large values of

where $\phi(\cdot)$ is the standard normal density. By performing algebreic manipulations, one can show that the likelihood ratio test equivalently rejects the null for large values of

The test statistic $T^{\text{P/N}}$ can be interpreted as the shortest distance between $\hat{\theta}$ and the null region, where the distance is inversely weighted by the asymptotic variances of the estimates, as illustrated in Figure 2.

Gail and Simon show that the test statistic $T^{\text{P/N}}$ can be calculated as

where $\mathds{1}(\cdot)$ denotes the indicator function. Furthermore, one can verify that for all $t>0$, the p-value is easily calculated as

The likelihood ratio test is quite intuitive and rejects the null when a positive estimate of association is observed in one population, a negative estimate is observed in another population, and both associations are statistically significant.

Now, we discuss approaches to test for absence/presence interactions. While one might be tempted to perform a likelihood ratio test for absence/presence interactions, the likelihood ratio test fails in the sense that it never can reject the null. To see this, we first recognize that, similar to the positive/negative interaction test, the absence/presence likelihood ratio test would reject for large values of

Again, the test statistic $T^{\text{A/P}}$ is the shortest distance between $\hat{\theta}$ and the null region $\Theta_{0}^{\text{A/P}}$. Because the alternative region $\Theta_{1}^{\text{A/P}}$ has zero area, $\hat{\theta}$ lies in the null region with probability one. Therefore, the test statistic is always 0, and the likelihood ratio test has no power.

One might alternatively attempt to test the absence/presence null by separately testing the null hypotheses $H_{0}^{g}:\theta_{g}=0$ for $g\in\{1,2\}$ and rejecting the absence/presence null when $H_{0}^{g}$ is rejected for one $g$ and not rejected for the other. To control the size of a test of this form, we need to simultaneously control the type-I error under two scenarios: (1) there is an association in both sub-populations, and (2) there is no association in either population. When there is an association in both populations, a type-I error occurs when we incorrectly fail to reject one of $H_{0}^{g}$. When there is no association in either population, we make a type-I error when we incorrectly reject one of $H_{0}^{g}$. Thus, controlling the size of the test for $H_{0}^{\text{A/P}}$ using this approach requires simultaneous control of the type-I error rate and type-II error rate of tests for $H_{0}^{g}$. If the tests for $H_{0}^{g}$ are consistent — that is, the type-II error rates tend to zero with sufficiently large samples — this approach is asymptotically valid, as only the type-I error rates for tests of $H_{0}^{g}$ need to be controlled. However, with even moderately large samples, we will not be able to correctly reject false $H_{0}^{g}$ with absolute certainty unless the true association is strong and hence easy to detect. This would make the test of $H_{0}^{\text{A/P}}$ unreliable in the presence of weak signal.

Specifically, we require $\theta_{1},\theta_{2}>o(n^{-1/2})$. To see this, we construct a more formal argument. For simplicity, suppose $\sigma_{1}=\sigma_{2}=1$. We consider tests of $H_{0}^{g}$ of the form

where $a$ is a constant that would be selected to control the size. The probability of rejecting the absence/presence null is

Suppose $\theta_{1}$ and $\theta_{2}$ are guaranteed to be greater than $o(n^{-1/2})$ if they are both nonzero. Then, for $\theta_{1},\theta_{2}\neq 0$, $\sqrt{n}|\hat{\theta}_{g}|\to\infty$, and $\mathbb{P}\left(\text{Reject }H_{0}^{\text{A/P}}\right)\to 0$. Therefore, to control the size, we are only required to select $\alpha$ so that the type-I error is controlled when $\theta_{1}=\theta_{2}=0$; this can be done by taking $a$ as the $(1-\alpha/4)$ quantile of the standard normal distribution. However, if we allow $\theta_{g}<o(n^{1/2})$, we can see a drastically inflated type-I error rate. For instance, for a small $\epsilon>0$, let $\theta_{1}=n^{-1/2+\epsilon}$, $\theta_{2}=n^{-1/2-\epsilon}$. Then $\sqrt{n}|\hat{\theta}_{1}|\to\infty>a$ while $\sqrt{n}|\hat{\theta}_{2}|\to 0<a$, so $\mathbb{P}\left(\text{Reject }H_{0}^{\text{A/P}}\right)\to 1$. Thus, when small signal is permitted, tests of this form will be asymptotically anti-conservative.

Both approaches discussed above for testing absence/presence interactions fail for a similar reason: it is difficult to gather evidence supporting that a measure of association is exactly equal to zero. This is captured by the alternative region having zero area, causing the failure of the first approach. In the second approach, to obtain evidence supporting that an association is zero, we require that $H_{0,g}$ is only accepted when $\theta_{g}=0$; for this, we rely upon a minimum signal strength condition to guarantee that any non-zero association is detected.

To mitigate the challenges described in Section 2, we consider a refinement of the absence/presence null hypothesis. The key idea is that in practice, absence/presence interactions can be approximated by considering the settings where an association is at least moderately large in one population and negligible or near zero in the other; or when one association is substantially stronger than the other. This means that we can expand the alternative region to include neighborhoods of zero in a way that the absence/presence interpretation is preserved.

Recall that when there exists an absence/presence interaction, the ratio of the maximum of the absolute value of the $\theta_{g}$ to the minimum is infinite. We cannot test that the ratio is infinite because we will never have evidence to support that the denominator is exactly zero. However, we can test that the ratio is large because we may have evidence to support that the denominator is very small. Motivated by this intuition, we propose to test whether the relative difference between sub-population measures of association is greater than a large pre-specified constant $\kappa>1$. Formally, let $\theta_{\max}=\max_{g}|\theta_{g}|$ and $\theta_{\min}=\min_{g}|\theta_{g}|$. We define the new relative difference null hypothesis $H_{0}^{\kappa}$ as

Equivalently,

The null region $\Theta_{0}^{\kappa}$, illustrated in Figure 1, is the union of four linear subspaces — each residing in a separate orthant of $\mathds{R}^{2}$; the boundary of each subspace is the union of the spans of vectors with absolute direction $(\kappa,1)^{\prime}$ and $(1,\kappa)^{\prime}$.

The relative difference null region can be viewed as a relaxation of the absence/presence null. For a large choice of $\kappa$, both our original and refined null hypotheses have the same interpretation: the greater measure of association is substantially larger than the lesser. However, we find it appealing that $H_{0}^{\kappa}$ has a reasonable interpretation for any choice of $\kappa$; that is, the multiplicative difference in strength of association is no larger than $\kappa$.

To motivate defining the refined null hypothesis in terms of relative differences rather than absolute differences, we argue that testing for relative differences has at least the following benefits. First, relative differences are unitless, so the relative difference null is compatible with unitless measures of association such as the Pearson correlation coefficient, which are often preferred in the analysis of biological data. Second, a reasonable $\kappa$ can be selected without prior knowledge of ranges of strength of association.

Of course, the relative difference null hypothesis depends on the choice of $\kappa$. For small values of $\kappa$, the relative difference null may be too dissimilar from the absence/presence null to retain its interpretation; for large $\kappa$, the alternative region becomes very small, and the hypothesis may be overly conservative. In the following subsections, we first construct a test of the relative difference null hypothesis for a pre-specified $\kappa$ and then describe an approach to identify the set of $\kappa$ such that the test rejects the null hypothesis, as to circumvent tuning parameter selection.

We now develop a testing procedure for the new relative difference null hypothesis. The relative difference null region, unlike the absence/presence null region, has non-zero area, so a likelihood ratio test will not fail in the same manner as the likelihood ratio test for absence/presence interactions. Similar to the previously discussed examples, the likelihood ratio test statistic is

and can be interpreted as the shortest (weighted) distance between $\hat{\theta}$ and the null region $\Theta_{0}^{\kappa}$. Clearly, the test statistic is zero whenever $\hat{\theta}$ lies in the null region. Otherwise, $T^{\kappa}$ is the shortest distance between $\hat{\theta}$ and the closest of the four linear subspaces that define $\Theta_{0}^{\kappa}$. The test statistic can be calculated as the distance between $(\hat{\theta}_{\max},\hat{\theta}_{\min})^{\prime}$ and its projection onto the span of the vector $(\kappa,1)^{\prime}$. The test statistic’s geometric interpretation is illustrated in Figure 2.

The likelihood ratio test statistic is straightforward to calculate. Let $\hat{\theta}_{\max}=\max_{g}|\hat{\theta}_{g}|$ and $\hat{\theta}_{\min}=\min_{g}|\hat{\theta}_{g}|$ be the strongest and weakest estimated absolute association. In the following lemma, we state that $T^{\kappa}$ is equal to the difference between $\hat{\theta}_{\max}$ and $\kappa\hat{\theta}_{\min}$ divided by a normalizing constant. Thus, the test statistic can be viewed as a plug-in of $\hat{\theta}$ into (3) with an additional normalizing constant and closely resembles the test statistic proposed by Fieller (1940) to conduct inference about ratios of means.

We now discuss how to obtain a p-value for the relative difference hypothesis. First, we obtain an observed test statistic $t$, a realization of $T^{\kappa}$ calculated from the data. Following the approach described in Section 2.2, we define the p-value as $\rho^{\kappa}(t)$, the largest of all asymptotic tail probabilities $\lim_{n\to\infty}\mathbb{P}(T^{\kappa}>t|\theta=\theta_{0})$ such that $\theta_{0}$ belongs to the null region. To determine the maximum tail probability, we characterize the limiting distribution of $T^{\kappa}$ assuming $\theta=\theta_{0}$ for all $\theta_{0}$ in the null region.

Though the null region contains an infinite number of values, $T^{\kappa}$ can only attain one of three limiting distributions corresponding to the following three cases:

1. 1.

   The true association is in the interior of the null region, i.e., $\theta_{\max}-\kappa\theta_{\min}<0$.

2. 2.

   The true association is on the boundary of the null region, but both associations are non-zero, i.e., $\theta_{\max}=\kappa\theta_{\min}>0$.

3. 3.

   The true association is zero in both sub-populations, i.e., $\theta_{1}=\theta_{2}=0$.

In Proposition 1, we describe the asymptotic behavior of $T^{\kappa}$ for cases 1 and 2 above. We provide here some intuition for the result and reserve a formal argument for the Appendix In case 1, it is easy to argue that because $\hat{\theta}$ is consistent, $T^{\kappa}$ is a negative number with probability tending to one, and therefore never provides evidence against the null. In case 2, $T^{\kappa}$ asymptotically follows a standard normal distribution. To see this, we note that because both associations are non-zero, consistency and asymptotic normality of $\hat{\theta}$ imply that, for large $n$, the sign and order of the estimates are deterministic. That the signs are asymptotically deterministic implies that $|\hat{\theta}|$ is asymptotically normal (speaking loosely, $|\hat{\theta}_{g}|\to\text{sign}(\theta_{g})\hat{\theta}_{g}$), and that order is asymptotically deterministic implies that $\hat{\theta}_{\max}$ and $\hat{\theta}_{\min}$ are asymptotically independent. Therefore, taking the difference between $\hat{\theta}_{\max}$ and $\hat{\theta}_{\min}$, suitably standardized, is asymptotically equivalent to taking a difference between two independent normal random variables with equal means. Dividing by the asymptotic variances gives the claimed result.

In case 3, when both associations are zero, the asymptotic distribution of the test statistic is more complicated. More specifically, $\theta_{g}=0$ implies that, asymptotically, $|\hat{\theta}_{g}|$ follows a half-normal distribution instead of a normal distribution. Moreover, the order of $|\hat{\theta}|$ remains random in the limit. This gives rise to a non-standard limiting distribution. In particular, unlike cases 1 and 2, the limiting distribution depends on the asymptotic variances of the sub-population estimates (and also the ratio of the sample sizes in the unbalanced case). We are nonetheless able to derive an analytic expression for the distribution function, stated in Proposition 2. We reserve this statement for the Appendix, as the expression is cumbersome.

By Propositions 1 and 2, we can calculate the p-value as the maximum of the tail probabilities in cases 2 and 3. That is,

where $\Phi(\cdot)$ denotes the standard normal distribution function and $F^{\kappa}(\cdot)\equiv\lim_{n\to\infty}\mathbb{P}(T^{\kappa}<t|\theta=\theta_{0})$ is the limiting distribution function for the test statistic when $\theta=0$. Though the limiting distribution under $\theta=0$ is non-standard, tail probabilities can be calculated easily, as we show in the Appendix (Remark 1). The p-value is therefore simple to calculate.

We have derived an analytic approximation for the power of the likelihood ratio test. In Proposition 2, we more generally characterize the limiting distribution of the test statistic under hypotheses of the form $(\theta_{1},\theta_{2})=n^{-1/2}(c_{1},c_{2})$. The local asymptotic power of the proposed test at the level $\alpha$ is available by considering $\beta^{\kappa}_{\alpha}(c_{1},c_{2})\equiv\lim_{n\to\infty}\mathbb{P}\left(T^{\kappa}>t^{*}_{1-\alpha}|\theta=n^{-1/2}(c_{1},c_{2})\right)$, where we define $t^{*}_{1-\alpha}$ as the maximum of the $1-\alpha$ quantiles of the limiting distribution of $T^{\kappa}$ under scenarios 2 and 3 described above. An analytic finite-sample approximation of the power can be calculated as $\beta^{\kappa}_{\alpha}\left(n^{1/2}\theta_{1},n^{1/2}\theta_{2}\right)$.

A contour plot of the local asymptotic power is given in Figure 3. We consider both cases of equal and unequal variance of the estimators. We observe that the likelihood ratio test has low power when the strongest effect $\max\left\{|c_{1}|,|c_{2}|\right\}$ is small, and power improves considerably when the strongest effect grows. Additionally, we find that in the presence of unequal variance, the test has greater power when the weakest effect is estimated with higher precision than the strongest effect.

Rather than perform a the likelihood ratio test for a pre-specified $\kappa$, one may prefer to directly estimate the relative difference in effect. Naïvely, one might consider estimating the relative difference in effect as $\hat{\theta}_{\max}/\hat{\theta}_{min}$. However, this estimate will behave poorly when $\theta_{1}=\theta_{2}=0$ and should therefore not be reported in practice.

To overcome this issue, we propose to quantify the relative difference in effect size by inverting the likelihood ratio test, similar to Fieller (1940). We define

as the largest $\kappa>1$ such that the likelihood ratio test rejects the null hypothesis at the $\alpha$ level. When the likelihood ratio test fails to reject for all $\kappa>1$, we will use the convention $\kappa_{\max}^{\alpha}=1$. We find it appealing that $\kappa^{\alpha}_{\max}$ converges to $1$ if $\theta_{\max}=\theta_{\min}$, and $\kappa_{\max}$ should approach but not exceed $\theta_{\max}/\theta_{\min}$ otherwise. We discuss calculation of $\kappa^{\alpha}_{\max}$ in the Appendix.

When it is of interest to identify both absence/presence and positive/negative qualitative interactions, it may be desirable to test for both simultaneously. In this section, we construct an omnibus test that achieves exact control of the size asymptotically.

We define the omnibus qualitative interaction null hypothesis as

The null region $\Theta_{0}^{\text{P/N},\kappa}$ is the intersection of the positive/negative and relative difference null regions, as depicted in Figure 4. We observe that as $\kappa\to\infty$, the omnibus null and alternative regions tend to the positive/negative null and alternative regions.

To construct the likelihood ratio test, we proceed using similar arguments to those presented in Section 3.2. The likelihood ratio statistic $T^{\text{P/N},\kappa}$ is the distance between the estimate $\hat{\theta}$ and its projection onto the null region $\Theta_{0}^{\text{P/N},\kappa}$, inversely weighted by the asymptotic variance of $\hat{\theta}$. A simple expression for $T^{\text{P/N},\kappa}$ is given in Lemma 2.

Unsurprisingly, the likelihood ratio statistic for the omnibus test approaches the likelihood ratio statistic for the Gail-Simon likelihood ratio statistic for positive/negative interactions in the limit of large $\kappa$. The tests will, therefore, be nearly identical for sufficiently large $\kappa$.

To characterize the asymptotic behavior of the omnibus test statistic at each location under the null, we use similar arguments to those in the previous section. If $\theta$ belongs to the interior of the null region, $T^{\text{P/N},\kappa}$ converges in probability to zero. If $\theta$ belongs to the boundary of the null region and is non-zero, $T^{\text{P/N},\kappa}$ converges weakly to a uniform mixture of zero and the chi-squared distribution with one degree of freedom. If $\theta$ is zero, the limiting distribution of $T^{\text{P/N},\kappa}$ is non-standard, though it can be characterized nonetheless. Formal statements of asymptotic properties of $T^{\text{P/N},\kappa}$ are given in Propositions 3 and Remark 2 (the statement of Remark 2 is also cumbersome, and is reserved for the Appendix).

Calculating the p-value for the omnibus test is no more difficult than calculating the p-value for the absence/presence test. Defining $F^{\text{P/N},\kappa}(t)\equiv\lim_{n\to\infty}\mathbb{P}(T^{\text{P/N}}<t|\theta=0)$, the p-value can be calculated as

where $t$ is the value of the test statistic $T^{\text{P/N},\kappa}$ calculated on the observed data. We characterize the local asymptotic power of the omnibus likelihood ratio test in Proposition 4 in the Appendix.

Our objective is to determine whether there are any pairs of genes that are much more strongly associated in one estrogen receptor group than the other. We consider the set of $p=145$ genes in the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000) breast cancer pathway and measure the association between gene expression levels using the Pearson correlation.

We test the relative difference null hypothesis for each pair of genes with $\kappa=2$. In Figure 7, we display the pairs of genes that are statistically significant at the $\alpha=.10$ level after a Bonferroni adjustment. We find that each of the genes progesterone (PGR), insulin-like growth factor 1 (IGF1R), and estrogen receptor 1 (ESR1) have multiple differential connections; each belongs to at least two pairs such that the association is twice as strong in the ER+ population than in the ER- population. These genes have been shown in the literature to be associated with sub-type and prognosis (Farabaugh et al., 2015; Reinert et al., 2019; Kurozumi et al., 2017).

The goal of this analysis is to assess whether any of the KEGG genes have a stronger association with time to death in one estrogen receptor group than in the other. For each gene, we fit a univariate Cox proportional-hazards model with time to death as the outcome in both of the estrogen receptor groups separately; we measure association using the log hazard ratio. A total of 64 deaths occurred in the ER+ group, and 33 deaths occurred the ER- group. We calculate $\kappa^{\alpha}_{\max}$ with $\alpha=.10$ for each gene.

In Figure 7, we compare the log hazard ratios of the ER+ and ER- groups in a scatterplot. Though the log hazard ratios for most genes are similar between subgroups, there are twelve genes with $\kappa^{\alpha}_{\max}$ larger than one. A complete list is available in Table 1. The two genes with the strongest interactions are Growth Factor Receptor-bound Protein 2 (GRB2; $\kappa^{\alpha}_{\max}=2.04$), which has a stronger association in the ER- group, and Adenomatous Polyposis Coli (APC; $\kappa^{\alpha}_{\max}=1.91$), which has a stronger association in the ER+ group. Both genes have been hypothesized to be associated with breast cancer carcinogenesis (Daly et al., 1994; Jin et al., 2001).