Let w ₂ = y ₁ y ₂ be the pairwise Hadamard product of y ₁ and y ₂. Then, if x _j represents an “intercept” column of the design matrix with all elements unity [i.e., if x _j = (1,…,1)′], we set s x j , w 2 = 0 . Otherwise, we proceed as follows: We first approximate the conditional mean and variance of y _k given x _j = x for k ∈ {1, 2} through a linear regression model as g ( y k | x ) = a k j + b k j x (5) and h ( y k | x ) = ( n − 1 ) ( s y k 2 − b k j s x j , y k ) / ( n − 2 ) (6) where b k j = s x j , y k / s x j 2 and a k j = y ̄ k − b k j x ̄ j . We note that this conditional variance will be constant at any value of x following from the linear regression assumption of homoscedasticity.

Then, we calculate the sample partial correlation of y ₁ and y ₂ controlling for x _j : r y 1 , y 2 . x j = r y 1 , y 2 − r x j , y 1 r x j , y 2 ( 1 − r x j , y 1 2 ) ( 1 − r x j , y 2 2 ) , (7) setting r y 1 , y 2 . x j = 0 if either r x j , y 1 or r x j , y 2 = 1 . As the expectation of the conditional correlation equals the partial correlation under the assumption of a multivariate linear relationship between ( y ₁, y ₂) and x _j (Baba et al., 2004), we use the partial correlation as an estimate of the conditional correlation of y ₁ and y ₂ at all possible values of x _j . So, we approximate the covariance of y ₁ and y ₂ conditional on x _j : h ( y 1 , y 2 | x ) = r y 1 , y 2 . x j h ( y 1 | x ) h ( y 2 | x ) (8)

These terms let us approximate the conditional mean of w ₂ at a given value x of x _j : g ( w 2 | x ) = g ( y 1 | x ) g ( y 2 | x ) + h ( y 1 , y 2 | x ) (9)

Then, letting f _j (x) be an assumed probability distribution/mass function for x _j with support S j [e.g., if x _j is a vector of minor allele counts with MAF p, letting f j ( x ) = ( 2 x ) p x ( 1 − p ) 2 − x and S j = { 0,1,2 } ) we approximate the sample covariance of x _j and w ₂: s x j , w 2 ≈ ∑ x ∈ S j f j ( x ) ( x − x ̄ j ) g ( w 2 | x ) , (10) swapping the sums for integrals across the support when appropriate.

We calculate the sample mean of w ₂ as w ̄ 2 = y ̄ 1 y ̄ 2 + s y 1 , y 2 ( n − 1 ) / n (11)

To approximate the variance, we first approximate the conditional variances of w ₂ at all levels of x _j : h ( w 2 | x ) = h ( y 1 | x ) h ( y 2 | x ) + g ( y 1 | x ) h ( y 2 | x ) + g ( y 2 | x ) h ( y 1 | x ) + g ( w 2 | x ) h ( y 1 , y 2 | x ) (12) And then approximate the sample variance as: s w 2 2 ≈ ∑ x ∈ S j n f j ( x ) − 1 h ( w 2 | x ) + n f j ( x ) g ( w 2 | x ) − w ̄ 2 2 / ( n − 1 ) (13) once again swapping the sum for an integral across S j when appropriate. This approach leads to a different variance estimate for each predictor x _j . We treat the median of these estimates across each j as the estimated variance.

Hence, taking the means, variances, and pairwise covariances of x _j , y ₁, and y ₂ and a distributional assumption about x _j , we approximate the covariance of one variable ( x _j ) with the product of the other two ( w ₂ = y ₁ y ₂) as well as the product’s mean and variance.

Repeating this algorithm for each predictor x _j and following the linear regression equations presented in Section 2.2 allows for calculation of covariate adjusted slope coefficients for the multiple regression model w 2 = X β ̂ + ϵ as well as the standard errors of these slope estimates.

