High-throughput mediation analysis of human proteome and metabolome identifies mediators of post-bariatric surgical diabetes control

To improve the power of mediation in high-throughput studies, here we introduce High-throughput mediation analysis (Hitman), which accounts for direction of mediation and applies empirical Bayesian linear modeling. We apply Hitman in a retrospective, exploratory analysis of the SLIMM-T2D clinical trial in which participants with type 2 diabetes were randomized to Roux-en-Y gastric bypass (RYGB) or nonsurgical diabetes/weight management, and fasting plasma proteome and metabolome were assayed up to 3 years. RYGB caused greater improvement in HbA1c, which was mediated by growth hormone receptor (GHR). GHR’s mediation is more significant than clinical mediators, including BMI. GHR decreases at 3 months postoperatively alongside increased insulin-like growth factor binding proteins IGFBP1/BP2; plasma GH increased at 1 year. Experimental validation indicates (1) hepatic GHR expression decreases in post-bariatric rats; (2) GHR knockdown in primary hepatocytes decreases gluconeogenic gene expression and glucose production. Thus, RYGB may induce resistance to diabetogenic effects of GH signaling. Trial Registration: Clinicaltrials.gov NCT01073020.


Model
As discussed in the manuscript, we need to know a priori whether or not our treatment or exposure (E) changes the outcome (Y ) and in what direction. We then test whether or not the observed causal effect of the exposure on the outcome in our study is significant and is in the direction known a priori. We represent this potential causal effect as E −→ Y. Lack of an arrow indicates no causal effect, 1 so the model assumes that Y does not affect E (as occurs when E is randomized).
The test of E −→ Y possibly includes covariates, but does not include the mediators. Our unmediated model (without covariates or error terms) is E λ1 −→ Y , where the structural parameter λ 1 represents the total effect (http://davidakenny.net/cm/mediate.htm). The corresponding structural equation model (SEM), with optional covariate X. If there are multiple covariates per sample, X is a vector per sample.
The direction of the observed total effect is the sign ofλ 1 , which is the estimate of λ 1 . sgn(λ 1 ) is predefined, and we should have sgn(λ 1 ) = sgn(λ 1 ). We test the null hypothesis λ 1 = 0, and if it is significant and sgn(λ 1 ) = sgn(λ 1 ), we infer that E causally affected Y. Once we have a causal effect, we can move on to Hitman to identify which analytes may be mediating it. We are then interested in the mediation effect (or indirect effect), represented by α 1 θ 2 . θ 1 is the direct effect per mediator. As lack of an arrow indicates no causal effect, this model assumes that neither M nor Y affect E (as occurs when E is randomized), and that Y does not affect M. Hitman does not test these assumptions.
Our model in panel a is based on a biological system, such as a cell, where an intervention (or exposure) like a gene knockout can have effects across all activities in the cell. The number of potential mediators of the exposure is P, where P is astronomically large. In an omics experiment, we measure p (with p much smaller than P, p << P) analytes across n samples, with n << p. The P mediators in Figure 1a have unknown relationships among themselves and on the outcome. Ideally, we would identify these relationships so that we can account for them, but with few samples and an astronomical number of mostly unmeasured mediators, we do not attempt that here, and instead pursue an exploratory, one-at-a-time analysis of our measured mediators, so that we can identify which mediators to pursue with further studies.
Our linear SEM per mediator with covariate X is: where U terms represent omitted factors that explain finite sources of variation. The SEM's coefficients represent finite structural parameters estimable from the data. These terms must be mutually independent for the validity of mediation analyses. When we randomize E, we know that U E is independent of the other U 's. However, U M is often dependent with U Y , due to another mediator (likely unmeasured) affecting both M and Y . So, "theoretical knowledge must be invoked to identify the sources of these correlations and control for common causes (so called 'confounders') of M and Y whenever they are observable". 2 Such confounders can be included in X. When X is a vector per sample, ν 1 , α 2 , and θ 3 are vectors.
Hitman is based on the causal steps approach. One variant of the causal steps approach is the joint significance test. Other variants of the causal steps approach test that the exposure affects the outcome, but the joint significance test does not 34 . The joint significance test has been shown to have more power than the product method and to control its false positive rate, because it is an intersection-union test. 5 The joint significance test calls as significant "consistent" mediators, which carry the causal effect, but it also calls as significant "inconsistent" mediators, which suppress the causal effect. However, the causal effect of the exposure on the outcome provides prior evidence and motivates the search only for consistent mediators.  4 We address this in Hitman.
For each mediator (M ), the null hypothesis is that E has no effect on M (α 1 = 0); or that M has no effect on Y (θ 2 = 0); or that the direction of mediation, sgn(α 1 θ 2 ), is not consistent with the direction of the total effect, sgn(λ 1 ), with sgn being the sign function, The alternative hypothesis is that E affects M and M affects Y , and that the direction of E −→ M and M −→ Y are consistent with that of the total effect of E on Y . This tests if M explains at least some of the dependence between E and Y ; M need not explain all of the dependence, but it could, in which case we have complete mediation and the direct effect vanishes.
Hitman implements tests with high-throughput mediators using the R/Bioconductor linear modeling package Limma. 6 Limma models the variance of features (e.g. proteins or metabolites) with an empirical Bayesian method, which exploits information about shared technical variance between features for improved power, especially when the sample size is small. 6 To model variance of feature abundance in linear modeling, Limma models feature abundance as the dependent variable.

Workflow Assumptions
A1. E is randomized (given X, if X is provided) A2. M and Y follow the linear SEM in equations 2 and 3 A3.

Before Hitman
We need to ensure that the assumptions hold. Herein we describe and validate Hitman/Lotman for the simple case where the SEM is linear and M and Y being normally distributed, but it is straightforward to extend our approach to E being randomized given propensity scores; the value of M being associated with its variance, as happens with RNA-seq data, which Limma accommodates with weights; 7 and Y being modeled by a generalized linear model, such as Y following a binomial distribution.
The assumption that most differentiates Hitman from other mediation methods is A5. To ensure A5 holds, we pre-define sgn (λ 1 ) and then we test the total effect λ 1 with H 0 : λ 1 = 0 against H a : λ 1 ̸ = 0. This test must be significant and the observed causal effect must be in the same direction as previously known. If the test is not significant, we cannot be confident in the direction of the total effect in our study. Examples where A5 would not hold is if surgery did not cause a significant change in HbA1c or if surgery worsened HbA1c, which would contradict prior knowledge. If this pre-test fails, then STOP, because Hitman would assign all mediators a p-value of one, which is not likely to be helpful.

Hitman
For each mediator M : 1. Measure effect of E and M using equation 2. Here and later, if X is absent, remove its term. If X is vector-valued per sample, α 2 is a vector. Test α 1 = 0 in Limma and define the resulting p-value p 1 .
2. Measure association of M and Y given E (and possibly X) using equation 3. We would like to test θ 2 = 0 using Limma. To model the variance of M with empirical Bayesian methods, we need to make M the dependent variable. So we use an approach similar to partial correlation.
i. Estimate the residuals of ii. Estimate the residuals of M = α 0 + α 1 E + α 2 X + e as e M .
iii. From the linear model e M =θ 2 e Y + ϵ, test θ 2 = 0 in Limma and define the resulting p-value p 2 . . Otherwise, the direction of the indirect effect is inconsistent with that of the total effect, and the final p-value = 1. The intuition behind the test is that knowing λ 1 allows for one-sided testing of the indirect effect.

Now define
To account for testing multiple mediators, false discovery rates (FDRs) or family-wise error rates (FWERs) can be calculated from the mediator p-values.

Simplified Hitman without Limma: Lotman
If we follow the Hitman workflow but apply linear regression models without Limma, in step 2 we test θ 2 = 0 from equation 3 directly. This method could be used when Limma is not appropriate, such as for low-throughput data, so we term it Low-throughput mediation analysis (Lotman), and it is available in the R package Hitman.

Mathematical proof
Here, we consider Hitman theoretically to show that it controls its false positive rate. Hitman's use of the partial correlation approach in step 2 is a valid alternative to multiple regression, and incorporating Limma maintains control of the false positive rate, because Limma is theoretically sound 8 and empirically validated. 6 We do not consider these aspects theoretically here, although Hitman controls its false positive rate in simulations below, providing evidence that these aspects do not corrupt Hitman's validity. Hitman without these aspects is Lotman, so we consider the validity of both methods by essentially considering Lotman, but we point out where results with Limma would differ. Lotman applies tests in steps 1 and 2 with ordinary least squares (OLS), which provides unbiased estimates with valid p-values.
For a single mediator, the null and alternative hypotheses are: where ∩ represents "and" and ∪ represents "or". The terms in H a are not independent, because the definition of S depends on α 1 , θ 2 , and λ 1 in equation 4. We have rejected the null hypothesis that λ 1 = 0 before applying Hitman, soλ 1 ̸ = 0, and the results of this test are convincing, so we assume that there is a true causal effect in the direction observed, i.e. sgn (λ 1 ) = sgn λ 1 ̸ = 0. Without loss of generality, we consider Thus, S = 0 =⇒ H 0 . Both α 1 =⇒ S = 0 and θ 2 =⇒ S = 0, so H 0 =⇒ S = 0: From equation 4, S = 0 ∩ sgn (λ 1 ) = 1 =⇒ sgn (α 1 θ 2 ) ̸ = 1. We can see from the Hitman algorithm and hypotheses that α 1 and θ 2 are treated symmetrically, so without loss of generality, we consider the null hypothesis space, where by assumption A5 we also haveλ 1 > 0. To demonstrate that Hitman yields valid p-values, we need to show that Hitman p-values, p HM , obey P (p HM ≤ u|H ⋆ 0 ,λ 1 > 0) ≤ u for each 0 ≤ u ≤ 1. Hitman is partially based on the joint significance test, and it has been shown that the joint significance p-values (p JS ) satisfy this property. 5 where 0 ≤ p JS 2 ≤ 1 2 . We can immediately see for u = 1 that P (p HM ≤ 1) ≤ 1, and for u = 0, P (p HM = 0) = 0, because p HM is continuous over 0, 1 2 . We now consider how p 1 and p 2 are calculated. We know by assumption A4 that 0 < σ α1 , σ θ2 < ∞. For arbitrary sample size N (sufficient to estimate parameters and variances) and full rank covariate matrix X of rank k, using OLS we have , where := refers to "defined by," and T 2 : Whereas if we consider the asymptotic case or the actual instead of the estimated standard deviation, T 1 and T 2 follow the standard normal distribution with density ϕ and cumulative distribution function Φ. The degrees of freedom of the t-statistics will be increased with Limma's empirical Bayesian variance estimate, especially when N − k is small.
Instead of relying on the t-statistics converging to the normal distribution asymptotically, our mathematical proof extends to small sample sizes, so we account for the two t-distributions' different degrees of freedom.
We define the T 1 as having cumulative distribution function F 1 and probability density function f 1 , whereas T 2 as having cumulative distribution function F 2 and probability density function f 2 . To explain properties that apply to t-distributions (i.e. independent of their degrees of freedom, so including the standard normal distribution as a special case), we use F or f with subscripts suppressed. For example, we use below the properties of the t-distribution that We also suppress subscripts for the joint distribution ofα 1 andθ 2 .

Overview of mathematical proof
We have shown above that P (p HM ≤ u|H ⋆ 0 ,λ 1 > 0) ≤ u for u = 0 and for u = 1, so we must show this for 0 < u < 1. Equation 7 says that we need to consider P (Ŝ = 1|H ⋆ 0 ,λ 1 > 0), so our next step is to This bound allows us to immediately show that P (p HM ≤ u|H ⋆ 0 ,λ 1 > 0) ≤ u for u ≥ 1 2 . However, the case of u < 1 2 is more challenging, so we tackle it in several steps. First we show that P (p HM ≤ u|H ⋆ 0 ,λ 1 > 0) ≤ u when either α 1 = 0 or θ 2 = 0, which is on the boundary of the space We then evaluate the false positive rate over the full space α 1 × θ 2 | H ⋆ 0 ,λ 1 > 0. To find the optima of this space, we calculate its derivative and identify where the derivative is zero. We then calculate the values on the optima and on the boundary. We find that the maxima of P (p HM ≤ u|H ⋆ 0 ,λ 1 > 0) is on the boundary, where we know P (p HM ≤ u|H ⋆ 0 ,λ 1 > 0) ≤ u, which completes the proof.

Maximizing probability estimated direction is consistent
Under H ⋆ 0 we also haveλ 1 > 0. However we suppress the notation of conditioning on these variables for brevity in the following equation. We also haveα 1 ⊥ ⊥θ 2 . 9 Then, 1 2 To find the upper bound of equation 9, we maximize equation 9 as which cannot be solved by convex programming. However, we can examine where one of v 1 or v 2 is fixed, giving a linear program. Consider the linear program for any fixed v 1 ∈ ( 1 2 , 1]: , and under the above constraints, max Case of u ≥ 1 2 A straightforward consequence of the bound on P (Ŝ = 1) is for 1 2 ≤ u < 1. Under H ⋆ 0 andλ 1 > 0, we see that

False positive rate on the boundary
We next consider the special case of α 1 × θ 2 on its boundary under the null. The boundary of H ⋆ 0 for α 1 × θ 2 is α 1 = 0 ∪ θ 2 = 0. A property of the boundary is that it corresponds to the null hypothesis of the joint significance test, where P (p JS ≤ u) ≤ u. We can use this fact to show p HM control its false positive rate under H ⋆ 0 on the boundary for 0 < u < 1 2 : In the special case α 1 = 0 and θ 2 → ∞, p JS = p 1 so: because p 1 ∼ U (0, 1). Analogously, in the special case α 1 → −∞ and θ 2 = 0, p JS = p 2 , and we find that P (p HM ≤ u|α 1 → −∞, θ 2 = 0) = u. So we see in these special cases that P (p HM ≤ u) reaches its maximum value of u.
When both parameters are zero,

False positive rate in general
As above, under H ⋆ 0 we haveλ 1 > 0. We now define a density and cumulative distribution function ofα 1 that suppresses notation: whereσα 1 ⊥ ⊥ α 1 ,α 1 . We have similar notation forθ 2 with density and distribution having subscript 2. We determine the false positive rate for fixed 0 < u < 1 2 through the law of total probability, where for fixed values a, t and u, the probabilities P (a ≤ F −1 1 (u)) and P (t ≥ −F −1 2 (u)) resolve to zero or one and are independent because a and t are realizations of random variables, rather than being random variables themselves. That is, We continue equation 14 to evaluate P (p HM ≤ u|H ⋆ 0 ,λ 1 > 0,σα 1 ,σθ 2 ) but suppress conditioning on H ⋆ 0 ,λ 1 > 0,σα 1 ,σθ 2 : When one of the parameters is on the boundary at zero, say α 1 = 0, then equation 14 simplifies to
So off the boundary there is the optimum, = lim α1→−∞,θ2→∞ which is a minimum. Thus, the false positive rate is controlled off the boundary and on the boundary, so the false positive rate is controlled in general.

Example of inconsistent mediator
We provide a clear example of an inconsistent mediator in the figure below, which includes regression lines and 95% confidence interval bands. Source data is available at src_data_hm_toy.csv. This figure is reminiscent of Simpson's Paradox. 10 The exposure increases the mediator and the exposure increases the outcome, so the direct effect is positive. But the mediator decreases the outcome, so the indirect effect is negative. Thus, we have inconsistent mediation, so Hitman conservatively estimates this mediator's p-value as 1, whereas the joint significance method would yield a p-value < 0.001 for this example.

Simulations following Barfield et al. (2017) & MacKinnon et al. (2002)
We validated our size and power following a simulation study. 11 This study's parameters were a subset of a previous simulation study (MacKinnon et al., 2002). Unlike the methods tested with these simulation scenarios, Hitman requires a predefined direction of the total effect, and includes a pre-test that the total effect is significant and that its direction agrees with the observed total effect direction. It would seem unfair to provide Hitman the true direction of the total effect and this also wouldn't be sufficient to increase the chances that the total effect is significant. Instead, we could address this with a robust direct effect, so that the observed effect will tend to be significant and in the same direction as the true total effect. Although Barfield et al. (2017) set the direct effect to be "small", θ 1 = 0.14, MacKinnon et al. (2002) included a scenario with a parameter that was "large" with value 0.59, so we chose to set θ 1 = 0.59.
To account for the high-throughput data Hitman is designed to be applied to, we simulated other, null mediators as M i = α 0 + α 2 X + e Mi for i = 2, 3, ..., 10, which are independent of the exposure and the outcome. These null mediators are utilized by Limma to estimate the shared technical variance, as would happen in an omics study. We do not test these other mediators.
We tested in what proportion of 10,000 simulations across 50 samples M achieved a p-value ≤ 0.05. Our results for these parameter values for Hitman, the joint significance test, and the mediate function from the R package mediation are shown in the table below. A statistical test's size is the probability of falsely rejecting the null hypothesis, which is the probability of a false positive or a Type 1 error. All the methods here control their size, as they maintain a false positive rate less than 5%.
A statistical test's power is the probability that the test correctly rejects the null hypothesis when the alternative hypothesis is true, and it's inversely related to the probability of making a Type 2 error. To test Hitman's power, we look at the cases where both of θ 2 and α 1 are positive. Here, Hitman's power is greater than the other methods.
Hitman's power is similar to Lotman because we have a relatively large sample size. However, for a smaller sample size, we can see a larger difference between Hitman and Lotman. Here we show the same

Omics simulations
We also validated our methods in an omics setting with 500 analytes (which can interchangably be called "features" or "genes"; g = 1, 2, ..., G with G = 500), where we address the multiplicity issue by controlling the false discovery rate (FDR). We simulated the data under 3 scenarios. In all scenarios, we simulated datasets with sample sizes of N=15 or N=50 and either 1, 5, or 25 consistent mediators. The exposure often affects many analytes, so simulated that the exposure additionally affects 200 analytes with the same effect size that it affects mediators. However, these 200 analytes are not simulated to affect the outcome, so they are not mediators.
In scenario 3, for each consistent mediator there was also an inconsistent mediator of the same effect size. As above, we set θ 0 = θ 3 = α 0 = α 2 = 0.14, and here we set β 1 = θ 2 = 2. To satisfy Hitman's requirement that the total effect be known a priori and observed in the dataset where Hitman is applied, we wanted θ 1 >> 2 * number of mediators, so we set θ 1 = 10 * number of mediators (including both consistent and inconsistent mediators).
We tested each gene for mediation using Hitman, Lotman, and the joint significance method. We corrected the p-values using the Benjamini-Hochberg method, thresholded the FDR at 15%, and calculated the number of potential mediators whose null is rejected, the proportion of true mediators whose null is rejected (power), and the proportion of false rejections among all rejections, which is the false discovery proportion (FDP). We conducted 1,000 such simulations and report the averages (means) of these statistics. Because FDR=E(FDP), the mean FDP is an estimate of the empirical FDR.
We still simulated following equation (1) here, with E ∼ N (0, 1), X ∼ N (0, 1), and ν 0 = ν 1 = 0. The matrix M with 500 rows and N columns is simulated as where α 0 , α 1 , α 2 ∈ R G ; α 1g = 1 if gene g is a consistent or inconsistent mediator or is otherwise simulated to be affected by the exposure and α 1g = 0 otherwise; and U M ∼ N G (0, Σ) with Σ a G-by-G positive definite matrix.
Our covariance matrix is derived from GTEx v8 RNA-seq gene read counts downloaded from https:// gtexportal.org/home/datasets. We choose to focus on skeletal muscle tissue, since it had the most samples. We filtered out genes with low expression, applied TMM normalization, 12 accounted for measured covariates known to affect expression and inferred technical factors with SVA, like, 13 calculated the covariance matrix of 500 randomly sampled genes, and then scaled this matrix so that the median variance would be one. The script implementing this process is gtex_muscle_cov_subset_scaled.R and the processed matrix is at gtex_muscle_cov_subset_scaled.csv.
When genes are independent (scenario 1), Σ is the matrix whose diagonal elements are the same as the diagonal of the scaled covariance matrix, and all off-diagonal elements are zero. When genes are dependent (scenarios 2 & 3) Σ is the scaled covariance matrix.
Finally, N (0, 1) and if gene g is an inconsistent mediator 0, if gene g is not a mediator 1, if gene g is a consistent mediator.
For scenario 3, the joint significance method has power to identify both consistent and inconsistent mediators, whereas Hitman and Lotman only have power for consistent mediators, so we had a choice as to assess the joint significance method's power and FDR with respect to only consistent mediators (like Hitman and Lotman) or with respect to both types of mediators. We chose to show both types of assessments: one as "Joint signif consistent" and the other as "Joint signif both".
The omics simulation results in the figure below shows the FDR and power of the considered methods in the three scenarios. Our code implementing the simulations is at simulations_omics.Rmd and the CSV tables with the simulation output are in the folder validation whose file names include "sc1", "sc2", or "sc3" for scenarios 1, 2, and 3, respectively.