Regression models for larger products of phenotypes can also be approximated by applying the established method recursively: first estimating the covariance of x _j and w ₂, then leveraging the covariance of x _j and w ₂ and x _j and y ₃ to estimate the covariance of x _j and w ₃, and so forth. This recursion procedure is described in more detail in the appendix and software to carry it out is discussed in Section 2.8.

Let w _l = y ₁ y ₂⋯ y _l = w _l−1 y _l . In order to estimate s x j , w l through our established method, we use x ̄ j , y ̄ l , w ̄ l − 1 , s x j 2 , s y l 2 , s w l − 1 2 , s x j , y l , s x j , w l − 1 , and s w l − 1 , y l as inputs to the method described in Section 2.3.1. That is, replacing y ₁ with w _l−1 and y ₂ with y _l . While x ̄ j , y ̄ l , s x j 2 , s y l 2 , and s x j , y l are assumed to be known, we must estimate s x j , w l − 1 and s w l − 1 , y l .

Continuation of the recursive process starting at w _l−1 and working down to w ₂ will yield an estimate for s x j , w l − 1 , or eventually the base case of s x j , w 2 .

To approximate s w l − 1 , y l , we re-express the term as s y l , w l − 2 y l − 1 and approximate the covariance of y _l and w _l−2 y _l−1 through the method described in Section 2.3.1.

A diagram of the start of this recursion is displayed in Supplementary Figure S1.

This recursive estimation is impacted by the order in which the phenotypes are multiplied. So, any set of more than two phenotypes will render m!/2 possible ways to estimate the regression model through this method. Hence, we approximate the covariances and means using all permutations of y ₁, …, y _m unique up to the order of the first two terms as the order of our phenotypes, and take the median estimate of each term across all permutations as its predicted value.

Our proposed method leverages several assumptions of the standard linear regression model in its approximations (namely those of homoscedasticity and linearity). As binary phenotypes notoriously violate both of these assumptions, we slightly adjust our approximation to better reflect this situation.

The covariance of two binary phenotypes is estimated using the same general framework as developed in Section 2.3.1. The only changes are to mean the variance estimates. All estimated means (conditional and otherwise) are initially estimated as proposed in Section 2.3.1 and then restricted to the open interval (0, 1) (e.g., letting g(y _k |x) = a _kj + b _kj x if 0 < a _kj + b _kj x < 1, ϵ if a _kj + b _kj x ≤ 0, and 1 − ϵ if 1 ≤ a _kj + b _kj  x for some small ϵ > 0).

Instead of estimating a phenotype’s conditional variance from a linear model’s residual variance, we estimate it as h ( y k | x ) = g ( y k | x ) ( 1 − g ( y k | x ) ) (14)

Further, we calculate the product’s sample variance as s w 2 2 = w ̄ 2 ( 1 − w ̄ 2 ) n / ( n − 1 ) (15)

Binary phenotypes are of particular importance because their products can be interpreted as logical combinations.

We can represent the logical conjunction y ₁ ∧ y ₂ (read as “ y ₁ and y ₂”) as the product y ₁ y ₂. Likewise, we express the logical disjunction y ₁ ∨ y ₂ (“ y ₁ or y ₂”) as 1 _n − ((1 _n − y ₁)(1 _n − y ₂)).

By framing both disjunctions and conjunctions in terms of phenotype multiplication, we can apply our established methods to approximate the covariances of these combinations with predictors and ultimately estimate linear models for these logical combinations.

While the case of the conjunction is a trivial application of the above methods of multiplying phenotypes, we will briefly describe how to model the disjunction. To do so, we consider the modified phenotypes y 1 ′ = 1 n − y 1 and y 2 ′ = 1 n − y 2 (these represent the statements “not y ₁” and “not y ₂”). This gives us y ₁ ∨  y ₂ = 1 _n − y ₁′ y ₂′. Then, y ̄ l ′ = 1 − y ̄ l , s x j , y l ′ = − s x j , y l , and s y k ′ , y l ′ = s y k , y l . If we set w 2 ′ = y ̄ 1 ′ y ̄ 2 ′ , our method allow us to estimate s x j , w 2 ′ for each x _j as well as w ̄ 2 ′ and s w ′ 2 . Leveraging these estimates, s x j , w 2 = − s x j , w 2 ′ , w ̄ 2 = 1 − w ̄ 2 ′ , and s w 2 2 = s w 2 ′ 2 , where w ₂ is equivalent to the disjunction y ₁ ∨  y ₂. Using these terms as inputs for the framework presented in Section 2.2 allow for coefficient and standard error estimation for the linear model y 1 ∨ y 2 = X β ̂ + ϵ.

To verify that our linear model with PCSS approach appropriately maintained the Type I error rate at a variety of α thresholds, we carried out a simulation under the null hypothesis that the predictor variant has no linear association with any of the phenotypes of interest. This null hypothesis represents a reasonable subset of the exact null hypothesis which is that the product of phenotypes has no linear relationship with the predictor. We carried out this simulation with varying sample size, MAF, phenotype means, phenotype correlations, and for continuous phenotypes, phenotype variances, for products of two binary phenotypes, two continuous phenotypes, and three continuous phenotypes. When simulating continuous phenotypes we assumed that Y i k ∼ N ( μ k , σ k 2 ) with Cor(Y _ik , Y _il ) = ρ _kl while when simulating binary phenotypes we let Y _ik ∼ Bernoulli(μ _k ) (again with Cor(Y _ik , Y _il ) = ρ _kl ). Simulation parameters (e.g., n, μ ₁, ρ ₁₂) were randomly sampled from various distributions (full details are available in Supplementary Table S1). We carried out 10⁸ simulations for each collection of continuous phenotypes and 10⁷ simulations for the case of binary phenotypes.

To evaluate our method’s ability to replicate the results of covariate adjusted linear models fit to IPD, we carried out three 2 ^k factorial simulations—one for the product of two binary phenotypes, one for the product of two positive continuous phenotypes, and one for the product of three positive continuous phenotypes. We carried out 1,000 simulations at each possible combination of parameters. In each simulation, we modeled the phenotype product as a function of a SNP and binary covariate. For the simulations with only two phenotypes, we also included a continuous covariate in our models.

In all simulations, we simulated n subjects’ SNP minor allele counts x ₁ at HWE with varying MAF. We simulated a binary covariate x ₂ ∼ Bernoulli(logit⁻¹α₂ x ₁). When generating sets of two phenotypes we also generated a continuous covariate x ₃ from a linear model with x ₁ such that x ̄ 3 = 0 , s x 3 2 = 1 , and r x 1 , x 3 = α 3 . This resulted in a SNP with two covariates (p = 3) in our two phenotype simulations, and a SNP with one covariate (p = 2) in our three phenotype simulation.

We generated individual phenotype measures through the model u ( y i k ) = β k 0 + ∑ j = 1 p x i j β k j + ϵ i k where u(y _ik ) = y _ik for continuous phenotypes, u(y _ik ) = logit(Pr(Y _ik = 1)) for binary phenotypes, and ϵ_i ^′ follows a multivariate normal distribution with μ = 0 and Σ _(i,j) = σ _i σ _j ρ _ij . Parameter values were selected such that, under optimal settings, empirical power was roughly 80–90% at a significance threshold of 10⁻⁸. Full details of simulation parameters are available in Supplementary Table S2.

In each simulation, we estimated coefficients, standard errors, t statistics, and two-sided p-values for the null hypothesis that there was no relationship between the product of phenotypes and the SNP ( x ₁) after adjusting for covariates both using IPD and using PCSS.

Additionally, when simulating two binary phenotypes we fit covariate-adjusted logistic regression models for the logged odds that y _1i y _2i = 1 using IPD and returned the relevant two-sided p-value to compare the results of the linear model fit using PCSS to the correctly specified logistic model.

To analyze the results of our Type I Error simulations we calculated the empirical Type I Error rate when approximating linear models using PCSS at each specified significance threshold.

For each of the three 2 ^k factorial simulations, we assessed the PCSS model’s accuracy relative to its IPD counterpart.

We calculated the bias and mean squared error when estimating the SNP’s slope coefficient, standard error, and absolute value test statistic. We also modeled errors estimating the slope coefficient, standard error, and test statistic through multiple linear regression models with logical indicators for each of the k parameter settings as predictors, testing at the Bonferroni adjusted significance threshold of 0.05/k to determine which simulation parameters affected our method’s accuracy.

We compared test decisions regarding the significance of the SNP when modeling the phenotype product after adjusting for covariates at significance thresholds 10⁻¹, 10⁻², …, 10⁻⁸. When analyzing binary phenotypes we also compared test decisions between the linear model fit using PCSS and the logistic regression model fit on IPD to demonstrate the robustness of linear models to model binary outcomes.

We measured our overall bias and mean squared error in slope, standard error, and absolute value test statistic estimates for each outcome modeled. We recorded test decisions for both the IPD and PCSS models and whether the two results agreed or disagreed regarding the statistical significance of each SNP. When one approach found a SNP to be significant and the other did not, we noted if the non-significant result was “borderline” significant (α ≤ p < 10α).

Empirical Type I error rates when using PCSS are displayed in Table 1. In all simulations, the approach’s empirical Type I error rate was below the tested significance threshold.

The PCSS method’s errors when approximating slope coefficients, their standard errors, and test statistics are available in Table 2. When aggregated over all simulation settings we observe slight positive bias both when estimating the slope and the absolute value of the t test statistic for each collection of phenotypes. The magnitude of the mean test statistic error is comparable across all three simulations. Figure 2 displays our PCSS method’s approximated slope coefficients compared to slope coefficients calculated using IPD for the SNP while modeling the phenotype product and adjusting for covariates. Similar graphical comparisons of standard error and test statistic estimates are available in Supplementary Figures S2, S3.

When modeling estimation errors for two continuous phenotypes through a linear regression model with indicator variables for all of the simulation settings (k = 12, n = 2 ^k × 10³), our model for the slope error found all settings except the residual phenotype variances, σ k 2 , to be significantly associated with the PCSS model’s slope estimate’s error at the adjusted significance threshold 0.05/k. All settings had significant associations with our error when estimating the standard error of the slope coefficient, or the test statistic. In the case of two binary phenotypes (k = 14, n = 2 ^k × 10³), we found all settings to have significant associations with the error in slope, standard error, and test statistic estimates. For three continuous phenotypes (k = 13, n = 2 ^k × 10³), we also found all settings to have significant associations with the error when predicting the slope coefficient, its standard error, and its test statistic.

Figure 3 shows comparisons of estimated and calculated p-values for a two-sided t test under the null hypothesis that the SNP had no linear association with the phenotype product after adjusting for covariates. (Figure 3 only includes p-values greater than 10⁻¹⁵ for the sake of visual clarity; Supplementary Figure S4 repeats this visualization without any restrictions.) Figure 4 shows various error rates rate between the IPD and PCSS models’ test decisions based on these p-values at differing significance thresholds. We see that all PCSS models overall disagreement rates to their IPD companions decrease as the significance threshold becomes more stringent. Likewise, when the IPD model rejected the null hypothesis, the PCSS model rarely failed to reject with error rates at most 13% which again decreased as the significance threshold decreased. When the IPD model failed to reject the null hypothesis, the PCSS approach’s conditional error rate varied by the model’s response. When modeling the product of two continuous or binary phenotypes, the error rate stayed relatively constant across all thresholds at around 3 and 15%, respectively. But, when modeling the product of three continuous phenotypes, the error rate increased as the significance threshold became more strict. Lastly, when compared to the test decisions of a covariate adjusted logistic regression model, our PCSS approximation of the related linear model tends to reach the same conclusions, with a moderate conservative tendency, especially at more strict significance thresholds.